Rotation errors in numerical manifold method and a correction based on large deformation theory

Abstract

Numerical manifold method (NMM) is an effective method for simulating block system, however, significant errors are found in its simulation of rotation problems. Three kinds of errors, as volume expansion, stress vibration, and attenuation of angular velocity, were observed in the original NMM. The first two kind errors are owing to the small deformation assumption and the last one is due to the numerical damping. A large deformation NMM is proposed based on large deformation theory. In this method, the governing equation is derived using Green strain, the large deformation iteration and the open-close iteration are combined, and an updating strategy is proposed. The proposed method is used to analyze block rotation, beam bending, and rock falling problems and the results prove that all three kinds of errors are eliminated in this method.

Keywords

numerical manifold method / rotation / large deformation / Green strain / open-close iteration

Introduction

The numerical manifold method (NMM) proposed by Shi [¹] is an effective method for simulations of jointed rock mass. As a generalization of the Discontinuous Deformation Analysis (DDA) [²] that dealing with block systems, NMM provides a unified framework for simulations of both continuous media and discontinuous media [³] and is a method capable for large displacement simulation and contact analysis.

So far, great progress has been made to improve the accuracy of the NMM. In NMM, the deformation field is described by a weighted cover system which is similar to the later concept of “partition of unity” [⁴] in Finite Element Method (FEM). Therefore, a higher order NMM can be easily established by using higher order polynomials in the local covers, in the weight functions or in both. A second order NMM with linear weight functions is provided in Ref. [⁵], further, a stiffness matrix automatically derived algorithm for higher order NMM was proposed in Ref. [⁶]. However, similar to the partition of unity based on FEM (PUFEM) [⁴,⁷], a numerical problem “rank deficiency” appeasers in the NMM procedure with non-constant local covers [⁸] and a special matrix solver is required. The concept “cover” in the NMM is more general than the element in PUFEM. 1) The shape of cover is independent of physical body. This could benefit to some moving boundary problem, e.g., crack growth and intersection [⁹–¹²]. 2) The cover in NMM could totally overlap or partly overlap [¹³,¹⁴]. 3) The concept “cover” generalizes the element used in numerical simulation, e.g., the element with a combined lumped mass matrix [¹⁵]. In addition, similar to the Extended Finite Element Method (XFEM), some special covers have been introduced into NMM for local enrichment, such as the cover for the field near crack tips (e.g., Ref. [⁹]) and the material discontinuities (e.g., Ref. [¹⁶]).

Promising developments besides the improvements of cover system also have been made, including: 1) 3D NMM [¹⁷–¹⁹]; 2) Better contact accuracy, such as replacing the penalty function contact algorithm (see Ref. [¹], Chapter 6) with a new scheme based on augmented Lagrangian method [²⁰] or complementarity theory [²¹], and the general solution of contact problem-contact theory [²²]; 3) Other developments, such as the recursion formula of the simplex integration [²³] and the explicit NMM [²⁴,²⁵].

This study is motivated by the phenomenon “unreasonable volume expansion after large rotation” in DDA procedures (e.g., Ref. [²⁶].). Some early researchers have attributed the volume expansion of DDA to the fact that the linear displacement field cannot precisely describe the rotation field, which contains trigonometric function (sin and cos) and can be expressed as,

(1)\begin{array}{l}{u}_t=u_{\mathrm{0}}-(y-y_{\mathrm{0}})\sin \alpha +(x-x_{\mathrm{0}})(\cos \alpha -\mathrm{1})\hfill \\ v_t=v_{\mathrm{0}}+(x-x_{\mathrm{0}})\sin \alpha +(y-y_{\mathrm{0}})(\cos \alpha -\mathrm{1})\hfill \end{array},

where u_t, v_t are the displacements of point (x, y), and Eq. (1) denotes the displacement of point (x, y) which is first rotating around (x_{\mathrm{0}},y_{\mathrm{0}}) with a rotation angle of \alpha and then has a parallel translation of (u_{\mathrm{0}},v_{\mathrm{0}}). Some remedial measures for DDA has been proposed, such as post modification after each step [²⁶], replacing \cos \alpha by \mathrm{1}-\mathrm{0.5}{\alpha }^{\mathrm{2}} [²⁷] or replacing \cos \alpha -\mathrm{1}by -{\sin }^{\mathrm{2}}\alpha /(\mathrm{1}+\cos \alpha ) [²⁸], and proved to be effective in decreasing the volume expansion. However, These measures have some limitations: 1) The volume expansion in NMM/DDA is due to the accumulating error from the small deformation assumption, that is, any rotation with zero small deformation strain will always lead to a volume expansion, and thus block volume expands step by step in dynamic procedure. Because strains are not included in Eq. (1), the corrections based on Eq. (1) are not complete; 2) The displacement function of DDA has similar expression with Eq. (1) and can be easily to be modified, yet, NMM not.

To correct the rotation errors for NMM, Fan et al. [²⁹] proposed a modified NMM using strain-rotation decomposition (similarly, a modified DDA is also proposed in Ref. [³⁰].) yet trigonometric functions are involved and approximation is required in this method. Wei and Jiang [³¹] proposed a modified NMM using large deformation theory. However, the contact issue is not addressed in Ref. [³¹]. In summary, the main limitation of previous researches will be following three questions: 1) Why volume expansion? How to quantitatively interpret (or predict) the volume expansion? 2) Is there any rotation error other than the volume expansion? 3) How to ensure the accuracy of both rotations and contacts? The answers will be addressed in this paper.

In this study, following the idea of the updated Lagrangian scheme used in FEM [³²,³³] and replacing the small strain with Green strain, a new NMM scheme is proposed to erase the rotation errors with considering the contact issue. The proposed large deformation NMM will be further examined by theoretical solutions and practical problems.

The volume expansion phenomenon in the NMM

A brief introduction of NMM

NMM proposed by Shi [¹] can be regarded as a superior FEM. It releases the constraints of physical boundaries using two subsets, as mathematical elements and physical elements to describe the approximation function and the physical body separately [³,³⁴], and can flexibly handle discontinuities using a robust contact algorithm named as “open-close” iteration.

NMM Element

A comparison between FEM, NMM, and DDA elements (with same order) is illustrated in Table 1. The major differences between these three elements are:

For FEM, the elements shapes are limited, either triangle or quadrilateral. The mesh is dependent on the boundaries.

For DDA, the element shape is the shape of block, i.e., each block is a DDA element.

For NMM, the element contains two subsets: 1) mathematic element defining the approximation function; 2) physical element representing the physical body. Therefore, each physical element is a subset of a mathematic element, e.g., the pentagon in Table 1 is covered by the triangle mathematic element. The local approximation function of a mathematical element can be defined as

(2)\left\{\begin{array}{l}u\hfill \\ v\hfill \end{array}\right\}={\displaystyle \mathrm{\sum }_{i=\mathrm{1},n}{w}_i(x,y)d_i},

