A modified numerical manifold method for simulation of finite deformation problem

Accepted Manuscript

A modified numerical manifold method for simulation of finite deformation problem Wei Wei , Qinghui Jiang PII: DOI: Reference:

S0307-904X(17)30298-6 10.1016/j.apm.2017.04.026 APM 11734

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

16 October 2016 25 March 2017 25 April 2017

Please cite this article as: Wei Wei , Qinghui Jiang , A modified numerical manifold method for simulation of finite deformation problem, Applied Mathematical Modelling (2017), doi: 10.1016/j.apm.2017.04.026

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Highlights  False solutions have been found when modeling finite deformation problems using the original NMM. A new numerical manifold method for simulation of finite deformation is proposed.



The proposed method eliminates the errors caused by large rotation and large strain in the original NMM.

AC

CE

PT

ED

M

AN US

CR IP T



ACCEPTED MANUSCRIPT

A modified numerical manifold method for simulation of finite deformation problem

Wei Wei 1

Qinghui Jiang 1, 2*

1. School of Civil Engineering, Wuhan University, Wuhan 430072, P. R. China

CR IP T

2. School of Civil Engineering and Architecture, Nanchang University, Nanchang 330031, P. R. China

AC

CE

PT

ED

M

AN US

*Corresponding author: E-mail: [email protected] Telephone: +86-27-68772221

Abstract: With many people contributing to its modifications and advancements, the numerical manifold method (NMM) is now recognized as an efficient tool to solve the continuum-discontinuum coupling problem in geotechnical engineering. However, false solutions have been found when modeling finite deformation problems using the

ACCEPTED MANUSCRIPT

original NMM. Based on the finite deformation theory, a modified version of NMM is derived from the weak form of conservation of momentum and the corresponding traction boundary condition. By taking the dual cover system as the displacement approximation, the governing equations of the modified NMM are formulated. A comparison of the governing equations of the original NMM and modified NMM illustrates the reason that the original NMM is not suitable for simulation of finite

CR IP T

deformation problems. Three numerical examples are investigated to verify the capability of proposed method to predict static and dynamic finite deformation response. Numerical results show that the modified NMM eliminates the errors caused by large rotation and large strain, and obtains a good agreement with analytical

AN US

solutions and the finite element method.

AC

CE

PT

ED

M

Keywords: Numerical manifold method; finite deformation; large strain; large rotation

ACCEPTED MANUSCRIPT

1 Introduction The numerical manifold method (NMM), which was initially proposed by Shi [1, 2], is a combination of the finite element method (FEM) and discontinuous deformation analysis (DDA) [3]. It provides a unified way to solve continuous and

CR IP T

discontinuous problems in geotechnical engineering. Over the past few decades, numerous modifications and enhancements [4] has been done to improve this method. Owing to its advantage in solving continuous and discontinuous deformations, most NMM researchers concentrated on crack propagation problems. Tsay et al. [5]

AN US

firstly employed NMM to investigate crack growth path using local remeshing technique. Chiou et al. [6] applied the NMM to predict mixed mode crack growth combined with the virtual crack extension model. Li et al. [7] developed the meshless manifold method to solve two-dimension crack problems. Terada et al. [8] extend the

M

finite cover method (FCM), an alias of NMM, to analyze the progressive failure processes using the cohesive fracture zone model in heterogeneous solids and

ED

structures. Kurumatani and Terada [9] extended the FCM to simulate the crack growth of quasi-brittle heterogeneous solids by introducing „multi-cover-layer modeling‟. Ma

PT

et al. [10] adopted NMM to model the complex crack problems such as multiple branched and intersecting cracks. Zheng and Xu [11] discussed the rank deficiency,

CE

the integrals with singularity of 1/r, kinked cracks and the mesh independency of NMM in the simulation of crack growth using NMM. Zheng and Liu [12] reduced the

AC

problem of growth of multiple cracks to a nonlinear complementarity problem (NCP) and solved it by the projection-contraction algorithm PCA Using the fracturing algorithm based on the Mohr–Coulomb criterion with a

tensile cutoff, Ning et al. [13], An et al. [14] and Wong [15] made preliminary studies to model progressive failure in rock slopes and analyzed the stability of slopes. By introducing a rock bolt element in the NMM, Wu and Wong [16] investigated the rock block falling case in reinforced tunnels. Wei et al. [17] established a numerical model

ACCEPTED MANUSCRIPT

for modelling the bolts across the joint surfaces in the NMM. Zhang et al. [18] proposed a method combining NMM and graph theory to assess the slope stability and obtain the safety factor and the critical slip surface of slopes. By developing a sequential excavation algorithm, Tal et al. [19] applied the improved NMM to analyze the stability of tunnels in the Jingping hydropower project in Sichuan Province, China. Except for the applications to structures and solids, the NMM is also employed to

CR IP T

conduct seepage analysis [20-22], the analysis of thin plate bending problems [23] and so on.

To improve the accuracy of solutions, high order approximations of field variants in the NMM are constructed [24, 25], which might lead to the issue of linear

AN US

dependency (LD). An et al. [26] dissected the origin of LD and proposed an approach to predict the rank deficiency of the global stiffness matrix. Ghasemzadeh et al. [27] developed dynamic high order NMM without linear dependence problems. Recently, Tian and Wu [28] proposed an effective high order scheme avoiding the linear

M

dependency problem in the XFEM setting, which is also applicable to the NMM. So far the NMM has been successfully applied to many aspects. However, the

ED

simulation of finite deformation problem using NMM is seldom reported. Fan and Zheng [29] proposed a new S-R-D Based NMM to simulate small or large

PT

deformation together with impact/contact. In this present study, a modified version of NMM based on finite deformation theory is proposed. By considering the

CE

conservation of mass, the derivation of discrete equation of the modified NMM started from the weak form of conservation of momentum and the corresponding

AC

traction boundary condition. Then the two cover systems were introduced to construct the approximation of displacement field. To illustrate the reason that the original NMM for calculating the finite deformation problems is inaccurate, the comparison of governing equations of the original and modified NMM is demonstrated. Finally, three numerical examples are investigated to verify the accuracy and validity of proposed method.

ACCEPTED MANUSCRIPT

2 Governing equations As shown in Fig.1, we consider a body which occupies a domain Ω0 with a boundary Γ0 in initial state at a time t=0. The domain Ω0 is called the initial configuration and the position of a material point in this configuration is denoted by X.

CR IP T

when forces act on the body, it starts movement and deformation. The body occupied a domain with a boundary in a certain time t. Similarly, the domain Ω is called the current configuration and the position of a material point in this configuration is denote by x. Therefore, the motion of the body is denoted by 𝐱 = 𝐮(𝐗, 𝑡) + 𝐗

AN US