where d_i represents the approximation function of i mathematical cover (i.e., the cover functions; if first order displacement function is used, a mathematical cover can be regarded as a FEM node) and is the unknown variables to be solved; w_i represents the weight of i mathematical cover (just like the shape function of i node, but it has many other options (e.g., Refs. [⁹,¹³,¹⁶]).

NMM first adopts such kind an idea that uses two subsets in element description. Later, similar idea was proposed by Hansbo and Hansbo [³⁵] and successfully applied to the simulation of contacts, material discontinuous and crack growth. In their method, the mathematic node that outside the physical body as phantom node. Thus, Hansbo and Hansbo’s method is also named as phantom node method by later researchers [³⁶,³⁷]. Basically, the main difference between phantom node method [³⁵] and the original NMM [¹] might be the contact matrix and the contact iteration.

Equilibrium equation

As same as that used in DDA, the equations of equilibrium in NMM is based on principle of minimum potential energy (or principle of virtual work), and is established on the deformed configuration of last time step:

(3)({\mathit{K}}_{\mathrm{0}}+{\mathit{K}}_{\mathrm{con}})\mathrm{\Delta }{\mathit{d}}_n+\mathit{M}{\ddot{\mathit{d}}}_n=\mathit{f}+{\mathit{f}}_{\sigma }+{\mathit{f}}_{\mathrm{con}},

where \mathrm{\Delta }\mathit{d} is the unknown displacement increment of the n\mathrm{th} time step; {\dot{\mathit{d}}}_{n-\mathrm{1}}is the velocity and {\ddot{\mathit{d}}}_n is the acceleration; {\mathit{K}}_{\mathrm{0}}^{} and {\mathit{K}}_{\mathrm{con}}^{} are the assembled stiffness matrices for the deformation and the contacts, respectively; \mathit{f}, {\mathit{f}}_{\sigma }, {\mathit{f}}_{\mathrm{con}} are the assembled force for external loads, initial stress and contacts, respectively. The time integration in the original NMM (and DDA) is performed by substituting {\ddot{\mathit{d}}}_n solved from Eq. (4) to the governing equation Eq. (3) (see Ref. [¹])

(4)\mathrm{\Delta }{\mathit{d}}_n={\dot{\mathit{d}}}_{n-\mathrm{1}}\mathrm{\Delta }t+\mathrm{0.5}{\ddot{\mathit{d}}}_n\mathrm{\Delta }{t}^{\mathrm{2}}\mathrm{.}

It can be noted that Eq. (4) is the same with the Newmark method with parameters \beta =\mathrm{1} and \mathit{\gamma }=\mathrm{0.5} [³⁸] which can be expressed as,

(5)\begin{array}{l}{\dot{\mathit{d}}}_n={\dot{\mathit{d}}}_{n-\mathrm{1}}+\left[(\mathrm{1}-\beta ){\ddot{\mathit{d}}}_{n-\mathrm{1}}+\beta {\ddot{\mathit{d}}}_n\right]\mathrm{\Delta }t,\hfill \\ \mathrm{\Delta }{\mathit{d}}_n={\dot{\mathit{d}}}_{n-\mathrm{1}}\mathrm{\Delta }t+\left[(\mathrm{0.5}-\mathit{\gamma }){\ddot{\mathit{d}}}_{n-\mathrm{1}}+\mathit{\gamma }{\ddot{\mathit{d}}}_n\right]\mathrm{\Delta }{t}^{\mathrm{2}}\hfill \end{array}\mathrm{.}

For Newmark method expressed by Eq. (5), \beta with a value larger than 0.5 will induce numerical damping [³³]. Compared Eq. (4) to Eq. (5), it can be noted that Eq. (4) uses constant acceleration within a time step interval \mathrm{\Delta }t. Such numerical damping allows contact force oscillation to dissipate rapidly and results in a stable contact state [³⁹,⁴⁰], however, also leads to an attenuation of angular velocity in rotation simulation (as to be verified by the example in Section 3.1).

Discontinuity and contact

Many efforts have been paid to unify the continuous and discontinuous computation [⁴¹]. For discontinuous computation, contacts among physical bodies must be carefully considered. In NMM, an updated Lagrangian scheme is used (see Eq. (3)) to accurately search the contact position after large displacements. Then, an implicit scheme, named as “open-close” iteration, is used to deal with contact issues in NMM (and DDA). Three contact states, as open, sliding, and locking, give different stiffness matrices and forces, which will be assembled to {\mathit{K}}_{\mathrm{con}} and {\mathit{f}}_{\mathrm{con}} (see Eq. (3)) in each open-close iteration.

Usually, NMM/DDA uses small time step interval and is inertia dominated (Eq. (3) and referring to Ref. [¹]). Such inertia dominated phenomenon can benefit to the convergence of open-close iteration. With the help of contact springs, the convergence of open-close iteration can be always ensured if the time step is small enough in NMM (referring to Ref. [¹], Chapter 3).

In short, NMM proposed by Shi [¹] has following features:

1) Its element contains two subsets, as mathematic element defining the approximation function and the physical element describing integration area (i.e., physical body).

2) Its reference configuration is updated after each small time step and constant acceleration is used in the time step integration.

3) Its convergence of contact is ensured by the “open-close” iteration with time step interval small enough.

Compared to many meshless methods capable for large displacement simulation, such as the Smoothing Particle Hydrodynamics (SPH) which is capable for the interaction simulation between fluid and solid [⁴²], the crack particles method that is efficient for crack growth simulation [⁴³], and the Peridynamics method that adopts a simplified scheme based on integral form rather than the partial differential form [⁴⁴,⁴⁵], NMM is more flexible and benefits from that the approximation function (mathematical element) is free from the physical boundaries. In NMM, the analysis precision is majorly determined by the order or approximation function of the mathematical element. Meanwhile, many old problems, such as the difficulty of the mesh generation for complex geometry or the simulation of accurate boundary conditions, are no longer exists.

However, the updated Lagrangian scheme, the implicit time integration (Eq. (3)) and the strict convergence criterion of the open-close iteration make NMM be a procedure trustworthy but less efficient.

The volume expansion phenomenon during rotation process in the NMM

The original NMM [¹] is an implicit dynamic procedure and follows the small assumption used in FEM. In each time step, the displacement increment must be small enough to ensure the accuracy of calculation and the convergence of contact iteration.

Small deformation assumption (or small strain assumption)

For describing the deformation field, the following approximation equation is usually used in the original NMM/DDA:

(6)\{\begin{array}{l}{\mathit{\epsilon }}_x\approx {\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{x}}}\approx {\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{X}}}\hfill \\ {\mathit{\epsilon }}_y\approx {\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{y}}}\approx {\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{Y}}}\hfill \\ {\mathit{\gamma }}_{xy}\approx {\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{x}}}+{\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{y}}}\approx {\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{X}}}+{\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{Y}}}\hfill \end{array},

where the strain \left\{{\mathit{\epsilon }}_x,{\mathit{\epsilon }}_y,{\mathit{\gamma }}_{xy}\right\} is called small deformation strain and thus Eq. (6) is named as the small deformation assumption; the lower case (\mathit{x},\mathit{y}) is the body coordinates of deformed configuration; the upper case (\mathit{X},\mathit{Y}) is the body coordinate of reference configuration (see Table 1) and (\mathit{u},\mathit{v}) is the displacement, i.e.:

(7)(\mathit{x},\mathit{y})=(\mathit{X},\mathit{Y})+(\mathit{u},\mathit{v})\mathrm{.}

The strain error in Eq. (6) is limited in each step if the real displacement is small enough, however, such error can be accumulated by steps, and can lead to significant rotation errors.

A simple rotation simulation using the original NMM

As shown in Fig. 1, an equilateral triangle ABC is fixed at vertex A and applied with an initial angular velocity of {\omega }_{\mathrm{0}}=\mathrm{2}\mathrm{\pi }/s. The time step interval in the computation is used as \mathrm{\Delta }t=\mathrm{0.01}\mathrm{s}. Therefore, the triangle is supposed to turn a circle every 100 time steps with constant volume.

In this test, the boundary condition of the fixed vertex A was imposed by eliminating the freedoms (just like that adopted in FEM) instead of the penalty method (see Ref. [¹]). Therefore, no errors are involved in the boundary condition.

Referring to Fig. 2, significant volume expansion was observed in the simulation using original NMM. Although the stain was zero in the whole process (each component less than 1e−7), about 43% volume error and −9.1% cumulative rotation angle error appeared after one circle. About 153% volume error and −23.0% cumulative rotation angle error appeared after three circles (see Fig. 2).

The rotation expansion associated with small deformation assumption

In fact, the rotation expansion phenomenon is associated with the small deformation assumption used in NMM. Since DDA use the same assumption (see Ref. [²]), such phenomenon also exists in DDA.

A graphic explanation of the rotation expansion phenomenon

A square block OABC (the dashed line) with unit length is shown in Fig. 3. Without loss of generality, the point O is assumed to be fixed and a rigid rotation of block OABC is considered. Let OA\prime B\prime C represents the deformed square, then an incremental rotation angle of \alpha is used in the simulation and no other deformation exists, the displacement at any point can be calculated by Eq. (1), as shown in Fig. 3(a).

The strain should be zero for such a rigid rotation problem in Fig. 3. If small deformation assumption is adopted (Eq. (6)), we have,

(8){\mathit{\epsilon }}_x={\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{x}}}=\mathrm{0},{\mathit{\epsilon }}_y={\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{y}}}=\mathrm{0},{\mathit{\gamma }}_{xy}={\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{x}}}+{\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{y}}}=\mathrm{0.}

Note that the displacement at the point A is \left({\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{x}}},{\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{y}}}\right), thus with Eq. (4), the point A\prime is just on the line AB; similarly, the point C\prime is on the extended line of BC (see Fig. 3(b)). Clearly, the edge length of deformed square is \sqrt{\mathrm{1}+{\tan }^{\mathrm{2}}\alpha }. Therefore in this condition, the calculated volume strain {\mathit{\epsilon }}_v, i.e., the strain error, is

(9){\mathit{\epsilon }}_v={\tan }^{\mathrm{2}}\alpha \approx {\alpha }^{\mathrm{2}}+o\left({\alpha }^{\mathrm{4}}\right)\mathrm{.}

The volume will expand during rotation and the volume expansion in each step is approximately equal to the square of incremental rotation angle.

A theoretical explanation of the rotation expansion phenomenon

Green strain can also be used to explain such rotation expansion phenomenon. Let \mathit{E} denotes the Green strain tensor [³²,⁴⁶], \mathit{e} denotes the small deformation strain tensor, \mathit{u} denotes the displacement and \mathit{X} denotes the body coordinate. We have

(10)E_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial u_i}{\partial X_j}}+{\displaystyle \frac{\partial u_j}{\partial X_i}}+{\displaystyle \frac{\partial u_i}{\partial X_k}}{\displaystyle \frac{\partial u_j}{\partial X_k}}\right)=e_{ij}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\frac{\partial u_i}{\partial X_k}{\displaystyle \frac{\partial u_j}{\partial X_k}}\mathrm{.}

In small deformation assumption, only the linear-part \mathit{e} is used, i.e., the nonlinear term {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\partial u_i}{\partial X_k}}{\displaystyle \frac{\partial u_j}{\partial X_k}} is the strain error. Specifically, for a two-dimensional analysis, let \mathit{\epsilon } represents the vector form of Green strain tensor \mathit{E}, {\mathit{\epsilon }}_{\mathrm{0}} the vector of small strain \mathit{e}, {\mathit{\epsilon }}_{\mathrm{N}} the vector of high order strain, i.e., nonlinear term. We have

(11)\mathit{\epsilon }={\left\{{\epsilon }_x,{\epsilon }_y,{\mathit{\gamma }}_{xy}\right\}}^{\mathrm{T}}={\mathit{\epsilon }}_{\mathrm{0}}+{\mathit{\epsilon }}_N,

(12)\begin{array}{l}{\mathit{\epsilon }}_{\mathrm{0}}={\left\{{\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{X}}},{\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{Y}}},{\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{X}}}+{\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{Y}}}\right\}}^{\mathrm{T}}\hfill \\ \\ {\mathit{\epsilon }}_{\mathrm{N}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\{{\left({\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{X}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{X}}}\right)}^{\mathrm{2}},{\left({\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{Y}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{Y}}}\right)}^{\mathrm{2}},}^{\mathrm{}}\hfill \end{array}{2\left({\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{X}}}\frac{\partial \mathit{u}}{\partial \mathit{Y}}+\frac{\partial \mathit{v}}{\partial \mathit{X}}{\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{Y}}}\right)\}}^{\mathrm{T}},

where (\mathit{X},\mathit{Y}) and (\mathit{u},\mathit{v}) represent the body coordinates and the displacement, respectively.

Referring to Fig. 3(b), the following relation is valid,

(13)\begin{array}{l}{\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{X}}}=-{\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{Y}}}=\tan \alpha ,\hfill \\ \\ \alpha \approx {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{X}}}-{\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{Y}}}\right)\mathrm{.}\hfill \end{array}

Therefore, if small deformation strain is zero, i.e., {\mathit{\epsilon }}_{\mathrm{0}}=\left\{\mathrm{0}\right\} (referring to Eq. (8)), the strain error, i.e., the nonlinear part of strain, is

(14){2\left({\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{X}}}\frac{\partial \mathit{u}}{\partial \mathit{Y}}+\frac{\partial \mathit{v}}{\partial \mathit{X}}{\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{Y}}}\right)\}}^{\mathrm{T}}{\mathit{\epsilon }}_{\mathrm{N}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\{{\left({\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{X}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{X}}}\right)}^{\mathrm{2}},{\left({\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{Y}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{Y}}}\right)}^{\mathrm{2}},}^{}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\tan }^{\mathrm{2}}\alpha ,{\tan }^{\mathrm{2}}\alpha ,\mathrm{0}\right)}^{\mathrm{T}},

and the real volume strain would be contributed from {\epsilon }_{\mathrm{N}}, as

(15){\mathit{\epsilon }}_v={\tan }^{\mathrm{2}}\alpha \approx {\alpha }^{\mathrm{2}}+o\left({\alpha }^{\mathrm{4}}\right)\mathrm{.}

In short, giving a small rotation angle \alpha and zero small deformation strain, the strain error of small deformation assumption is \left\{{\mathit{\epsilon }}_x,{\mathit{\epsilon }}_y,{\mathit{\gamma }}_{xy}\right\}\approx \left\{\mathrm{0.5}{\alpha }^{\mathrm{2}},\mathrm{0.5}{\alpha }^{\mathrm{2}},\mathrm{0}\right\}, and therefore the volume strain error is {\mathit{\epsilon }}_V\approx {\alpha }^{\mathrm{2}}. False volume expansion will happen in the rotation process if small strain assumption is adopted.

Rotation errors in the NMM for the rigid rotation problem

Volume error and angle error in the rotation simulation by NMM

In this section, the results of numerical simulation are compared with a theoretical solution to verify whether the volume expansion in the NMM is owing to small deformation assumption. Based on Eq. (15), the volume under small deformation assumption can be calculated as,

(16)V_n\approx (\mathrm{1}+{\alpha }_{n-\mathrm{1}}^{\mathrm{2}})V_{n-\mathrm{1}}\approx (\mathrm{1}+{\alpha }_{n-\mathrm{1}}^{\mathrm{2}})(\mathrm{1}+{\alpha }_{n-\mathrm{2}}^{\mathrm{2}})\mathrm{\cdots }(\mathrm{1}+{\alpha }_{\mathrm{0}}^{\mathrm{2}})V_{\mathrm{0}},

where {\alpha }_i^{} represent the incremental rotation angle at i\mathrm{th} time step and V_i represent the volume.

While in the NMM simulation, the computation settings are listed in Table 2. Five tests, with time step interval \mathrm{\Delta }t ranges from 0.01–0.0001s (i.e., theoretical incremental rotation angle \alpha ranges from 3.6° to 0.036°), were investigated.

The results of five tests are summarized in Table 3. Both the volume error and the rotation angle error are linear increase with the theoretical incremental rotation angle \alpha (see Fig. 4).

The volume error can be calculated as V_n/V_{\mathrm{0}}-\mathrm{1} (referring to Eq. (16)). As shown in Fig. 2(b) and Table 3, the real volume error does fit the volume error calculated by Eq. (16). Thus, this result verified that small deformation assumption is the reason of the rotation expansion phenomenon.

Both the volume error and the angle error have strong linear relationship (R^{\mathrm{2}}>\mathrm{0.997}) with the rotation angle increment used in the simulation (see Fig. 4). However, it is time-consuming to maintain an extremely low rotation angle increment. Even the incremental rotation angle is taken as \alpha =\mathrm{0.036}°, i.e., 10000 steps for a rotation of one circle, the volume error still has 0.4%.

The linear relationship between the volume errors and the incremental rotation angle α can be explained. If incremental rotation angle \alpha is much less than 1, we have

(17)V_n=(\mathrm{1}+{\alpha }_{n-\mathrm{1}}^{\mathrm{2}})(\mathrm{1}+{\alpha }_{n-\mathrm{2}}^{\mathrm{2}})\mathrm{\cdots }(\mathrm{1}+{\alpha }_{\mathrm{0}}^{\mathrm{2}})V_{\mathrm{0}}\approx (\mathrm{1}+{\alpha }_{n-\mathrm{1}}^{\mathrm{2}}+{\alpha }_{n-\mathrm{2}}^{\mathrm{2}}+\mathrm{\cdots }+{\alpha }_{\mathrm{0}}^{\mathrm{2}})V_{\mathrm{0}}\approx (\mathrm{1}+n{\alpha }_{}^{\mathrm{2}})V_{\mathrm{0}},

and n is the step number for a circle rotation, thus

(18)n\approx \mathrm{2}\mathrm{\pi }/\alpha .

Substitute Eq. (18) into Eq. (17), the block volume after a rotation of one circle can be estimated by

(19)

where V_{\mathrm{circle}} is the block volume after one circle, and thus \mathrm{2}\mathrm{\pi }\alpha is a estimation of the volume error.

Besides the volume error, a negative rotation angle error is found in Table 3. Use test 2 as an example, the curves of angular velocity and the rotation angle versus time are plotted in Fig. 5, where the simulation time is 10 s and the block should finish ten circles in theory. It can be found that the angular velocity was decreasing in the simulation and the block only finishes 7.8 circulars after 10 s. Such attenuation of angular velocity is supposed to be owing to the numerical damping in the original NMM scheme (referring to Section 2.1). This issue will be further discussed and corrected later in this paper.

Stress vibration phenomenon during rotation process

Besides the volume error and the rotation angle error (attenuation of angular velocity), stress vibration is also found in upper simulation using NMM.