where, 𝐮 is the displacement of the material point 𝐗.

(1)

The governing equations for the body is given as follows:

(2)

𝛔∇ = 𝐂 ∇ ∶ 𝐃

(3)

𝐃 = 𝑠𝑦𝑚(𝛁 ∙ 𝐮̇ )

(4)

M

𝛁 ∙ 𝛔 + 𝜌𝐛 = 𝜌𝐮̈

here, 𝛁 is the gradient operator, 𝛔 is the Cauchy stress, 𝜌 is the density of material,

ED

𝐛 is the body acceleration vector. 𝐮, 𝐮̇ , 𝐮̈ denote the displacement, velocity and acceleration vector, respectively. 𝛔∇ represents any objective rate of the Cauchy stress.

PT

𝐂 ∇ is the tangential modulus corresponding to 𝛔∇ , here the superscript “∇” refers to an objective rate. D is the rate of deformation.

CE

A compressible Neo-Hookean type is selected to model the constitutive behavior

AC

and the strain function can be expressed as 1

𝜓(𝐼, 𝐽) = 2 𝜆0 (ln𝐽)2 − 𝜇0 (𝐼 − 3 − 2ln𝐽)

(5)

where I is the first invariants of the deformation, J is the determinant of deformation gradient. 𝜆0 and 𝜇0 are the lame constants in the initial configuration. In the current configuration, the tangential modulus 𝐂 ∇T corresponding to the Truesdell rate of Cauchy stress 𝛔∇T in two-dimensional plain stain is written as 𝐂

∇T

=𝐽

−1

𝜆 + 2𝜇 [ 𝜆 0

𝜆 𝜆 + 2𝜇 0

0 0] 𝜇

(6)

ACCEPTED MANUSCRIPT

with 𝜆 = 𝜆0 , 𝜇 = 𝜇0 − 𝜆ln𝐽

(7)

As shown in Eq. (6), the tangential modulus 𝐂 ∇T has a similar form as in Hooke‟s law for small strain problem. The displacement boundary condition on Γ𝑢 and the traction boundary condition on Γ𝑡 are given as follow: on

Γ𝑢

𝛔𝐧 = 𝐭̅

on

Γ𝑡

(8)

CR IP T

𝐮=𝐮 ̅

(9)

The initial condition for the problem is provided as follow: 𝐮(𝐗, 𝑡0 ) = 𝐮0 (𝐗)

𝐮̇ (𝐗, 𝑡0 ) = 𝐮̇ 0 (𝐗)

(10)

the principle of mass conservation

AN US

In addition, during the process of deformation and movement, the body should satisfy

𝜌(𝐗, 𝑡)𝐽(𝐗, 𝑡) = 𝜌0 (𝐗)

(11)

where 𝜌0 is the initial density in the initial configuration, and 𝜌 is the density in the

M

current configuration.

ED

3 Fundamental of NMM

PT

The numerical manifold method (NMM), originally proposed by G.H. Shi [1], provide a unified framework for solving continuous and discontinuous problems in

CE

geotechnical engineering. In the NMM, two cover systems, the mathematical cover (MC) and physical cover (PC) respectively, are introduced to discretize the problem

AC

domain. The MC consists of a set of overlapped small patches and should be large enough to cover the physical domain. These small patches can span discontinues and the boundary of physical domain. Each small patch in MC is termed as a mathematical patch and denoted by 𝑀𝐼 (I =1~𝑛𝑀 ), where 𝑛𝑀 is the total number of mathematical patches in MC. The PC is the union of all physical patches which are generated by intersection of mathematical patches and physical domain. That is, the external boundary and/or internal discontinuities may intersect each mathematical

ACCEPTED MANUSCRIPT

patch into several isolated pieces. Each pieces, if it is within the physical domain, is 𝑗

termed as a physical patch, denoted by 𝑃𝐼 (j = 1~ 𝑚𝐼 ), where 𝑚𝐼 is the number of physical patches generated from 𝑀𝐼 . The manifold element (ME) is the common region of several overlapped physical patches. On each mathematical patch, a weight function is defined such that 𝜑𝐼 (𝐗) ≥ 0, 𝐗 ∈ 𝑀𝐼

CR IP T

(12a)

𝜑𝐼 (𝐗) = 0, 𝐗 ∉ 𝑀𝐼

(12b)

∑ 𝜑𝐼 (𝐗) = 1, 𝐗 ∈ Ω

(12c)

where Eqs. (12a) and (12b) indicate that the weight function is non-zero values only on its corresponding mathematical patches, but zero elsewhere; while Eq. (12c) is the

AN US

partition of unity property to assure a conforming approximation. The weight function 𝑗

𝜑𝐼 (𝐗) associated with the 𝑀𝐼 will be accordingly transferred to 𝑃𝐼 , any of the PCs 𝑗

𝑗

𝑗

𝑃𝐼 in 𝑀𝐼 , which is expressed as 𝜑𝐼 (𝐗) on 𝑃𝐼 . For convenience, a single index i is

M

adopted to denote two indices of the physical patch, namely 𝐼−1

ED

𝑖(𝐼, 𝑗) = ∑ 𝑚𝑙 + 𝑗

(13)

𝑙=1

where, 𝑚𝑙 is the number of the physical patches formed in mathematical patch l.

PT

On each physical patch 𝑃𝑖 , a local approximation 𝐮𝑖 (𝐗, 𝑡) reflecting the local characteristic of the field is defined. For two-dimensional problem, a convenient way

CE

to construct a basis of local approximation spaces is through using polynomial

AC

functions, namely,

𝑆 = *1, 𝑋, 𝑌, … , 𝑋 𝑝 , 𝑋 𝑝−1 𝑌, … , 𝑋𝑌 𝑝−1 , 𝑌 𝑝 +

(14)

Then, the local approximation defined in physical patch i can be expressed as 𝐮𝑖 (𝐗, 𝑡) = 𝐩𝑇 (𝐗) ∙ 𝐝𝑖 (𝑡)

(15)

where 𝐝𝑖 is the vector of unknowns to be calculated, and 𝐩𝑇 is the matrix of polynomial base which may be constant, linear- or higher-order terms. The general form of 𝐩𝑇 is expressed as 𝐩𝑇 (𝐗) = [

1 0 𝑋 0 1 0

0 𝑌 𝑋 0

0 𝑌

… …

𝑌𝑝 0 ] 0 𝑌𝑝

(16)

ACCEPTED MANUSCRIPT

However, high-order terms adopted in NMM will lead to Linear Dependence Problem, which have been discussed by An et al. [26], Ghasemzadeh et al. [27], Tian and Wen [28]. For simplicity, the constant term is adopted in this paper, with 1 0 ] 0 1 𝑑 𝐝𝑖 (𝑡) = [ 𝑖1 ] 𝑑𝑖2