For the simple rotation example illustrated in Fig. 1, the stress ({\mathit{\sigma }}_x, {\mathit{\sigma }}_y or {\mathit{\tau }}_{xy}) will be uniform for the single element used and should complete two periods after a circle. For example, {\mathit{\sigma }}_x here should reach maximum value at rotation angle of −30° (that is, 330°, where the center of gravity is on the positive axis for x) and 150° because of the centrifugal force in the rotation process. Figure 6 shows the stress variation with time of test 1. (The wave of {\mathit{\sigma }}_y is similar to that of {\mathit{\sigma }}_x, thus only {\mathit{\sigma }}_x and {\mathit{\tau }}_{xy} are illustrated.)

Referring to Fig. 6, it can be found that the stress containing two waves, as

1) The low frequency one, which approximately completes two periods in 1 s and is as expected.

2) The high frequency one, which has a remarkable amplitude. It can be found that the stress vibration with high frequency has a period equal to 2Dt, i.e., 2 times of time step interval. As shown in Fig. 7, when Dt is changing from 0.01 to 0.005 s, the period of stress vibration changes from 0.02 to 0.01 s. Consequently, the stress vibration cannot be eliminated by decreasing the time step interval \mathrm{\Delta }t.

The low frequency waves versus the time and the rotation angle are plotted in Fig. 8. The phases of waves differ with each other in time domain (Fig. 8(a)) due to the decrease of incremental rotation angle induced by numerical damping, and as expected, the phase differences disappear in the rotation angle domain (Fig. 8(b)), nevertheless, significant amplitude differences are still observed among the three waves. According to the later discussion in Section 5.2, the low frequency wave with small \mathrm{\Delta }t = 0.0001 s is nearly as same as that obtained in large deformation analysis. The stress calculation with large \mathrm{\Delta }t is not correct.

In summary, an unreasonable stress vibration is observed during the rotation simulation. Such vibration neither disappeared due to the increasing of time, nor eliminated by using smaller time step interval (e.g., Fig. 7). Such stress vibration may cause significant error in the stress calculation and is not acceptable.

A summary of rotation error

Three types of rotation errors exist in the original NMM, and are listed below:

1) Volume expansion. Volume expansion is a direct result of small deformation assumption. The volumetric strain error within a time step is approximately equal to the square of rotation angle increment and the error will be cumulative in the simulation.

2) Speed slowing. Because of the numerical damping, the angular velocity decreases with the simulation steps.

3) Stress vibration. A stress vibration with remarkable amplitude and a period of \mathrm{2}\mathrm{\Delta }t is observed during the rotation simulation. The vibration cannot be eliminated by decreasing the time step interval and will cause unacceptable accuracy of stress calculation.

The implementation of a large deformation NMM

To deal with the volume expansion phenomenon, many commercial FEM softwares, such as Abaqus and Comsol, have provided an additional option using large deformation theory. Meanwhile, some DDA researchers have found the rotation expansion phenomenon and provide some remediation methods, such as [²⁷,⁴⁷]. However, as indicated in Sections 1 and 3, a correction still use small deformation assumption is insufficient and cannot erase the stress vibration.

To erase the rotation errors in NMM, a large deformation NMM (LDNMM) is proposed and its basic idea is as follows: Green strain is used instead of the linear small deformation strain and a large deformation iteration is introduced to deal with the nonlinearity. There are two difficulties to be overcome as:

1)The derivation of governing equation with the use of Green strain for NMM;

2)The coupling of large deformation iteration with contact iteration in NMM.

The governing equation of the large deformation NMM

Giving a virtual displacement \delta \mathit{d} and based on the principle of virtual work, the governing equation of can be written as,

(20){\displaystyle \int \delta {\mathit{\epsilon }}^{\mathrm{T}}(\mathit{\sigma }+{\mathit{\sigma }}_{\mathrm{0}})\mathrm{d}A}+\mathrm{\cdots }=\mathit{f}\delta \mathit{d},

where {\sigma }_{\mathrm{0}} is the initial stress; A is the integration area; f is the extern force and {\displaystyle \int \delta {\mathit{\epsilon }}^{\mathrm{T}}(\mathit{\sigma }+{\mathit{\sigma }}_{\mathrm{0}})dA} is the virtual work from deformation, which can be divided into 2 parts

(21){\displaystyle \int \delta {\mathit{\epsilon }}^{\mathrm{T}}D\mathit{\epsilon }\mathrm{d}A}\mathrm{and}{\displaystyle \int \delta {\mathit{\epsilon }}^{\mathrm{T}}{\mathit{\sigma }}_{\mathrm{0}}\mathrm{d}A},

where the first part gives the deformation stiffness matrix and the second part refers to the initial stress matrix. The only things necessary is to derive the expression of \delta \mathit{\epsilon } and \mathit{\epsilon }.

If the Green strain (Eq. (10)) is used instead of the linear small strain, the linear part is retaining as same as that in the original NMM and the nonlinear part of strain should be added. The nonlinear part of strain is

(22){\mathit{\epsilon }}_{\mathrm{N}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\{{\left({\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{X}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{X}}}\right)}^{\mathrm{2}},{\left({\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{Y}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{Y}}}\right)}^{\mathrm{2}},}^{\mathrm{}}{\mathrm{2}\left(\frac{\partial \mathit{u}}{\partial \mathit{X}}{\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{Y}}}+\frac{\partial \mathit{v}}{\partial \mathit{X}}{\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{Y}}}\right)\}}^{\mathrm{T}},

and can be also expressed as

(23){\mathit{\epsilon }}_{\mathrm{N}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\begin{array}{llll}{\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{X}}}\hfill & {\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{X}}}\hfill & \mathrm{0}\hfill & \mathrm{0}\hfill \\ \mathrm{0}\hfill & \mathrm{0}\hfill & {\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{Y}}}\hfill & {\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{Y}}}\hfill \\ {\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{Y}}}\hfill & {\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{Y}}}\hfill & {\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{X}}}\hfill & {\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{X}}}\hfill \end{array}\right]\left\{\begin{array}{l}{\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{X}}}\hfill \\ {\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{X}}}\hfill \\ {\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{Y}}}\hfill \\ {\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{Y}}}\hfill \end{array}\right\}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{A}\theta ,

where \theta ={\left\{{\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{X}}},{\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{X}}},{\displaystyle \frac{\partial \mathit{u}}{\partial \mathit{Y}}},{\displaystyle \frac{\partial \mathit{v}}{\partial \mathit{Y}}}\right\}}^{\mathrm{T}}. It is noted that the matrix \mathit{A} is determined by \theta, and then the differential of nonlinear strain can be derived as