𝐩𝑇 (𝐗) = [

(17) (18)

CR IP T

With the above concepts, the displacement approximation of a manifold element (ME) is obtained by pasting the local approximation of physical patches with corresponding weight functions. Considering the local approximation of physical patch is constant, the displacement approximation of ME can be expressed as 𝑛

𝐮

𝑡) = ∑ 𝜑𝑖 (𝐗) ∙ 𝐮𝑖 (𝐗, 𝑡) = 𝐓𝐝𝑒

AN US

ℎ (𝐗,

(19)

𝑖=1

where

𝜑 0 𝜑2 𝐓=[ 1 0 𝜑1 0 𝑑12

𝑑21

𝑑22

…

𝑑𝑛1

𝑑𝑛2 -𝑇

(20) (21)

M

𝐝𝑒 = ,𝑑11

0 … 𝜑𝑛 0 𝜑2 … 0 𝜑𝑛 ]

ED

and n is the number of physical covers sharing manifold element ME.

PT

4 Numerical algorithm for finite deformation 4.1 The semi-discrete momentum equation

CE

In this section, a numerical algorithm for analysis of finite deformation problem

using NMM is proposed. Using the variation principle, we obtain the weak form of Eq.

AC

(2) and Eq. (9) as

∫ 𝝈 ∶ 𝛿𝐃𝑑Ω − ∫ 𝜌𝐛 ∶ 𝛿𝐮𝑑Ω − ∫ 𝐭̅ ∶ 𝛿𝐮𝑑Ω + ∫ 𝜌𝐮̈ ∶ 𝛿𝐮𝑑Ω = 0 Ω

Ω

Γ𝑡

(22)

Ω

Where 𝛿 is the first order variation. In Eq. (22), the first two terms and the last term on the left are derived from Eq. (2), while the third term is derived from Eq.(9). In the NMM, finite-dimensional subspaces, Vℎ ⊂ V and V0ℎ ⊂ V0 , are used as the approximating trial and test space. Using the Galerkin method, the weak form for

ACCEPTED MANUSCRIPT

discrete problem can be stated as: Find 𝒖ℎ ∈ Vℎ ⊂ V such that ∫ 𝝈(𝐮ℎ ) ∶ 𝐃(𝛿𝐮ℎ )𝑑Ω − ∫ 𝜌𝐛 ∶ 𝛿𝐮ℎ 𝑑Ω − ∫ 𝐭̅ ∶ 𝛿𝐮ℎ 𝑑Ω + ∫ 𝜌𝐮̈ ℎ ∶ 𝛿𝐮ℎ 𝑑Ω = 0 (23) Ω

Ω

Γ𝑡

Ω

The trial function 𝐮ℎ and the test function 𝛿𝐮ℎ are expressed as 𝐮ℎ = ∑ 𝜑𝑖 (𝐗) ∙ (𝐩𝑇 (𝐗) ∙ 𝐝𝑖 ) = 𝐓𝐝

CR IP T

𝑖

(24)

𝛿𝐮ℎ = ∑ 𝜑𝑖 (𝐗) ∙ (𝐩𝑇 (𝐗) ∙ 𝛿𝐝𝑖 ) = 𝐓𝛿𝐝 𝑖

(25)

Substituting Eq. (24) of trial function and Eq. (25) of test functions into Eq. (23) with

AN US

the arbitrariness of test function, the discrete approximation to the weak form (22) is written as

𝐌𝐝̈ℎ + 𝐟 𝑖𝑛𝑡 = 𝐟 𝑒𝑥𝑡

(26)

where M is the mass matrix, and 𝐟 𝑖𝑛𝑡 , 𝐟 𝑒𝑥𝑡 are the internal nodal force and the

M

external nodal force, respectively. The element contribution to M, 𝐟 𝑖𝑛𝑡 , 𝐟 𝑒𝑥𝑡 are as follows

(27)

𝐟𝒆𝑖𝑛𝑡 = ∫ 𝐁 𝑇 𝝈𝑑Ω

(28)

ED

𝐌𝑒 = ∫ 𝜌𝐓 𝑇 𝐓𝑑Ω = ∫ 𝜌0 𝐓 𝑇 𝐓𝑑Ω Ω𝑒0

PT

Ω𝑒

Ω𝑒

CE

̅ Ω 𝐟𝒆𝑒𝑥𝑡 = ∫ 𝜌𝐓 𝑇 𝐛𝑑Ω + ∫ 𝐓 𝑇 𝐭𝑑 Ω𝑒

(29)

Ω𝑒

AC

with

𝐁=

𝜕𝜑1 𝜕𝑥 0 𝜕𝜑1 [ 𝜕𝑦

0 𝜕𝜑1 𝜕𝑦 𝜕𝜑1 𝜕𝑥

𝜕𝜑2 𝜕𝑥 0 𝜕𝜑2 𝜕𝑦

𝜕𝜑𝑛 𝜕𝑥

0 𝜕𝜑2 𝜕𝑦 𝜕𝜑2 𝜕𝑥

⋯

0 𝜕𝜑𝑛 𝜕𝑦

0 𝜕𝜑𝑛 𝜕𝑦 𝜕𝜑𝑛 𝜕𝑥 ]

(30)

4.2 The time integration scheme Eq. (26) is a semi-discrete momentum equations since that they have not been discretized in time. Similar to the time integration of the conventional NMM, the implicit time integrator called the Newmark β-method is implemented to discretize Eq.

ACCEPTED MANUSCRIPT

(26) in time. In the time integration formula, the updated displacements and velocities are obtained by 𝐝𝑛+1 = 𝐝̃𝑛+1 + 𝛽Δ𝑡 2 𝐝̈𝑛+1 𝐝̃𝑛+1 = 𝐝𝑛 + Δ𝑡𝐝̇𝑛+1 +

(31𝑎)

Δ𝑡 2 (1 − 2𝛽)𝐝̈𝑛 2

(31𝑏)

with 𝐝̇𝑛+1 = 𝐝̃̇𝑛+1 + γΔ𝑡𝐝̈𝑛+1

CR IP T

𝐝̃̇ 𝑛+1 = 𝐝̇𝑛 + (1 − γ)Δ𝑡𝐝̈𝑛

(32𝑎) (32𝑏)

here, the parameter 𝛽 determine the accuracy of time integration. The parameter 𝛾 is an algorithm damping introduced in Newmark β-method, which can weaken the

AN US

oscillation of calculating process and converge to equilibrium state quickly. The algorithm damping is equal to zero with 𝛾 = 0.5 and nonzero with 𝛾 > 0.5. In the original NMM, the parameters 𝛽 and 𝛾 are equal to 0.5 and 1.0, respectively. Therefore, the algorithm damping is introduced in the original NMM. More details about algorithm damping refer to the work by Jiang et al. [30]. In Eqs. (31) and (32),

M

we separate the historical terms of nodal variables calculated in time step n, 𝐝̃𝑛+1 and

ED

𝐝̃̇ 𝑛+1 , which is convenient for the following algebraic operation and for the construction of implicit time integration procedures.

PT

Substituting Eq. (31) into Eq. (26) at time step n+1 gives 1 𝐌(𝐝𝑛+1 − 𝐝̃𝑛+1 ) − 𝐟 𝑒𝑥𝑡 (𝐝𝑛+1 ) + 𝐟 𝑖𝑛𝑡 (𝐝𝑛+1 ) = 𝟎 𝛽∆𝑡 2

(33)

CE

𝐫=

which is a set of nonlinear algebraic equations in the nodal displacements 𝐝𝑛+1 . The

AC

Newton-Raphson method can be adopted to solve Eq.(33), namely

Here, 𝐝𝑛+1 𝑣

𝜕𝐫(𝐝𝑛+1 𝑣 ) ∆𝐝 (34) 𝜕𝐝 is the displacement obtained by v-step iterations in time step n+1. The 𝑛+1 𝐫(𝐝𝑛+1 𝑣+1 ) ≈ 𝐫(𝐝𝑣 ) +

matrix 𝜕𝐫/𝜕𝐝 is the Jacobian matrix (or effective tangent stiffness matrix), denoted by A 𝐀(𝐝) ≈

𝜕𝐫(𝐝) 1 𝜕𝐟 𝑖𝑛𝑡 (𝐝) 𝜕𝐟 𝑒𝑥𝑡 (𝐝) = 𝐌 + − 𝜕𝐝 𝛽∆𝑡 2 𝜕𝐝 𝜕𝐝

(35)

Considering that the external forces including traction force along boundary and

ACCEPTED MANUSCRIPT

volume force are constant in geotechnical engineering, the external nodal force is independent with the deformation of the body. Therefore, the tangent stiffness matrix of the external nodal force is zero, namely 𝜕𝐟 𝑒𝑥𝑡 (𝐝) =0 𝜕𝐝 The tangent stiffness matrix of the internal nodal force is expressed as 𝐊 𝑒𝑥𝑡 =

𝜕𝐟 𝑖𝑛𝑡 (𝐝) = = 𝐊 𝑚𝑎𝑡 + 𝐊 𝑔𝑒𝑜 𝜕𝐝

(37)

𝐊 𝑚𝑎𝑡 = ∫ 𝐁 T 𝐂 𝜎𝑇 𝐁𝑑Ω

(38)

CR IP T

𝐊

𝑖𝑛𝑡

with Ω

𝑇 ,𝝈-𝐁 𝐊 𝑔𝑒𝑜 = ∫ 𝐁𝑁𝐿 ̂ 𝑁𝐿 𝑑Ω Ω

AN US

and

M

𝜎𝑥𝑥 𝜏𝑥𝑦 0 𝜏𝑦𝑥 𝜎𝑦𝑦 0 ,𝝈 ̂- = [ 𝜎𝑥𝑥 0 0 𝜏𝑦𝑥 0 0 𝜕𝜑1 𝜕𝜑2 0 0 𝜕𝑥 𝜕𝑥 𝜕𝜑1 𝜕𝜑2 0 0 𝜕𝑦 𝜕𝑦 = ⋯ 𝜕𝜑1 𝜕𝜑2 0 0 𝜕𝑥 𝜕𝑥 𝜕𝜑1 𝜕𝜑2 0 0 [ 𝜕𝑦 𝜕𝑦

PT

ED

𝐁𝑁𝐿

(36)

0 0 𝜏𝑥𝑦 ] 𝜎𝑦𝑦 𝜕𝜑𝑛 0 𝜕𝑥 𝜕𝜑𝑛 0 𝜕𝑦 𝜕𝜑𝑛 0 𝜕𝑥 𝜕𝜑𝑛 0 𝜕𝑥 ]

(39)

(40)

(41)

where 𝐊 𝑚𝑎𝑡 is the material tangent stiffness depending on the material response (or

CE

the constitutive of material). 𝐊 𝑔𝑒𝑜 is the geometric tangent stiffness involving the current state of stress.

AC

In each iteration step v, the stress increment is obtained by 𝑛+1 𝑛+1 𝑇 𝑛+1 ∆,𝝈-𝑣 = ∆,𝝈-△𝑇 𝑣 + 𝐋 ∙ ,𝝈-(𝒅𝒗 ) + ,𝝈-(𝒅𝒗 ) ∙ 𝐋 − 𝑡𝑟(𝐋),𝝈-(𝒅𝒗 )

(42)

where,

𝐋 = 𝐖(𝐓∆𝒅𝑣 )𝑇

(43)

𝜕 𝜕𝑥

𝜕 𝑇 ] 𝜕𝑦