(24)\delta {\mathit{\epsilon }}_{\mathrm{N}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\delta \mathit{A}\theta +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{A}\delta \theta =\mathit{A}\delta \theta ,

and \theta can be calculated from the displacement as

(25)\theta =\mathit{G}\mathrm{\Delta }{\mathit{d}}^e,

where \mathrm{\Delta }{\mathit{d}}^e is the unknown incremental displacement of a NMM element e; \mathit{G} is a coefficient matrix which can be derived from Eq. (2) and just similar to the coefficient matrix for linear strain {\mathit{B}}_0 (see Ref. [¹], Section 5.1) G is independent of \mathrm{\Delta }\mathit{d} and can be regarded as a constant matrix in the derivation.

Substituting Eq. (25) into Eqs. (23) and (24), we have

(26)\begin{array}{l}{\mathit{\epsilon }}_{\mathrm{N}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{A}\mathit{G}\mathrm{\Delta }{\mathit{d}}^e={\mathit{B}}_{\mathrm{N}}\mathrm{\Delta }{\mathit{d}}^e\hfill \\ \delta {\mathit{\epsilon }}_{\mathrm{N}}=\mathit{A}\mathit{G}\mathrm{\Delta }{\mathit{d}}^e=\mathrm{2}{\mathit{B}}_{\mathrm{N}}\mathrm{\Delta }{\mathit{d}}^e\hfill \end{array},

where {\mathit{B}}_{\mathrm{N}} is the nonlinear part of strain matrix, as

(27){\mathit{B}}_{\mathrm{N}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{A}\mathit{G}\mathrm{.}

Therefore, substituting Eq. (26) into Eq. (21), the following relation valids:

(28)$$\begin{gathered} \int {\delta \mathit{\epsilon }}^{\mathrm{T}}\mathit{D}\mathit{\epsilon }\mathrm{d}A=\delta {\left({\mathrm{\Delta }\mathit{d}}^e\right)}^{\mathrm{T}}\int {({\mathit{B}}_0+{2\mathit{B}}_{\mathrm{N}})}^{\mathrm{T}}\mathit{D}({\mathit{B}}_0+{\mathit{B}}_{\mathrm{N}})\mathrm{d}A\mathrm{\Delta }{\mathit{d}}^e \\ \int {\delta \mathit{\epsilon }}^{\mathrm{T}}{\mathit{\sigma }}_0\mathrm{d}A=\delta {\left({\mathrm{\Delta }\mathit{d}}^e\right)}^{\mathrm{T}}\int {({\mathit{B}}_0+{2\mathit{B}}_{\mathrm{N}})}^{\mathrm{T}}{\sigma }_0\mathrm{d}A, \end{gathered}$$

where {\mathit{B}}_{\mathrm{0}} is the coefficient matrix of the linear strain and is as same as that used in original NMM (see Ref. [¹]). Hence, the stiffness matrix of the element e is

(29){\mathit{K}}^e={\displaystyle \int {({\mathit{B}}_{\mathrm{0}}+\mathrm{2}{\mathit{B}}_{\mathrm{N}})}^{\mathrm{T}}\mathit{D}({\mathit{B}}_{\mathrm{0}}+{\mathit{B}}_{\mathrm{N}})dA}={\mathit{K}}_{\mathrm{0}}^e+{\mathit{K}}_{\mathrm{N}}^e,

(30)\{\begin{array}{c}{\mathit{K}}_0^e=\int {{\mathit{B}}_0}^{\mathrm{T}}\mathit{D}{\mathit{B}}_0\mathrm{d}A\hfill \\ {\mathit{K}}_{\mathrm{N}}^e=\int {{2\mathit{B}}_{\mathrm{N}}}^{\mathrm{T}}\mathit{D}{\mathit{B}}_0+{2\mathit{B}}_{\mathrm{N}}\mathit{D}{\mathit{B}}_{\mathrm{N}}+{{\mathit{B}}_0}^{\mathrm{T}}\mathit{D}{\mathit{B}}_{\mathrm{N}}\mathrm{d}A\hfill \end{array},

and the force from initial stress of element e is

(31){\mathit{f}}_{\sigma }^e=-{\displaystyle \int {({\mathit{B}}_{\mathrm{0}}+\mathrm{2}{\mathit{B}}_{\mathrm{N}})}^{\mathrm{T}}{\sigma }_{\mathrm{0}}\mathrm{d}A}\mathrm{.}

Thus, the governing equation in matrix form at n\mathrm{th} time step is

(32)\left({\mathit{K}}_{\mathrm{0}}+{\mathit{K}}_{\mathrm{N}}+{\mathit{K}}_{\mathrm{con}}\right)\mathrm{\Delta }{\mathit{d}}_n+\mathit{M}{\ddot{\mathit{d}}}_n=\mathit{f}+{\mathit{f}}_{\sigma }+{\mathit{f}}_{\mathrm{con}},

where {\mathit{K}}_{\mathrm{0}}^{}, {\mathit{K}}_{\mathrm{con}}^{}\mathit{f}, {\mathit{f}}_{\sigma }, and {\mathit{f}}_{\mathrm{con}} are as same as those involved in original NMM (Eq. (3)); {\mathit{K}}_{\mathrm{N}} and {\mathit{f}}_{\sigma } are two newly introduced terms to consider the nonlinear strain. Different with that {\mathit{B}}_{\mathrm{0}} is a constant matrix, {\mathit{B}}_{\mathrm{N}} is not a constant matrix (Eq. (27)). Hence, iteration is required to achieve the convergence of nonlinear strain in the large deformation analysis.

The large deformation iteration

In the large deformation iteration of Eq. (32), \theta is chosen as the key iteration variable and the iteration procedure is implemented as follows:

1) An initial value {\theta }_{\mathrm{0}} is given, e.g., taken as the convergent value obtained in the last time step.

2) Matrix A is calculated from {\theta }_i, then {\mathit{B}}_{\mathrm{N}}, {\mathit{K}}_{\mathrm{N}}, and {\mathit{f}}_{\sigma } can be obtained using Eqs. (27), (30), and (31), respectively, where i is the number of iteration step and is counting from 0.

3) The governing equation (Eq. (32)) is solved, resulting in a displacement field, i.e., \mathrm{\Delta }\mathit{d}.

4) A new value of {\theta }_{i+\mathrm{1}} is calculated by Eq. (25) with using coefficient matrix \mathit{G} and \mathrm{\Delta }\mathit{d}.

5) The difference between {\theta }_{i+\mathrm{1}} and {\theta }_i are checked. If the error is larger than the tolerance, {\theta }_{i+\mathrm{1}} will be substituted into step 2 to replace {\theta }_i and the steps 2 to 5 will be repeated.

The upper iteration is in fact a fixed-point iteration which adopts the secant matrix instead of the tangent matrix that is used in Newton Raphson iteration. It is noted that the contact iteration [¹] (i.e., open-close iteration) also is a fixed-point iteration, thus the large deformation can be easily coupled with the contact iteration.

To present a clear vision of the large deformation iteration, both the dynamic and static simulation are discussed as follows:

Case 1: static simulation without contact

For a static problem, the govern equation (Eq. (32)) can be simplified as,

(33)\left({\mathit{K}}_{\mathrm{0}}+{\mathit{K}}_{\mathrm{N}}\right)d=f\mathrm{.}

The time step for such problem is meaningless. The large deformation iteration can be implemented as follows:

Step 1: Let i=\mathrm{0}, the initial value of \theta is given as {\theta }_{\mathrm{0}}^{}=(\mathrm{0},\mathrm{0},\mathrm{0},\mathrm{0});

Step 2: With {\theta }_i^{}, Eq. (27) is used to calculate matrix {\mathit{B}}_{\mathrm{N}}, Eq. (30) is used to calculate {\mathit{K}}_{\mathrm{N}}^{} and Eq. (33) is used to solve \mathit{d};

Step 3: Substitute d into Eq. (25) to compute {\theta }_{i+\mathrm{1}}^{};

Step 4: Let {\theta }_i^{} equals to {\theta }_{i+\mathrm{1}}^{} and repeat step 2 to 4 until the difference between {\theta }_{i+\mathrm{1}}^{} and {\theta }_i^{} is smaller than the tolerance.

Case 2: dynamic simulation without contacts

For a dynamic problem without contacts, the governing equation (Eq. (32)) can be simplified as,

(34)\left({\mathit{K}}_{\mathrm{0}}+{\mathit{K}}_{\mathrm{N}}\right)\mathrm{\Delta }{\mathit{d}}_n+\mathit{M}{\ddot{\mathit{d}}}_n=\mathit{f}+{\mathit{f}}_{\sigma }\mathrm{.}

In each time step, it is effective to use the same procedure for static problem (Case 1). However, a better initial value can help to make the iteration more efficient. The converged \theta at last time step is a good candidate, i.e.,

(35){\theta }_{n,0}^{}={\theta }_{n-\mathrm{1}}^{},

where {\theta }_{n,0}^{} is the initial value of \theta in n\mathrm{th} time step; {\theta }_{n-\mathrm{1}}^{} is the converged \theta in (n-1)th time step. If time interval do not change in n\mathrm{th} time step, The \theta of n\mathrm{th} time step should be close to that of (n-1)th due to system inertia.

The coupling of large deformation iteration and contact iteration

If contacts exist in the system, the large deformation iteration must be coupled with the open-close iteration in a practical NMM simulation. In this study, the open close iteration is embedded in the large deformation iteration (As shown in the flow chart Fig. 9) where a double nested loop is used with the outer loop for large deformation iteration and the inner loop for open-close iteration.

Case 3: practical NMM simulation involving contacts

A general simulation of NMM is dynamic and involving contacts. It follows the governing equation expressed in Eq. (32). When contacts are involved, open-close iteration is used and often needs to select a time step interval small enough to ensure the convergence. Hence, the time step interval is often non-constant. In this condition, system inertia can be considered by the velocity as follows:

(36)v_n^{\theta }={\displaystyle \frac{{\theta }_n^{}}{\mathrm{\Delta }{t}_n}},

and

(37)v_{n,\mathrm{0}}^{\theta }=v_{n-\mathrm{1}}^{\theta },

where v_{n,\mathrm{0}}^{\theta } is the initial value in n\mathrm{th} time step; v_{n-\mathrm{1}}^{\theta } is the converged value in (n-1)th time step. Then, the iteration processes can be implemented as follows:

Step 1: Let i=\mathrm{0}, and give initial value v_{n,\mathrm{0}}^{\theta }.

Step 2: Calculate {\theta }_{n,i}^{} by Eq. (36). Then, use Eqs. (27), (30), and (31) to calculate {\mathit{B}}_{\mathrm{N}}, {\mathit{K}}_{\mathrm{N}}^{}, and {\mathit{f}}_{\sigma }, respectively.

Step 3: Carry out open-close iteration. If the open-close iteration is not converged within 8 iterations, time step interval will be reduced by 0.35 times and go back to Step 2. After the open-close iteration is converged, go to next step.

Step 4: Solve Eq. (32) and Substitute the solved displacement \mathrm{\Delta }\mathit{d} into Eq. (25) to recalculate {\theta }_{n,i+\mathrm{1}}^{}.

Step 5: Substitute {\theta }_{n,i+\mathrm{1}}^{} into Eq. (36) to calculate v_{n,i\mathrm{+1}}^{\theta } and check the convergence. If the difference between v_{n,i\mathrm{+1}}^{\theta } and v_{n,i}^{\theta } for any element is larger than the tolerance, set i=i+\mathrm{1} and return to Step 2.

So far, a large deformation NMM is established using Green strain. In the new scheme, the large deformation iteration is coupled with the open-close iteration. In addition, because of the inertia effect, a better selection of initial value can benefit to the convergence of large deformation iteration.

Special discussion on the explicit time integration

Obviously implicit time integration is used in upper large deformation iteration. If an explicit time integration scheme is used, the governing equation can be

(38)M_{\mathrm{lumped}}\mathrm{\Delta }d=f_{\mathrm{res}},

where M_{\mathrm{lumped}} and f_{\mathrm{res}} represents the assembled lumped mass matrix (which is a diagonal matrix) and the assembled resultant force, respectively. Subsequently in the large deformation iteration, the element force from deformation can be calculated either by,

(39)f_{\mathrm{def}}^{\mathrm{e}}=-{\displaystyle {\int }_{}^{}(B_{\mathrm{0}}+\mathrm{2}{B}_{\mathrm{N}})\sigma }\mathrm{d}A,

or by

(40)f_{\mathrm{def}}^{\mathrm{e}}=f_{\sigma }^{\mathrm{e}}-(K_{\mathrm{0}}^{\mathrm{e}}+K_{\mathrm{N}}^{\mathrm{e}})\mathrm{\Delta }{d}^{\mathrm{e}},

where f_{\mathrm{def}}^{\mathrm{e}} is the element force from deformation and should be assembled into the resultant force f_{\mathrm{res}}.

An explicit time integration algorithm usually requires more iteration numbers. Generally, a small time step interval (such as 1e−5 s, or even as small as 1e−8 s), a large tolerance of unbalanced force and a lower stiffness for contact spring can benefit to the computational stability and accuracy. However, because only the diagonal matrix M_{\mathrm{lumped}} is involved in the explicit time integration, the computation cost can be speed up significantly. In this study, only the implicit algorithm is addressed. An explicit large deformation NMM can be a valuable work in future and it seems to be more efficient with using parallel computation and suitable for the simulation with large number of freedoms [²⁴,²⁵].

The updating rule in large deformation NMM

As same as DDA, NMM also adopts updated Lagrangian scheme, i.e., the variables are defined on the deformed configuration of last time step. Not only the coordinates of substance (i.e., mesh) should be updated, but also all the variables associated with physical body must be updated before carrying out a new step simulation. Such updated Lagrangian scheme has two benefits: 1) it can be easily coping with the deformation with large displacement; 2) it can help to accurately search the contact position.

In original NMM, only the coordinate of mesh and the density of an element are updated. In NMM, the physical body is served as an integration area. If the integration area is deformed and updated, not only the density, but also all the variables depending on the integration area should be updated either. For example, the initial stress should be updated because that {\sigma }_n^{} is defined on the physical body at n\mathrm{th} step, not on the deformed physical body at (n+1)th step. In large deformation theory, the deformation gradient \mathit{F} is used to represent the mapping rule from the reference configuration to the deformed configuration [³³]. Let u denotes the incremental displacement, \mathit{X} the initial position of n\mathrm{th} time step and \mathit{x} the deformed position (in other words, \mathit{x} is the initial position of (n+1)th time step), we have

(41)\mathit{x}=\mathit{X}+\mathit{u}\mathrm{.}

Then, the deformation gradient \mathit{F} is written as [³³],

(42)\left[F_{ij}\right]=\left[{\displaystyle \frac{\partial x_i}{\partial X_i}}\right]\left[\begin{array}{cc}\hfill {\displaystyle \frac{\partial u}{\partial X}}+\mathrm{1}\hfill & \hfill {\displaystyle \frac{\partial u}{\partial Y}}\hfill \\ \hfill {\displaystyle \frac{\partial v}{\partial X}}\hfill & \hfill {\displaystyle \frac{\partial v}{\partial Y}}+\mathrm{1}\hfill \end{array}\right]\mathrm{.}

Clearly, deformation gradient F can be obtained by \theta. Then, the mass density of an element should be updated as

(43){\rho }_{n+\mathrm{1}}={\displaystyle \frac{{\rho }_n}{\left|\mathit{F}\right|}},

where \rho is the density of an element e; \left|\mathit{F}\right| is the determinant of \mathit{F} and Eq. (43) can give the same result with that used in original NMM as

(44){\rho }_{n+\mathrm{1}}={\displaystyle \frac{{\rho }_n{A}_n}{A_{n+\mathrm{1}}}},

where A represents the area of an element e.

The stress obtained in n\mathrm{th} time step is

(45){\sigma }_n={\sigma }_n^{\mathrm{ini}}+D({\mathit{B}}_{\mathrm{0}}+{\mathit{B}}_{\mathrm{N}})\mathrm{\Delta }{\mathit{d}}_n^e,

then updated as

(46){\sigma }_{n+\mathrm{1}}^{\mathrm{ini}}={\displaystyle \frac{\mathrm{1}}{\left|\mathit{F}\right|}}\mathit{F}{\sigma }_n^{}{\mathit{F}}^{\mathrm{T}},

{\sigma }_{n+\mathrm{1}}^{\mathrm{ini}} is the initial stress of (n+1)th time step. The detail derivation of Eq. (46) could be found in Ref. [³³], section 6.2.2, where {\sigma }_n^{} is the second Piola – Kirchoff stress and {\sigma }_{n+\mathrm{1}}^{\mathrm{ini}} the Cauchy stress. If under small rotation or strain, {\sigma }_n^{} could be close to {\sigma }_{n+\mathrm{1}}^{\mathrm{ini}} and the error is limited (see Eqs. (42) and (46)). For large deformation problems, the upper treatment is necessary.

The updating rule in large deformation NMM can be summarized as:

Update the coordinates of mesh.

Update mass density using Eq. (43) and then update stress using Eqs. (45) and (46).

Update other variables applying on the integration area if exists.

Verifications of large deformation NMM

Verification of the volume error and rotation angle error

The simple rotation problem illustrated in Fig. 1 is simulated by the proposed large deformation NMM (LDNMM). The time step is adopted as 0.01s and is corresponding to an incremental rotation angle of 3.6° in each step. The simulation results of NMM and LDNMM are compared in Table 4 and illustrated in Fig. 10. Obviously, the volume expansion existed in NMM (T1 and T2, referring to Fig. 10) is eliminated in the LDNMM (T3 and T4) where the volume error is negligible (Table 4).

As the rotation angle error is owing to numerical damping, it can be eliminated by adjust the parameter involved in Newmark method [³⁸]. As introduced in section 2.1, there are two parameters, as \mathit{\gamma } and \beta(see Eq. (5)), in Newmark method. All the tests in Table 4 adopt same g value of 0.5. As discussed in Ref [³⁸], \beta with value larger than 0.5 introduces numerical damping and a \mathit{\gamma }larger than \mathrm{0.25}(\beta +\mathrm{0.5}) can ensure the numerical stability. By changing \beta value from 1.0 to 0.5, the rotation angle errors existed in T1 and T3 are eliminated in T2 and T4 where the rotation angle error is limited.

Based on the results in Table 4 and Fig. 10, it can be concluded that: 1) the small deformation assumption will lead to significant volume errors, and this volume errors are eliminated in large deformation NMM; 2) numerical damping will lead to some rotation angle error, and a simple remedy is to use \beta =\mathrm{0.5} instead of 1.0.

In addition, the volume error and the rotation angle error in T4 remain the same after five circles, as volume error less than 1e−6% and rotation angle error about −0.033%. The errors are much smaller than that obtained by the original NMM (referring to the Table 3) using 10000 times steps (T4 involves only 100 time steps), where the volume error is 0.395% and the rotation angle error is −0.099%, i.e., LDNMM is superior to NMM for the rotation simulation.

Verification of the stress error

The stress waves in the simple rotation tests are plotted in Fig. 11. Referring to Fig. 11, the stress from large deformation NMM is of high accuracy. The max error of new scheme with 3.6°/step is 0.8% of the amplitude, close to 0.9%, which is the error of original NMM with extreme small incremental rotation angle 0.036°/step and removing high frequency wave. Besides, the troublesome rotation stress vibration phenomenon occurred in the original NMM was eliminated.

Verification of the convergence of large deformation analysis