(44)

𝜎𝑥𝑥 ,𝝈- = [𝜏

𝜏𝑥𝑦 𝜎𝑦𝑦 ]

(45)

𝐖=[

𝑦𝑥

ACCEPTED MANUSCRIPT

4.3 The algorithm for modified NMM The algorithm for the modified NMM is shown in Table 1. The calculation begins with the imposition of the initial conditions. The displacement 𝐝𝑛+1 at time step n+1 is computed by full Newton-Raphson iterative procedure, where the Jacobian matrix A is inverted in each iteration of the procedure. At the beginning of iterative procedure, an estimate value of 𝐝𝑛+1 is needed, which is usually equal to 0

CR IP T

𝐝𝑛 . And then, the residual nodal force 𝐫(𝐝𝑛+1 𝒗 ) is calculated for this starting value, which depends on the internal nodal force, external nodal force and inertia nodal force. The Jacobian matrix A in this algorithm is calculated based on the latest state of the body. Simple essential boundary conditions, such as homogeneous displacement

AN US

conditions, can be enforced by modifying the Jacobian matrix. The equation corresponding to the vanishing displacement component is either omitted or replaced by a dummy equation that the component vanishes by zeroing the corresponding row and column and putting a one on the diagonal of the Jacobian. In each time step, the

M

displacement and the stress are updated until the convergence criterion is met. The whole simulation is completed when the time step n equal to N.

original NMM

ED

4.4 Discussion on the iterative algorithm between the modified NMM and

PT

Table 2 presents the iterative algorithms of the original NMM and modified NMM in the time step n. Compared with the algorithm of the original NMM, the

CE

differences from that of modified NMM are as follows: (1) The effective tangential stiffness matrix A is comprised of the following three

AC

terms: material stiffness 𝐊 𝑚𝑎𝑡 , geometric stiffness 𝐊 𝑔𝑒𝑜 and external stiffness 𝐊 𝑒𝑥𝑡 . Each term in A is nonlinear algebraic matrix. Therefore, the displacement 𝒅𝑛+1 at time step n+1 should be computed through full Newton-Raphson iterative procedure. In the original NMM, the tangential stiffness matrix is the material stiffness 𝐊 𝑚𝑎𝑡 and assumed to be constant in each time step. While in the modified NMM, the material stiffness 𝐊 𝑚𝑎𝑡 is not a constant since that the Truesdell modulus 𝐂 ∇T is changed with the movement and deformation of the body. (2) The integration field of mass matrix in the original NMM is based on the

ACCEPTED MANUSCRIPT

current configuration. However, the density 𝜌 is equal to the initial density 𝜌0 in the simulation, which would result in considerable error when the body undergoes a large volume deformation. (3) When the body undergoes large deformation or large rotation, the rotation of element in the body will result in the rotation of principal stress axis, which implies that the components of the stresses will change in the global coordinate system. In the

CR IP T

modified NMM, as shown in Eq. (42), the incremental stress ∆𝝈𝑣 is corrected by adding three stress rotation corrections. While, in the original NMM, the stress increment is calculated by 𝝈𝑛+1 = 𝝈𝑛 + 𝐂𝐁∆𝐝

(46)

AN US

Here, the 𝐂 is the elastic constitutive. Eq. (46) indicates that the stress increment is added to 𝝈𝑛+1 without stress rotation correction. Moreover, the tangential stiffness matrix 𝐊 is constant in each time step, which will result in large error with large

ED

5 Numerical examples

M

rotation.

5.1 Static large displacement analysis of a cantilever

PT

The first example is used to verify the capability of the modified NMM for solving static small strain but large rotation problems. In this example, a cantilever

CE

illustrated in Fig.3 subject to a uniformly distributed pressure is investigated. The length and height of the cantilever is 100mm and 10mm, respectively. A

AC

Neo-Hookean elastic constitutive model is considered for the cantilever, with Young‟s modulus E=120MPa and Poisson ratio v=0.2. The cantilever is fixed on the left end and the uniform distributed load P is applied to the upper and lower surface of the cantilever. The load retains its vertical direction during the loading process for the purpose that it is independent from the deformation of the cantilever. The manifold element mesh implemented in this example is shown in Fig. 4. Fig. 5 shows displacement-load curves obtained by the analytical solution,

ACCEPTED MANUSCRIPT

original NMM and modified NMM, respectively. As can be seen, the curve predicted by the original NMM is a straight line. Whereas, the curve predicted by the modified NMM are nonlinear and show a good agreement with the analytical solution proposed by Holden [31]. The deformed configuration and contour of y-directional displacement of the cantilever with P=50kN/m is illustrated in Fig.6. The section of deformed cantilever predicted by the original NMM is expanded, which is obviously

CR IP T

unreasonable. As a comparison, the section of deformed cantilever illustrated in Fig. 6(b) predicted by the modified NMM did not expand, which is consistent to that predicted by Holden. In addition, as shown in Fig.5 and Fig.6, when P is equal to 50kN/m, the y-directional displacement of monitoring point A calculated by the

AN US

original NMM is 11.6cm, which is 4.7cm larger than that obtained by analytical solution or the modified NMM and the relative error is 68%.

5.2 Dynamic large-displacement rotation of a square block

The second example is used to verify the capability of the proposed method to

M

solve dynamic small strain but large rotation problem. As depicted in Fig. 7, a squared block, with 1cm length, rotated around the lower left corner without initial velocity. In

ED

the process of rotating of the block, the gravity force is applied to the block and the value of gravitational acceleration is 10m2⁄s. The Neo-Hookean material behavior is

PT

assumed for the squared block and the following parameters are used: Young‟s modulus E=1GPa, Possion‟s ratio v=0.2 and the density ρ=2500kg/m3. Fig.8 shows

CE

the mathematical cover used to solve the problem. The time span for numerical simulation is 1.2s and is divided into 12000 time steps. Therefore, the increment of

AC

time step is △t =10-4 s in the calculation. Different from the problem of small displacement pendulum, there is not an analytical solution to the problem of large rotating pendulum studied in this example. Consequently, as a comparison, the FEM software ABAQUS was selected to predict dynamic response of the block. The FEM mesh with 64 elements and 81 nodes is illustrated in Fig. 9. The point A on the top right of block is chosen as the monitoring point. In order to be consistence with the time integration scheme adopted by the ABAQUS software, the time integral parameters of the modified NMM are taken as

ACCEPTED MANUSCRIPT

𝛽 = 0.25 and 𝛾 = 0.5, which means the algorithm damping is zero. Fig.10 shows the comparison of dynamic response obtained by the modified NMM and FEM. As can be seen, the curves predicted by the modified NMM show a good agreement with those by the ABAQUS. The x- and y- displacements versus time curves of monitor point A are presented in Fig.11 and Fig.12. As can be seen from Figs. 11-12, the curves obtained by the

CR IP T

modified NMM are almost identical to those obtained by the original NMM in the first half cycle. However, when the block swings back to the initial position, the y-direction displacement of monitor point A is 0.09cm larger than the initial position and the cycle length of this block increases. When swinging two and half cycle, the

AN US

block is overturning and rotates clockwise around axis O. Furthermore, as shown in Fig.13, the volume of the block gets smaller gradually during rotation. When the calculation is completed, the volume of the block obtained by the original NMM is reduced to 0.83 cm2, which is 0.17 cm2 less than the initial volume of the block, and

M

the relative error is 17%. The phenomenon shown in Fig. 11 and Fig. 13 indicated that the algorithm proposed by the original NMM can not meet conservation of energy and

ED

mass for the large rotation problem. For the modified NMM, the algorithm damping is included in the time integral scheme. Therefore, energy dissipation is observed when

PT

the block is swinging, and the amplitude of the curve is gradually reduced to zero until the block returns to the final equilibrium position.

CE

5.3 Static and dynamic analysis of a rubber sheet with a hole In the third example, the static and dynamic analysis of a rubber sheet with a

AC

hole, as illustrated in Fig.14, is studied. The purpose of this example is to demonstrate the capability of predicting the static and dynamic large strain response using the modified NMM. The material behavior of this rubber sheet was considered to be of Neo-Hookean type, with the parameters E=120MPa and v=0.2. The geometry of this problem is depicted in Fig.14 with the width of the squared sheet W=1m and the radius of the hole R=0.3 m. The manifold element mesh implemented in this example is shown in Fig.15. As a comparison, the solution of FEM provided by the ABAQUS software is also

ACCEPTED MANUSCRIPT

presented. The static load-displacement curves for three typical monitoring points on the sheet are depicted in Fig.16. It is noted that the displacement of point A and B are in the x direction and point C is in the y direction. As can be seen from Fig.16, the load-displacement curves obtained by the improved NMM are nonlinear and have a good agreement with those of ABAQUS. While, the curves predicted by the original NMM is linear. Furthermore, compared with the displacements of monitoring points

CR IP T

A, B and C obtained by the modified NMM or ABAQUS, the absolute errors of displacements calculated by the original NMM are -0.044m, -0.117m and 0.045m, respectively, when the uniform distributed load q reaches 20KN/m2. Accordingly, the relative errors of displacement are 10.5%, 18.8% and 28.8%, respectively. The

AN US

deformation configuration and the contour of x-directional displacement obtained by those three methods are shown in Fig.17. It can be observed that the deformed configuration predicted by the modified NMM is identical to that by ABAQUS. For the dynamic analysis, the time span is 0.1s and has been divided into 10000

M

times steps. Namely, the time step △t is taken as 10-2 ms. Similar to the second example, the time integral parameters of the modified NMM are taken as 𝛽 = 0.25

ED

and 𝛾 = 0.5, which are consistent with those used in the FEM. Fig. 18 show the displacement response predicted by the modified NMM and ABAQUS software. It

PT

can be seen that the dynamic response predicted by the modified NMM have a good correspondence with that by the FEM method.

CE

For the comparison between the modified NMM and original NMM, the time integral parameters are taken as 𝛽 = 0.5 and 𝛾 = 1 . Fig.19 compared the

AC

displacement response predicted using the modified NMM and original NMM. Since the time step size adopted in this example is small, the decrease of amplitude caused by algorithm damping is not significant in the first phase. Therefore, the curve predicted by the original NMM is similar to that by the modified NMM in the first half cycle. However, the curve predicted by the original NMM gradually deviates from that obtained by the modified NMM in the next half cycle. As discussed in Section 4.4, the accumulative error of stresses calculated by the original NMM increases gradually with increasing time step, and the internal nodal force obtained by

ACCEPTED MANUSCRIPT

the original NMM is less than that by the modified NMM in each time step. When the sheet reaches an ultimate state, for the modified NMM, the internal nodal force is equal to the sum of the external nodal force and the inertia nodal force. However, for the original NMM, the nodal force is less than the sum of the external nodal force and the inertia nodal force. Therefore, the sheet continues to deform and the curves

CR IP T

gradually deviate from those obtained by the modified NMM.

6 Conclusion

In the present study, a modified NMM based on finite deformation theory is

AN US

proposed to solve large strain and large rotation problem. Using the Neo-Hookean constitutive and the approximation of manifold cover, the governing equations of the modified NMM is derived based on the integral weak form of the momentum conservation equation and the stress boundary conditions referring to the current

M

configuration. Compared with the governing equations of the original NMM, the stiffness matrix in the modified NMM is nonlinear in each time step and a stress

ED

rotation correction is necessary to correct the stress increment. Three numerical examples are presented to verify the accuracy and validity of

PT

the proposed method in solving finite deformation problems. For the small strain problem with large rotation, the static analysis of a cantilever beam and the dynamic

CE

analysis of a square block rotation are implemented respectively. For the large strain problem, a static and dynamic analyses of a sheet with a hole is subsequently

AC

simulated. It is shown that the results predicted by the original NMM are consistent with the analytical or other numerical solution only in the case that the structures undergo small stain with small rotation, while have a great deviation from those reference solutions when the structures undergo large rotation or large strain deformation. However, the solutions obtained by the modified NMM present a good correspondence with analytical or numerical results in both cases. With those modifications in the present study, three issues in the original NMM

ACCEPTED MANUSCRIPT

has been resolved when simulating the large deformation and rotation of structures, namely, false volume contraction, and disobeying conversation of energy and conservation of mass. Future work will be extended to improve contact formulations for finite deformation problems.

The work reported in this paper has received financial

CR IP T

Acknowledgements.

support from the National Natural Science Foundation of China (Nos. 51679173) and Natural Science Foundation of Hubei Province (No. 2016CFA083).

AN US

Reference

[1] Shi G. Manifold method of material analysis. In Transactions of the Ninth Army Conference on Applied Mathematics and Computing. DTIC Document: Minneapolis, 1992; 51–76. [2] Shi G. Modeling rock joints and blocks by manifold method." In Rock Mechanics

M

Proceedings of the 33rd US Symposium, Santa Fe, New Mexico, 1992; 639–648. [3] Shi G. Discontinuous deformation analysis: a new numerical model for the statics and

ED

dynamics of block systems. Ph.D. Thesis, University of California, Berkeley, 1988. [4] Ma GW, An XM, He L. The numerical manifold method: a review. International Journal of

PT

Computational Methods 2010; 07(1): 1-32. [5] Tsay RJ, Chiou YJ, Chuang WL. Crack growth prediction by manifold method. Journal of

CE

engineering mechanics 1999; 125(8):884-890.

[6] Chiou YJ, Lee YM, Tsay RJ. Mixed mode fracture propagation by manifold method.

AC

International Journal of Fracture 2002; 114(4): 327-347.

[7] Li SC, Cheng YM. Enriched meshless manifold method for two-dimensional crack modeling. Theoretical and applied fracture mechanics 2005; 44(3): 234-248.

[8] Terada K, Ishii T, Kyoya T, Kishino Y. Finite cover method for progressive failure with cohesive zone fracture in heterogeneous solids and structures. Computational Mechanics 2007; 39(39):191-210. [9] Kurumatani, M, Terada K. Finite cover method with multi-cover layers for the analysis of

ACCEPTED MANUSCRIPT

evolving discontinuities in heterogeneous media. International Journal for Numerical Methods in Engineering 2009; 79(1): 1-24. [10] Ma GW, An XM, Zhang HH. Modeling complex crack problems using the numerical manifold method. International Journal of Fracture 2009; 156(1):21-35. [11] Zheng H, Xu DD. New strategies for some issues of numerical manifold method in simulation of crack propagation. International Journal for Numerical Methods in

CR IP T

Engineering 2014; 97(97):986-1010.

[12] Zheng H, Liu F, Du X. Complementarity problem arising from static growth of multiple cracks and MLS-based numerical manifold method. Computer Methods in Applied Mechanics and Engineering, 2015, 295: 150-171.

AN US

[13] Ning, YJ, An XM, Ma GW. Footwall slope stability analysis with the numerical manifold method. International Journal of Rock Mechanics and Mining Sciences 2011. 48(6): 964-975.

[14] An XM, Ning YJ, Ma GW, He L. Modeling progressive failures in rock slopes with

M

non-persistent joints using the numerical manifold method. International Journal for Numerical and Analytical Methods in Geomechanics 2014; 38(7):679–701.

ED

[15] Wong LNY, Wu ZJ. Application of the numerical manifold method to model progressive failure in rock slopes. Engineering Fracture Mechanics 2014; 119:1-20.

PT

[16] Wu ZJ, Wong LNY. Underground rockfall stability analysis using the numerical manifold method. Advances in Engineering Software 2014; 76:69-85.

CE

[17] Wei W, Jiang QH, Peng J. New rock bolt model and numerical implementation in numerical manifold

method.

International

Journal

of

Geomechanics

2016;

DOI:

AC

10.1061/(ASCE)GM.1943-5622.0000669.

[18] Zheng W, Zhuang XY, Tannant DD, Cai Y, Nunoo S. Unified continuum/discontinuum modeling framework for slope stability assessment. Engineering Geology 2014; 179: 90-101.

[19] Tal Y, Hatzor YH, Feng XT. An improved numerical manifold method for simulation of sequential excavation in fractured rocks. International Journal of Rock Mechanics and Mining Sciences 2014; 65: 116-128. [20] Jiang QH, Deng SS, Zhou CB, Lu WB. Modeling unconfined seepage flow using three-dimensional numerical manifold method.

Journal of Hydrodynamics 2010;

ACCEPTED MANUSCRIPT

22(4):554-561 [21] Hu MS, Wang Y, Rutqvist J. An effective approach for modeling fluid flow in heterogeneous media using numerical manifold method. International Journal for Numerical Methods in Fluids 2015; 77 (8): 459-476. [22] Zheng H, Liu F, Li CG. Primal mixed solution to unconfined seepage flow in porous media with numerical manifold method. Applied Mathematical Modelling 2015; 39(2): 794-808.

CR IP T

[23] Zheng H , Liu ZJ, Ge XR. Numerical manifold space of Hermitian form and application to Kirchhoff's thin plate problems. International Journal for Numerical Methods in Engineering 2013; 95(9): 721-739.

[24] Chen GQ, Ohnishi Y, Ito T. Development of high-order manifold method. International

AN US

Journal for Numerical Methods in Engineering 1998; 43(4): 685-712.

[25] Zhang GX, Sugiura Y, Hasegawa H, Wang G. The second order manifold method with six node triangle mesh. Proceedings of the Japan Society of Civil Engineers 2002, 19(1), 1S-9S. [26] An XM, Li LX, Ma GW, Zhang HH. Prediction of rank deficiency in partition of unity-based

M

methods with plane triangular or quadrilateral meshes. Computer Methods in Applied Mechanics and Engineering 2011; 200(5): 665-674.

ED

[27] Ghasemzadeh H, Ramezanpour MA, Bodaghpour S. Dynamic high order numerical manifold method based on weighted residual method. International Journal for Numerical Methods in

PT

Engineering 2014; 100(8): 596-619.

[28] Tian R, Wen LF. Improved XFEM-An extra-dof free, well-conditioning, and interpolating

CE

XFEM. Computer Methods in Applied Mechanics and Engineering 2015; 285: 639-658. [29] Huo F, Zheng H, and He S. S–R decomposition based numerical manifold method. Computer

AC

Methods in Applied Mechanics and Engineering 2016;304:452-478.

[30] Jiang QH, Chen YF, Zhou CB, Yeung, MCR. Kinetic energy dissipation and convergence criterion of discontinuous deformations analysis (DDA) for geotechnical engineering. Rock mechanics and rock engineering 2013; 46(6): 1443-1460. [31] Holden JT. On the finite deflections of thin beams. International Journal of Solids and Structures 1972; 8(8): 1051-1055.

ACCEPTED MANUSCRIPT Figure list:

AC

CE

PT

ED

M

AN US

CR IP T

Figure 1. Configurations of body Figure 2. The forming of manifold element Figure 3. The geometry of cantilever Figure 4. The NMM discretization of cantilever Figure 5. The y directional displacement-load curve of monitor point A Figure 6. Contour of y directional displacement of cantilever with load P=5KN/m Figure 7. The geometry of block Figure 8. The NMM discretization of block Figure 9. The FEM discretization of block Figure 10. x and y displacement-time curve of monitor point A Figure 11. y displacement-time curve of monitor point A Figure 12. x displacement-time curve of monitor point A Figure 13. Block volume-time curve Figure 14. The geometry of a rubber sheet with a hole Figure 15. The NMM discretization of a sheet with a hole Figure 16. Static load-displacement curve for a sheet with hole Figure 17. The deformed configuration and contour of x directional displacement of the sheet Figure 18. The displacement-time curves of monitor points compared between modified NMM and ABAQUS Figure 19. The displacement-time curves of monitor points compared between modified NMM and Original NMM

ACCEPTED MANUSCRIPT

y, Y Initial configuration

u(X , t) Ω0

u Γ0

X

O

CE

PT

ED

M

Figure 1. Configurations of body

AC

Ω

AN US

x

CR IP T

Current configuration

Γ

x, X

ACCEPTED MANUSCRIPT

2

PK

MK

1

PK

2

MJ

CR IP T

PJ

1

PJ

e

MI Math cover and Physical grid

1

PI Physical cover

AC

CE

PT

ED

M

AN US

Figure 2. The forming of manifold element

Manifold element

ACCEPTED MANUSCRIPT

P/2 A

h

P/2

CR IP T

L

AC

CE

PT

ED

M

AN US

Figure 3. The geometry of cantilever

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Figure 4. The NMM discretization of cantilever

ACCEPTED MANUSCRIPT

-1

u/10 m 1.2

1.0

CR IP T

0.8

0.6 0.4

Modified NMM original NMM analytical solution

0.0

0

AN US

0.2

1

2

3

4

P/KN

5

AC

CE

PT

ED

M

Figure 5. The y directional displacement-load curve of monitor point A

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

(a) original NMM (b) modified NMM Figure 6. Contour of y directional displacement of cantilever with load P=5KN/m

ACCEPTED MANUSCRIPT

CR IP1cm T

A

AN US

O y

1cm

M

x

AC

CE

PT

ED

Figure 7. The geometry of block

ACCEPTED MANUSCRIPT

AN US

CR IP T

A

M

O

AC

CE

PT

ED

Figure 8. The NMM discretization of block

ACCEPTED MANUSCRIPT

AN US

CR IP T

A

M

O

AC

CE

PT

ED

Figure 9. The FEM discretization of block

ACCEPTED MANUSCRIPT

displacement/cm

initial position

0.5 0.0

CR IP T

-0.5 -1.0

y displ (modified NMM) x displ (modified NMM) y displ (Abaqus) x displ (Abaqus)

-1.5 -2.0

time/t

0.0

AN US

-2.5 0.2

0.4

0.6

0.8

AC

CE

PT

ED

M

Figure 10. x and y displacement-time curve of monitor point A

ACCEPTED MANUSCRIPT

The block is overturning

displacement/cm

0.5

initial position

0.0 -0.5

-1.5

CR IP T

-1.0 modified NMM original NMM

-2.0

time/t

-2.5 0.0

0.2

0.4

0.6

0.8

1.0

AC

CE

PT

ED

M

AN US

Figure 11. y displacement-time curve of monitor point A

1.2

ACCEPTED MANUSCRIPT

displacement/cm

0.5

initial position

0.0 -0.5

CR IP T

-1.0 modified NMM original NMM

-1.5 -2.0

time/t

-2.5 0.0

0.2

0.4

0.6

0.8

1.0

AC

CE

PT

ED

M

AN US

Figure 12. x displacement-time curve of monitor point A

1.2

ACCEPTED MANUSCRIPT

area/cm2 1.0

modified NMM original NMM

time/t

0.8

0.0

CR IP T

0.9

0.4

0.8

AC

CE

PT

ED

M

AN US

Figure 13. Block volume-time curve

1.2

ACCEPTED MANUSCRIPT

1

CR IP T

C

q

y

AN US

1

B Unit:m x

A

AC

CE

PT

ED

M

Figure 14. The geometry of a rubber sheet with a hole

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

Figure 15. The NMM discretization of a sheet with a hole

ACCEPTED MANUSCRIPT

q/kN·m-2 20 A (Modified NMM) B(Modified NMM) C (Modified NMM) A (Original NMM) B (Original NMM) C (Original NMM) A (ABAQUS) B (ABAQUS) C (ABAQUS)

CR IP T

15 10

AN US

5 0 0

0.2

0.4

displacement/m

0.6

AC

CE

PT

ED

M

Figure 16. Static load-displacement curve for a sheet with hole

ACCEPTED MANUSCRIPT

1

XDisp 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05

y/m

0.8 0.6 0.4

0 0

0.5

1

x/m (a) modified NMM

0.8

y/m

1.5

AN US

1

CR IP T

0.2

0.6 0.4

M

0.2

ED

0

0.5

x/m

1

1.5

(b) ABAQUS

PT

0

CE

1

AC

y/m

0.8

XDisp 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05

0.6

0.4 0.2 0

0

0.5

1

1.5

x/m (c) original NMM Figure 17. The deformed configuration and contour of x directional displacement of the sheet

ACCEPTED MANUSCRIPT

displacement/m 1.2

A (modified NMM) B (modified NMM) C (modified NMM) A (ABAQUS) B (ABAQUS) C (ABAQUS)

1 0.8 0.6

0.2 0

CR IP T

0.4

time/ms

-0.2

0

5

10

15

20

25

30

35

AC

CE

PT

ED

M

AN US

Figure 18. The displacement-time curves of monitor points compared between modified NMM and ABAQUS

ACCEPTED MANUSCRIPT

displacement/m A (modified NMM) B (modified NMM) C (modified NMM) A (Original NMM) B (Original NMM) C (Original NMM)

1.2 1 0.8 0.6

0.2 0

CR IP T

0.4

time/ms

-0.2

0

5

10

15

20

25

30

35

AN US

Figure 19 The displacement-time curves of monitor points compared between modified NMM and Original NMM

M

Table 1 the algorithm for modified NMM

PT

ED

1. initial conditions and the initialization of parameters; a. set 𝒗0 , 𝛔0; 𝒅0 = 0 b. compute the mass matrix M by Eq.(27) 2. LOOP for time step n, n = 0, 1, …, N-1 2.1 get initial nodal force 𝒇(𝒅𝑛 ) = 𝒇𝑒𝑥𝑡 (𝒅𝑛 ) − 𝒇𝑖𝑛𝑡 (𝒅𝑛 ) 2.2 compute accelerations 𝒅̈𝑛 = 𝑴−1 𝒇𝑛

AC

CE

2.3 estimate the initial displacement at n+1 time step 𝒅𝑛+1 = 𝒅𝑛 0 2.4 LOOP for Newton iterations v in at n+1 time step, v = 0,1…,convergence. a. compute the nodal force 𝒇(𝒅𝑛+1 𝒗 ) 1 ̃ 𝑛+1 ) b. compute the acceleration 𝒅̈𝑛+1 = 𝛽∆𝑡 2 (𝒅𝑛+1 −𝒅 𝒗 𝒗

̈ 𝑛+1 − 𝒇(𝒅𝑛+1 c. computer the residual 𝒓(𝒅𝑛+1 𝒗 )= 𝑴𝒅𝒗 𝒗 ) 𝑛+1 ) d. compute Jacobian 𝐀(𝒅𝑣 by Eq.(35) and be modified by essential boundary conditions e. get incremental displacement Δ𝒅𝒗 = 𝑨−1 𝒓 𝑛+1 f. 𝒅𝑛+1 + ∆𝒅𝑣 𝒗+𝟏 = 𝒅𝑣 h. update stress 𝝈(𝒅𝑛+1 𝒗+𝟏 ) by Eq.(42) 2.5 END LOOP if the convergence criterion (|∆𝒅𝑣 | < 𝑒) is meet 3. ENDLOOP if the simulation complete

ACCEPTED MANUSCRIPT

Table 2 the iterative algorithms of the original NMM and modified NMM in the time step n method

The iterative algorithm 𝐝𝑛+1 = 𝐝𝑛 0

The modified NMM

𝑛+1 𝑛+1 𝑛+1 𝐀(𝐝𝑛+1 + ∆𝐝𝑣 𝑣 )∆𝐝𝑣 = 𝐫(𝐝𝑣 ) 𝐝𝑣+1 = 𝐝𝑣 𝑛+1 𝑛+1 𝑇 𝑛+1 ∆,𝝈-𝑣 = ∆,𝝈-△𝑇 𝑣 + 𝐋 ∙ ,𝝈-(𝒅𝒗 ) + ,𝝈-(𝒅𝒗 ) ∙ 𝐋 − 𝑡𝑟(𝐋),𝝈-(𝒅𝒗 )

𝐝𝑛+1 = 𝐝𝑛+1 𝑣+1 𝑤ℎ𝑒𝑛 ‖∆𝐝𝑣 ‖ < 𝑒

The original NMM

CR IP T

2𝐌 2𝐌 𝑛 ( 2 + 𝐊(𝐝𝑛 )) ∆𝐝 = 𝐟 𝑒𝑥𝑡 (𝐝𝑛+1 ) − 𝐟 𝑖𝑛𝑡 (𝐝𝑛+1 ) + 𝐝̇ ∆𝑡 ∆𝑡 𝐝𝑛+1 = 𝐝𝑛 + ∆𝐝

AC

CE

PT

ED

M

AN US

𝝈𝑛+1 = 𝝈𝑛 + 𝐂𝐁 △ 𝒅