The large deformation iteration in the new scheme is outside the open-close iteration. Since the open-close iteration used here is as same as that used in the original NMM and has been proved to be converged, only the convergence of large deformation iteration should be verified. A static analysis with large angle rotation is used to check the convergence of large deformation analysis.

As shown in Fig. 12, the deflection of cantilever beam ABCD with point forces P=\mathrm{500}\mathrm{N} applied at vertex C and D were simulated by the proposed large deformation NMM. Following are the settings: beam length L= 1000 mm, section size \mathrm{20}\mathrm{mm}\times \mathrm{100}\mathrm{mm}, elastic mesh with Young’s modulus E = 200 Mpa and Possion’s ratio \mu =\mathrm{0.3}.

Using the theories of small deflection beam, it is obtained that vertical displacement at the end is v=\mathrm{1000}\mathrm{mm}and horizontal displacement is u=\mathrm{0}. While using the theoretical solution of large deflection beam [⁴⁸], the displacement are v=\mathrm{603}\mathrm{mm}and u=-\mathrm{253}\mathrm{mm}. Clearly, the solution of small deflection beam is unreal and large deformation theory is necessary.

As shown in Fig. 13(a), the error of the horizontal and the vertical displacement are 1.9% and 0.8%, respectively. The iteration process in large deformation NMM is plotted in Fig. 13(b), where the max error of \theta decrease 1 order of magnitude about every 6 iterations and the vertical displacement is converged in 15 iterations. The convergence of large deformation iteration can be assured in several iterations, even for the case that large rotation occurs in one step.

Numerical examples

A simple rock fall simulation

Rock fall [⁴⁹] is a common disaster in geotechnical engineering and is lack of a good numerical simulation tool. A simple rock falling problem is simulated here to demonstrate the applicability of the proposed large deformation NMM.

As shown in Fig. 14, a square block with edge length of 0.05 m falls down to a slope. The slope is 0.5 m in height, with slope angle of 45° and connecting to horizontal ground at slope foot. The computation settings were: elastic plain strain with E=\mathrm{10}\mathrm{GPa}and \mu =\mathrm{0.3}, Gravity g=-\mathrm{10}\mathrm{m}/{\mathrm{s}}^{\mathrm{2}}, and default time step interval \mathrm{\Delta }{t}_{\mathrm{0}}= 0.01 s.

The block tracks from the original NMM and that from the proposed large deformation NMM were illustrated in Fig. 14. In LDNMM, the volume change of new scheme was −0.00024%, while that of original NMM was 62%. Besides, the two methods offer two different traces. In the simulation using the original NMM, the block starts to rotate after it contacted with the slope at point A (Fig. 14) and its volume will expand. Such rotation expansion after the first impact will cause a tiny error of the contact position and time during the open-close iteration in the original NMM. After the first impact, the block will keep rotating and expanding. Thus, the second impact between the block and slope happens earlier.

In short, the rotation errors in the original NMM could lead to the wrong solution of large deformation problems, such as impact, toppling, and rolling. A large deformation NMM is required to deal with such kind of simulation.

A simple sliding test on an arc surface

As shown in Fig. 15, a sliding test is simulated where a rectangle block is sliding along a frictionless, arc shaped slope (the arc is approximated by a polyline with 20 points). The block is supposed to enter the slope foot with an initial velocity of v, then move up to the highest point of the arc and finally return back to the horizontal ground at a speed of -v.

The NMM mesh is shown in Fig. 16. The gravity is set as g=-\mathrm{10}{\mathrm{m/s}}^2. The initial velocity v entering the arc should be v=\sqrt{\mathrm{2}gh}=\sqrt{\mathrm{20}}\mathrm{m/s} and is realized by imposing a horizontal acceleration \mathrm{4}\sqrt{\mathrm{20}}{\mathrm{m/s}}^{\mathrm{2}} from 0 to 0.25 s. The computation settings were: elastic plain strain with E=\mathrm{30}\mathrm{GPa} and \mu =\mathrm{0.3}, default time step interval \mathrm{\Delta }{t}_{\mathrm{0}}=\mathrm{0.005}\mathrm{s} and Newmark parameters [0.5, 1.0].

The block tracks simulated by LDNMM were illustrated in Fig. 17. It is found the block reaches the highest point at 0.921 s with the center elevation of 0.957 m (error −4.3%). The block volume V vs time t is plotted in Fig. 18, which indicates that: the volume in the original NMM is increasing step by step when block is sliding on the arc surface, and this unreasonable volume expansion is eliminated in the proposed LDNMM.

The recorded velocity versus time t is plotted in Fig. 19. Referring Fig. 19, if \beta =\mathrm{1}, the block speed after sliding back to the horizontal surface has lost 3.85% (the original NMM) and 4.05% (the proposed LDNMM), respectively. Such loss of velocity could be induced by the numerical damping.

During the sliding process, the energy of the sliding block is calculated by Eq. (42) as

(47)E=\frac{1}{2}{{mv}_c}^2+\frac{1}{2}{{I\omega }_c}^2+mg\mathrm{\Delta }h,

where m denotes the density, I is the rotational inertia, v_{\mathrm{c}} is the velocity of the block center, {\omega }_{\mathrm{c}} is the angular velocity of block center (which is measured using the velocities of two points on the block), and \mathrm{\Delta }h is the lifting height of the block.

The energy from different \beta is compared in Fig. 20. It can be noted that energy loss happens when the block is sliding on the arc surface. A small \beta can decrease such energy loss, however, may lead to unreasonable energy oscillation. Referring to Fig. 20, the energy oscillation is negligible when the test using a larger \beta such as 0.8 (i.e., introduce more numerical damping). A possible explanation can be that the numerical damping decreases the oscillation of the energy storage in contact spring which is not existed in nature. Another possible explanation is reported in Ref. [³⁹], as that the numerical damping allows contact force oscillation to dissipate rapidly and results in a stable contact state. Thus, although numerical damping will lead to energy loss during rotation process, the numerical damping is still necessary for the simulation involving contacts.

Conclusions

Through the theoretical and numerical analysis in this study, it can be revealed that volume expansion and stress vibration will happen in the rotation simulation by original NMM. Such errors are owing to the small deformation assumption adopted in the original NMM and can be eliminated by replacing the linear strain with Green strain.

A large deformation NMM is proposed by introducing Green strain and the large deformation iteration. In the proposed LDNMM, open-close iteration is embedded in the large deformation iteration and a updating strategy for variables defining on the physical body is included.

The proposed method was used to analyze block rotation, beam bending, and rock falling problems. The numerical test results demonstrate that the proposed LDNMM is proved to be converged, serves well for rotation problems, and can eliminate these rotation errors.

Numerical test results also indicate that the time step integration in NMM is in fact a special case of Newmark method which use b as 1.0. Numerical damping exists in NMM or LDNMM if the Newmark parameter b is larger than 0.5. Such numerical damping will lead to a decreasing of angular velocity; however, can help to avoid unreasonable oscillations of the energy in the simulation involving contacts. A further study for the influence of numerical damping on contacts is required.

In addition, although many efforts have been carried out to deal with crack issues over the past decades and great achievements have been achieved (e.g., Refs. [⁹,¹⁰,¹²,⁵⁰–⁵²], the issues of crack propagation induced (or accompanied with) large displacement and rotation, such as the bending crack, hydraulic fracturing and topping, etc., need urgent solutions. The coupling simulation of cracks, contacts and large deformation theory will be a promising work in in future.

References

[1] Shi G H. Manifold method. In: Proceeding of the 1st International Forum on DDA Simulation of Discontinuous Media(ICADD-1). New Mexico: TSI Press, 1996, 52–204

[2] Shi G H. Discontinuous deformation analysis: A new numerical model for the statics and dynamics of deformable block structures. Engineering Computations, 1992, 9(2): 157–168

[3] Ma G W, An X M, He L. The numerical manifold method: a review. International Journal of Computational Methods, 2010, 7(1): 1–32

[4] Babuska I, Melenk J M.The partition of unity method. International Journal for Numerical Methods in Engineering, 1997, 40(4): 727–758

[5] Chen G, Ohnishi Y, Ito T. Development of high-order manifold method. International Journal for Numerical Methods in Engineering, 1998, 43(4): 685–712

[6] Su H D, Xie X L, Liang Q Y. Automatic programming for high-order numerical manifold method. In: Proceeding of the 6th International Conference Analysis of Discontinuous Deformation (ICADD-6). Trondheim: A. A. Balkema publishers, 2003, 153–160

[7] Strouboulis T, Babuska I, Copps K. The design and analysis of the Generalized Finite Element Method. Computer Methods in Applied Mechanics and Engineering, 2000, 181(1–3): 43–69

[8] An X M, Li L X, Ma G W, Zhang H H. Prediction of rank deficiency in partition of unity-based methods with plane triangular or quadrilateral meshes. Computer Methods in Applied Mechanics and Engineering, 2011, 200(5–8): 665–674

[9] Zheng H, Xu D. New strategies for some issues of numerical manifold method in simulation of crack propagation. International Journal for Numerical Methods in Engineering, 2014, 97(13): 986–1010

[10] Zhang G X, Li X, Li H F. Simulation of hydraulic fracture utilizing numerical manifold method. Science China. Technological Sciences, 2015, 58(9): 1542–1557

[11] Yang Y, Tang X, Zheng H, Liu Q, He L. Three-dimensional fracture propagation with numerical manifold method. Engineering Analysis with Boundary Elements, 2016, 72: 65–77

[12] Ma G W, An X M, Zhang H H, Li L X. Modeling complex crack problems using the numerical manifold method. International Journal of Fracture, 2009, 156(1): 21–35

[13] Su H, Qi Y, Gong Y, Cui J H. Preliminary research of Numerical Manifold Method based on partly overlapping rectangular covers. In: DDA Commission of International Society for Rock Mechanics, Proceedings of the 11th International Conference on Analysis of Discontinuous Deformation (ICADD11). Fukuoka: CRC Press, 2013

[14] Liu Z, Zheng H. Two-dimensional numerical manifold method with multilayer covers. Science China. Technological Sciences, 2016, 59(4): 515–530

[15] Yang Y, Zheng H, Sivaselvan M V. A rigorous and unified mass lumping scheme for higher-order elements. Computer Methods in Applied Mechanics and Engineering, 2017, 319: 491–514

[16] An X, Ma G, Cai Y, Zhu H. A new way to treat material discontinuities in the numerical manifold method. Computer Methods in Applied Mechanics and Engineering, 2011, 200(47–48): 3296–3308

[17] Jiang Q, Zhou C, Li D. A three-dimensional numerical manifold method based on tetrahedral meshes. Computers & Structures, 2009, 87(13–14): 880–889

[18] He L, Ma G W. Development of 3D numerical manifold method. International Journal of Computational Methods, 2010, 7(1): 107–129

[19] Wu Y, Chen G, Jiang Z, Zhang L, Zhang H, Fan F, Han Z, Zou Z, Chang L, Li L. Research on fault cutting algorithm of the three-dimensional numerical manifold method. International Journal of Geomechanics, 2017, 17(5): E4016003

[20] Lin C T, Amadei B, Sture S, Jung J. Using an Augmented Lagrangian Method and Block Fracturing in the DDA Method (No. SAND-93-0817C; CONF-940506-3). Albuquerque, NM: Sandia National Labs, 1994

[21] Li X, Zheng H. Condensed form of complementarity formulation for discontinuous deformation analysis. Science China Technological Sciences, 2015, 58(9): 1509–1519

[22] Shi G. Contact theory. Science China Technological Sciences, 2015, 58(9): 1450–1496

[23] Lin S Z, Qi Y F, Su H D. Formulation of high-order numerical manifold method and fast simplex integration based on special matrix operations. In: Proceeding of the 7th International Conference Analysis of Discontinuous Deformation (ICADD-7). Honolulu: CRC Press, 2005

[24] Qu X L, Fu G Y, Ma G W. An explicit time integration scheme of numerical manifold method. Engineering Analysis with Boundary Elements, 2014, 48(6): 53–62

[25] Qu X L, Wang Y, Fu G Y, Ma G W. Efficiency and accuracy verification of the explicit numerical manifold method for dynamic problems. Rock Mechanics and Rock Engineering, 2015, 48(3): 1131–1142

[26] Koo C Y, Chern J C. Modification of the DDA method for rigid block Problems. Journal of Rock Mechanics and Mining Sciences, 1998, 35(6): 683–693

[27] MacLaughlin M M, Sitar N. Rigid body rotations in DDA. In: Proceedings of the 1st international forum on discontinuous deformation analysis (DDA) and simulation of discontinuous media. New Mexico: TSI Press, 1996

[28] Cheng Y M, Zhang Y H. Rigid body rotation and block internal discretization in DDA analysis. International Journal for Numerical and Analytical Methods in Geomechanics, 2000, 24(6): 567–578

[29] Fan H, Zheng H, He S. S-R decomposition based numerical manifold method. Computer Methods in Applied Mechanics and Engineering, 2016, 304: 452–478

[30] Fan H, Zheng H, Zhao J. Discontinuous deformation analysis based on strain-rotation decomposition. International Journal of Rock Mechanics and Mining Sciences, 2017, 92: 19–29

[31] Wei W, Jiang Q. A modified numerical manifold method for simulation of finite deformation problem. Applied Mathematical Modelling, 2017, 48: 673–687

[32] Bathe K J, Ramm E, Wilson E L. Finite element formulations for large deformation dynamic analysis. International Journal for Numerical Methods in Engineering, 1975, 9(2): 353–386

[33] Bathe K J. Finite Element Procedures. New Jersey: Prentice Hall, 1996

[34] Terada K, Asai M, Yamagishi M. Finite cover method for linear and non-linear analyses of heterogeneous solids. International Journal for Numerical Methods in Engineering, 2003, 58(9): 1321–1346

[35] Hansbo A, Hansbo P. A finite element method for the simulation of strong and weak discontinuities in solid mechanics. Computer Methods in Applied Mechanics and Engineering, 2004, 193(33–35): 3523–3540

[36] Rabczuk T, Zi G, Gerstenberger A, Wall W A. A new crack tip element for the phantom-node method with arbitrary cohesive cracks. International Journal for Numerical Methods in Engineering, 2008, 75(5): 577–599

[37] Chau-Dinh T, Zi G, Lee P S, Rabczuk T, Song J H. Phantom-node method for shell models with arbitrary cracks. Computers & Structures, 2012, 92–93(3): 242–256

[38] Wang C Y, Chuang C, Sheng J. Time integration theories for the DDA method with finite element meshes. In: Proceedings of the 1st International Forum on Discontinuous Deformation Analysis (DDA) and Simulations of Discontinuous Media. Albuquerque, NM: TSI Press, 1996

[39] Doolin D M, Sitar N. Time integration in discontinuous deformation analysis. Journal of Engineering Mechanics, 2004, 130(3): 249–258

[40] Jiang Q, Chen Y, Zhou C, Yeung M R. Kinetic energy dissipation and convergence criterion of discontinuous deformations analysis (DDA) for geotechnical engineering. Rock Mechanics and Rock Engineering, 2013, 46(6): 1443–1460

[41] Tang C A, Tang S B, Gong B, Bai H M. Discontinuous deformation and displacement analysis: From continuous to discontinuous. Science China. Technological Sciences, 2015, 58(9): 1567–1574

[42] Antoci C, Gallati M, Sibilla S. Numerical simulation of fluid-structure interaction by SPH. Computers & Structures, 2007, 85(11–14): 879–890

[43] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A simple and robust three-dimensional cracking-particle method without enrichment. Computer Methods in Applied Mechanics and Engineering, 2010, 199(37–40): 2437–2455

[44] Ren H, Zhuang X, Cai Y, Rabczuk T. Dual-horizon peridynamics. International Journal for Numerical Methods in Engineering, 2016, 108(12): 1451–1476

[45] Ren H, Zhuang X, Rabczuk T. Dual-horizon peridynamics: A stable solution to varying horizons. Computer Methods in Applied Mechanics and Engineering, 2017, 318: 762–782

[46] Zienkiewicz O C, Taylor R L, Fox D D. The finite element method for Solid & Structural Mechanics. Singapore: Elsevier, 2015

[47] Jiang W, Zheng H. An efficient remedy for the false volume expansion of DDA when simulating large rotation. Computers and Geotechnics, 2015, 70: 18–23

[48] Bisshopp K E, Drucker D C. Large deflection of cantilever beams. Mathematics and Fly Fishing, 1945, 3(3): 272–275

[49] Macciotta R, Martin C D, Cruden D M, Hendry M, Edwards T. Rock fall hazard control along a section of railway based on quantified risk. Assessment and Management of Risk for Engineered Systems and Geohazards, 2017, 11(3): 272–284

[50] Rabczuk T, Belytschko T. Cracking particles: A simplified meshfree method for arbitrary evolving cracks. International Journal for Numerical Methods in Engineering, 2004, 61(13): 2316–2343

[51] Rabczuk T, Belytschko T. A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Computer Methods in Applied Mechanics and Engineering, 2007, 196(29–30): 2777–2799

[52] Rabczuk T, Zi G, Gerstenberger A, Wall W A. A new crack tip element for the phantom-node method with arbitrary cohesive cracks. International Journal for Numerical Methods in Engineering, 2008, 75(5): 577–599

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature