Complementarity problem arising from static growth of multiple cracks and MLS-based numerical manifold method

Available online at www.sciencedirect.com

ScienceDirect Comput. Methods Appl. Mech. Engrg. 295 (2015) 150–171 www.elsevier.com/locate/cma

Complementarity problem arising from static growth of multiple cracks and MLS-based numerical manifold method Hong Zheng a,∗ , Feng Liu b , Xiuli Du a a Key Laboratory of Urban Security and Disaster Engineering, Ministry of Education, Beijing University of Technology, Beijing 100124, China b State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences,

Wuhan 430071, China Received 23 January 2015; received in revised form 4 June 2015; accepted 6 July 2015 Available online 10 July 2015

Highlights • • • • •

The growth of multiple cracks can be reduced to a nonlinear complementarity problem (NCP). The NCP can be solved by the projection-contraction algorithm PCA. A consensus in fracture mechanics on brittle fracture is subverted. Some interesting and profound phenomena in brittle fracture are revealed. The collapse load of a brittle structure is unique and independent of the cracking paths.

Abstract The full response of a brittle structure containing multiple cracks to loading under the servo control is of vital importance in the evaluation of properties of the structure. During crack growth, the fracture toughness condition at all the crack tips as well as the equilibrium condition should be obeyed, leading to a nonlinear complementarity problem (NCP). The vector-valued function in the NCP depends implicitly on the cracking increments which in turn determine the stress field. The stress field can be obtained through solving a mixed variational problem. Since the degrees of freedom in the discrete variational problem vary with cracking not only in magnitude but also in number, the Jacobian matrix of the NCP is hard to compute and, it cannot be expected to be solved by the Newton methods that involve the calculation of Jacobian matrices. Therefore, a well-scaled projection-contraction algorithm is designed. The proposed procedure is able to simulate growth of multiple cracks in a natural way, allowing the crack tips to stop anywhere with no sensitivity to node configuration or cracking increments. Through the analysis of some examples that have been widely tested, many interesting and profound phenomena are found which have never been revealed. c 2015 Elsevier B.V. All rights reserved. ⃝

Keywords: Brittle fracture; Crack initiation and growth; Nonlinear complementarity problem; Numerical manifold method; Moving least squares interpolation

∗ Corresponding author. Tel.: +86 27 8719 9226; fax: +86 27 8719 7386.

E-mail address: [email protected] (H. Zheng). http://dx.doi.org/10.1016/j.cma.2015.07.001 c 2015 Elsevier B.V. All rights reserved. 0045-7825/⃝

H. Zheng et al. / Comput. Methods Appl. Mech. Engrg. 295 (2015) 150–171

151

1. Introduction Failures of brittle structures are mainly caused by crack initiation and propagation. In order to simulate crack propagation more efficiently, several sophisticated numerical methods have been developed over the past two decades, among which are the element-free Galerkin method (EFGM, [1]), the partition of unity method (PUM, [2]), the extended finite element method (XFEM, [3]), the generalized finite element method (GFEM, [4]), and many others. Compared with all these methods, the numerical manifold method (NMM, [5]) has its own advantages and features, although it also falls into the category of the partition of unity. Investigations and applications of NMM are becoming increasingly active. For example, only in the first ten months of the last year had the number of international journal papers on NMM exceeded forty, excluding those on-line. NMM is easy to implement the p-adaptive analysis by using higher degree polynomials as local approximations, see [6,7], for example. Similar to GFEM, however, it might incur the issue of linear dependency. An et al. proposed an algorithm for predicting the rank deficiency of the stiffness matrix [8]. Recently, Ghasemzadeh et al. developed a scheme with high order polynomials as the local approximations that has no linear dependency [9]. Tian and Wen proposed a more effective procedure that enjoys high precision but gets rid of linear dependency [10]. Although formulated in the XFEM setting, it is applicable to NMM as well. By adding the feature of excavation simulation into the NMM code, Tal et al. modeled the stability of underground openings embedded in discontinuous rock masses under high in-situ stress conditions [11]. With incorporation of a rockbolt model into the NMM code, Wu and Wong simulated the reinforcement of tunnels [12] where rockfall might poses a great hazard. Zheng et al. proposed a method coupling NMM and the graph theory allowing for stability analysis of both rock and soil slopes within the same framework [13]. Wu and Fan incorporated the firstorder Higdon absorbing boundary into the numerical manifold method (NMM) to reduce reflections from artificial boundaries induced by truncating infinite media [14]. With no need to introduce any micro-mechanic parameters, NMM is capable of treating initiation and propagation of multiple cracks in a natural way. So far, most applications of NMM just focus on this field. Great efforts have been made to investigate the initiation and propagation of cracks in geomaterials or concretes, see [15,16] for example, where the propagation of cracks is due to tension and shear. Terada et al. presented the finite cover method [17], alias of NMM, for progressive failure analyses with cohesive fracture zone in heterogeneous solids and structures, where both the fracture energy and the tensile strength of the material have to be given. Zhang and Ma conducted the analysis of very complicated cracks in functionally graded materials [18]. By adopting the Mohr–Coulomb criterion with a tensile cut-off, progressive failures of rock slopes are investigated in [19–21], where frictional contact of crack surfaces is involved. NMM always advocates the mathematical covers are generated from uniform meshes in which all elements are congruent. In the presence of singularities such as crack tips, Yang et al. developed a procedure for local refinement [22]. Zhang et al. solved thermo-mechanical fracture problems in two-dimensional setting [23], where the singularity of temperature as well as displacement is reflected. Recently, another role of NMM is found in the analysis of thin plate bending problems [24]: NMM can be used to rescue those incompatible elements that developed in the finite element history, which require very stringent mesh configuration to pass the patch test condition. In addition to the application to solids and structures, NMM is also used to conduct seepage analyses, see [25–27] for details. Most simulations of crack propagation based on NMM are conducted in the setting of static loading. If the cracking criterion is based on fracture mechanics, the static loading condition requires that the fracture toughness condition as well as the equilibrium condition be obeyed. If the cracking criterion is controlled by the material strength, the strength criterion instead of the fracture toughness condition should be satisfied. During crack growth, however, it is assumed in most applications of NMM and other PU-based methods such as [28–30] that a crack tip always resides at an element boundary. Such a simplification would violate the fracture toughness condition or strength criterion, and lead more or less to the mesh dependency, which goes against the grain of NMM and XFEM. In this study, it is required that at any crack tip K eq ≤ K IC , where K eq represents the equivalent stress intensity factor determined by the maximum hoop stress criterion, and K IC the fracture toughness. The paper is organized as follows.

152

H. Zheng et al. / Comput. Methods Appl. Mech. Engrg. 295 (2015) 150–171

Section 2 encapsulates nonlinear complementarity problems to be involved in this study. In Section 3, a brief introduction to NMM is given. Section 4 discusses the fracture criterion for static cracking. In Section 5, the complementarity problem arising from the simulation of static growth of multiple cracks is formulated, in which the vector-valued function is calculated only after solving a mixed variational problem fitted to the servo control way. Since the complementarity problem cannot be solved by the Newton methods, in Section 6 a projection-contract algorithm tailored for the problem is developed. Some typical numerical tests are conducted in Section 7, where some interesting and profound phenomena are revealed, which have never been reported in the literature. Section 8 concludes this study. 2. Nonlinear complementarity problems Through investigations and applications over five decades, a huge volume of theories and algorithms on the nonlinear complementarity problems (NCPs) has been accumulated, which can be found in the monograph by Facchinei and Pang [31]. Here, we only touch upon some directly relative to this study. n → R n , the nonlinear complementarity problem associated with the Definition 1. Given a nonlinear mapping F: R+ n satisfying vector-valued function F, abbreviated as NCP(F), is to find a vector x ∈ R+

0 ≤ x ⊥ F (x) ≥ 0.

(1) n R+

Rn

Rn

Here, represents the n-dimensional Euclidean space; is the subspace of that is composed of vectors with nonnegative components. By x ≥ 0 we mean all components of x is nonnegative. The orthogonality condition of x and F (x) suggests that xT F (x) = 0, or xi Fi (x) = 0, ∀i = 1, . . . , n, in terms of the component-wise products, which is justified because x and F (x) are both nonnegative vectors. In applications, many problems can be reduced to the following mixed complementarity problem. m into R n and R m , respectively. The mixed complementarity Definition 2. Let G and H be two mappings from R n × R+ m such that problem, abbreviated as MiCP (G, H), is to find a pair of vectors (u, v) belonging to R n × R+  G (u, v) = 0 (2) 0 ≤ v ⊥ H (u, v) ≥ 0.

Here, vector u is called the free variable, and v the complementarity variable. Usually, system G (u, v) = 0 represents the equilibrium condition; while the complementarity relationship in system (2) represents the constraints [32,33]. 3. Brief introduction to NMM With the aim to solve in a unified way geotechnical problems from continua to discontinua, NMM introduces two cover systems, namely, the mathematical cover and the physical cover. In this section, we only give a brief introduction to NMM, more details can be found in [34]. The mathematical cover is constituted of simply-connected domains, denoted by Mi , i = 1, . . . , M; each Mi is referred to as a mathematical patch, and contains no geometric feature of the problem domain Ω . Here, M is the number of the mathematical patches. The union of all Mi ’s is required to cover Ω . The mathematical cover determines the solution precision, and can be deployed as the user wishes. Associated with each Mi is a weight function wi (x, y) continuous all over Ω . All wi (x, y)’s constitute the partition of unity subordinate to the mathematical cover {Mi }, satisfying wi (x, y) ≤ 1,

∀ (x, y) ∈ Mi ;

(3.1)

wi (x, y) = 0,

∀ (x, y) ̸∈ Mi ;

(3.2)

M  i=1

wi (x, y) = 1,

∀ (x, y) ∈ Ω .

(3.3)

H. Zheng et al. / Comput. Methods Appl. Mech. Engrg. 295 (2015) 150–171

153

Fig. 1. Problem domain (thick lines) and mathematical cover (fine lines).

a b c

Fig. 2. Physical patches. (a) Physical patches P11 and P12 generated from mathematical patch M1 ; (b) Physical patch P2 generated from mathematical patch M2 ; and (c) Physical patch P3 generated from mathematical patch M3 .

It is remarked that the continuity of wi (x, y) over Ω implies wi (x, y) vanishes on the boundary of Mi . Illustrated in Fig. 1 is such an example, where three mathematical patches, M1 —a big circular disc, M2 —a rectangle, and M3 —a small circular disc, cover the problem domain Ω containing the crack Γ . By cutting the mathematical patches one after another with the components of the problem domain, including the domain boundaries, the material interfaces, and the discontinuities, the physical patches are generated, as shown in j Fig. 2, where Pi represents the jth physical patch generated from mathematical patch Mi , j = 1, . . . , m i , with m i being the number of the physical patches generated from the same mathematical patch Mi . No superscript, such as P3 j in Fig. 2(c), indicates only one physical patch is generated from M3 , etc. All Pi ’s constitute the physical cover, which covers the problem domain exactly, namely, mi M   i

j

Pi = Ω .

j=1 j

Each Pi contains the local geometric features of the problem domain and can be assigned other given information. For example, P2 contains the crack tip of crack Γ (see Fig. 2(b)), and P3 is acted by a concentrated force f at point A (see Fig. 2(c)). j By exploiting the geometric and mechanical features of Pi , in principle one can construct a good enough local j j approximation to the solution over Pi , denoted by ui (x, y), which can always be expressed as j

j

j

ui (x, y) = Ti (x, y) di ,

j

(x, y) ∈ Pi .

(4)

j

j

Here, vector di is composed of the degrees of freedom for physical patch Pi ; for two-dimensional elasticity problems, j j j j the dimension of di is usually an even number, denoted by n i . The 2 × n i matrix Ti (x, y) is composed of some given j functions, which are able to reflect the local behaviors of the solution over Pi , see [34] for details. j

j

Each physical patch Pi inherits a weight function wi (x, y) from the corresponding mathematical patch Mi , all  j j wi (x, y) form the partition of unity subordinate to the physical cover Pi . Then, the addition of all the local approximations multiplied by the weight functions over the physical patches gives the global approximation to the

154

H. Zheng et al. / Comput. Methods Appl. Mech. Engrg. 295 (2015) 150–171

solution, namely, u (x, y) =

mi M  

j

j

wi (x, y) ui (x, y) ,

for (x, y) ∈ Ω .

(5)

i=1 j=1

To this point, it is clear that it is the cutting operation to the mathematical cover that makes possible the reproduction of weak and strong discontinuity in solution, and that it is because the local approximations over the physical patches can be constructed independently that a better global solution can be made. See [34] for details. j By numbering all the physical patches with single indices, say PI , etc., I = 1, . . . , P, instead of Pi , where P is the number of all physical patches, P=

M 

mi ,

i=1

and by substituting Eq. (4) into Eq. (5), Eq. (5) turns into u (x, y) = Nd.

(6)

Here, vector d consists of the components of degrees of freedom vector d I of all physical patches with a dimension of n as n=

P 

nI ,

I =1

with n I equal to the dimension of degrees of freedom vector d I of physical patch PI ; N in Eq. (6) is the 2 × n matrix, expressed by   N = N1 , · · · , N P , (7) where N I , I = 1, . . . , P, is known as the shape function matrix of physical patch PI , N I (x, y) = w I (x, y) T I (x, y) ,

(8)

with the dimension of 2 × n I . In principle, the mathematical cover can be constructed arbitrarily as long as it covers the problem domain. There are systematic schemes, see [35] for example, to generate the partition of unity subordinate to the mathematical cover, which can be made arbitrarily smooth. To date, however, a majority of applications of NMM have been selecting the finite element mesh as the mathematical cover, where a node corresponds to a mathematical patch that is composed of all elements linked at the node. During crack growth, the physical cover at the end of each load step will serve as the mathematical cover of the next load step, which in turn is cut by the cracks to be newly grown. As a result, the simulation of multiple crack growth needs to cut frequently the mathematical cover. If the mathematical cover is generated from the finite element mesh, sophisticated algorithms are necessary for very complicated cutting operations, which have been discussed in [19,36]. Instead, if the mathematical patches are specified as the influence domains of nodes scattered over the problem domain, such as the nodes in the MLS interpolation [34], cutting operations on the mathematical patches can be considerably simplified. In summary, the MLS-based NMM combines the advantages of both NMM and EFGM. 4. Fracture criterion for static cracking On the basis of energy dissipation, two types of fracture are distinguished: brittle fracture and ductile fracture. Ductile fracture is characterized by high energy dissipation due to plastic deformation and slip in the neighborhood of crack tips. This high energy dissipation is absent in brittle fracture. For this type of fracture to be investigated in this study, almost all the dissipated energy is due to the formation of new crack surfaces. On the basis of loading speed, two types of crack propagation can be classified: dynamic cracking and static cracking. For the type of dynamic cracking, a crack grows under the condition that the equivalent stress intensity

H. Zheng et al. / Comput. Methods Appl. Mech. Engrg. 295 (2015) 150–171

155

factor of the crack tip, denoted by K eq and to be explained subsequently, is greater than fracture toughness of the material, denoted by K IC . Belytschko et al. suggested to take advantage of the level set method to obtain the new position of the crack tip [37]. The key assumption is that the hyperbolicity indicator vanishes at the crack tip, leading to a hyperbolic partial differential equation with regard to the cracking speed. Due to the hyperbolic character of the partial differential equation, the standard Galerkin formulation is apt to incur instabilities and a stabilization technique, such as upwinding, Petrov–Galerkin, Galerkin least squares, streamline upwind Petrov–Galerkin, or finite increment calculus stabilization, is needed, as pointed out by Zhuang et al. in [38]. For the type of static cracking, which can be implemented only under the servo control, a crack grows under the condition of K eq = K IC ,

(9)

with K eq being the equivalent stress intensity factor at the crack tip. If the crack tip keeps still, then K eq < K IC .

(10)

The status of K eq > K IC is not existent for static cracking. In this study, the equivalent stress intensity factor K eq is determined by the maximum hoop stress criterion    √ 3 2 2 4 2K I I K I + 3 K I + 8K I I K eq =  3/2 ,  K I2

+ 12K I2I

− KI

K I2

(11)

+ 8K I2I

where K I and K I I represent the stress intensity factors of mode-I and mode-II, respectively. If Eq. (9) is satisfied at a crack tip, the crack will grow in the direction given by the argument  K I − K I2 + 8K I2I θC = 2 tan−1 4K I I

(12)

with the crack tip as the origin of the polar system and the crack line on the line of θ = −π. 5. Complementarity problem for static growth of multiple cracks In this section, we first formulate the NCP induced by the requirement that the fracture toughness conditions of all the crack tips be satisfied. Since the function in the NCP depends on the stress field, then we formulate the boundary value problem fitted to the static growth of cracks. 5.1. NCP induced by fracture toughness conditions Now we assume that at time t0 there are m crack tips preexistent in the 2-dimensional problem domain Ω0 in equilibrium, as shown in Fig. 3(a). At time t, a moment closely subsequent to time t0 , the m crack tips grow statically by length li , i = 1, . . . , m, and Ω0 turns into Ω , as illustrated in Fig. 3(b). According to the fracture criterion i at the ith crack tip has the discussed in the previous section, the cracking increment li and the quantity K IC − K eq complementarity relationship,   i i ≥ 0, li K IC − K eq = 0, li ≥ 0, K IC − K eq i denotes the equivalent stress intensity factor at the ith crack tip, see Eq. (11) for its for i = 1, . . . , m; where K eq definition. By assembling the complementarity relationships of all the m crack tips and then scaling, we have

0 ≤ a ⊥ H (a) ≥ 0,

(13)

namely, NCP(H). Here, vector a is called the cracking vector defined by a = (a1 , . . . , am )T ,

(14)

156

H. Zheng et al. / Comput. Methods Appl. Mech. Engrg. 295 (2015) 150–171

a

b

Fig. 3. Problem domain at time t0 (a), and time t (b), where l3 is the cracking increment of the 3rd crack tip, etc.

with li , (15) l where l is a length specified by the user to control the biggest cracking increment; and the vector-valued function, ai =

H (a) = (H1 (a) , . . . , Hm (a))T ,

(16)

with Hi (a) the ith component function Hi (a) = 1 −

i (a) K eq

K IC

,

(17)

i = 1, . . . , m. As will be seen shortly, it is justified to regard function H (a) as an implicit function of the cracking vector a. Both a and H range in the interval of [0, 1]. Hence, NCP(H) is well-scaled. 5.2. Boundary value problem fitted to the servo control way In order that the cracks in Ω statically grow, the servo control must be applied to Ω , either in a load-controlled way or in a displacement-controlled way. In order to obtain the complete stress–strain curve of a rock sample, as is well known, the testing machine must adopt the servo control, by which the axial load acting on the sample can be increased or decreased such that the sample is always in equilibrium to avoid collapse. There are two control ways: the load-controlled way and the displacement-controlled way. By a load-controlled way, we mean that the traction boundary S p of Ω is loaded by the servo control, while the displacement boundary Su is fixed, suggesting the boundary conditions ∂n σ = λl pc ,

on S p ;

(18)

and u = 0,

on Su ,

(19)

hold during the loading process. Demonstrated in Fig. 3(b) is a load-controlled way. Here, σ is the stress vector with components σx , σ y , and τx y ; ∂n represents the matrix   cos α 0 sin α ∂n = , (20) 0 sin α cos α with α the angle of the outward unit normal n of S p ; pc is the given traction vector; λl is the load factor to be determined under the condition that both the equilibrium and fracture toughness conditions are satisfied. u denotes the displacement vector with components u (horizontal) and v (vertical).

H. Zheng et al. / Comput. Methods Appl. Mech. Engrg. 295 (2015) 150–171

157

By a displacement-controlled way, we mean that the displacement boundary Su of Ω is servo-controlled, while the traction boundary S p is assumed to be free for the simplicity of presentation, indicating the boundary conditions ∂n σ = 0,

on S p ,

(21)

u = λd u, ¯

on Su ,

(22)

and

hold during the process of loading. Here, u¯ is the given displacement vector; λd is the displacement factor to be determined according to the same principle as the load-controlled way. In order to obtain the stress field, one can solve the mixed variational problem    T T  δu p = λl δuT pc   δε σ + Ω Su Sp  (23)  T   δp (u − λd u) ¯ = 0, Su

with regard to the unknown variables u and p, for the problem domain Ω at time t, at which the cracking vector a is assumed known. The load-controlled way corresponds to the condition that λl is variable but λd = 0; while the displacement-controlled way corresponds to the condition that λd is variable but λl = 0. Here, δu and δε the virtual displacement and strain, respectively; the vector p represents the unknown traction vector acting on Su . System (23) is derived from the Lagrange multiplier method that converts the essential boundary condition on Su to the natural boundary condition, where the Lagrange multiple vector is identified as the traction vector p on Su . Besides, it should be emphasized that system (23) corresponds to the problem domain Ω with the given cracking vector a, where the fundamental variables u and p are both total quantities. This is justified under the condition that all the cracks are subject to tension and shear during propagation. For the growth of cracks in compression and shear, the problem must be formulated in the incremental form, which is so complicated and to be reported in the future. To this point, it is justified to regard Hi (a) in Eq. (17) and accordingly H (a) in Eq. (16) as implicit functions of vector a: if a is given, the stress field σ can be obtained through solving the variational principle (23), and Hi (a) is hence obtained according to Eqs. (11) and (17). 6. Projection-contraction algorithm (PCA) By the MLS-based NMM [34], we can have the discrete form of problem (23) as follows:      K C d λl qc = , λd qu CT 0 pu

(24)

which is a KKT system and can be solved without difficulties by some specific algorithms, such as [39]. For the load-controlled way, λl is to be determined while λd is equal to zero. For the displacement-controlled way, λd is to be determined while λl is equal to zero. In Eq. (24), qc is the load vector equivalent to the traction vector pc on S p ,  N T pc ; (25) qc = Sp

pu is the Lagrange multiplier vector consisting of traction components of the nodes deployed along boundary Su , through which the traction vector, p = N p pu ,

(26)

on Su , can be interpolated; where N p is the shape function matrix consisting of the linear Lagrange interpolations to the traction vector p on Su . qu is the load vector corresponding to u, ¯  qu = N Tp u; ¯ (27) Su

158

H. Zheng et al. / Comput. Methods Appl. Mech. Engrg. 295 (2015) 150–171

C is the constraint matrix defined as  C= NT N p ;

(28)

Su

K is the stiffness matrix  K= BT DB,

(29)

Ω

and matrix B is defined by B = [B1 , . . . , B P ] ,

(30)

B I = ∂N I ,

(31)

with

I = 1, . . . , P, and  ∂ ∂x   ∂= 0   ∂ ∂y

0



 ∂   . ∂y  ∂  ∂x

(32)

Now for some given factor λl (load-control), or λd (displacement-control), we have the unknowns of a, d and pu . According to Definition 2, system (24) and complementarity relationship (13) collectively constitute a mixed complementarity problem, with d and pu as the free variables and a as the complementarity variable. In general, the most efficient algorithms for NCPs are non-smooth Newton methods, for example those proposed in [31]. But none of them can be expected to solve the mixed complementarity problem posed here, because not only vector d but also its dimension n depend implicitly on vector a, leading to huge difficulties in computing the Jacobian of mapping H (a), which is involved in the Newton methods. Fortunately, the difficulties can be overcome by fixed-point algorithms, which do not need Jacobians. In this study, we select the projection-contraction algorithm, abbreviated as PCA, to solve NCP(H) defined in Eq. (13). PCA is suitable to solving NCPs, with no need to compute the Jacobian [40]. Let us suppose at time t0 we have the cracked problem domain Ω0 in equilibrium as shown in Fig. 3(a). Now we find the new configuration Ω at time t under the load level λl pc on S p (load-control), as illustrated in Fig. 3(b); or the given displacement λd u¯ onto Su (displacement-control). 6.1. Determination of load factor To begin with, we evaluate a proper load factor λl such that the cracks grow a proper length at time t. Here, the load-controlled way is assumed with λd = 0 in system (24). For the displacement-controlled way, the treatment is the same, with λl replaced by λd and λl = 0. j j A crack tip at time t0 , say the jth crack tip, has either the status of K eq (t0 ) < K IC or the status of K eq (t0 ) ≈ K IC . Here, “≈” is in the sense of relaxation,    j  (33) K eq (t0 ) − K IC  < ε K K IC , with ε K being a user specified tolerance. For the convenience of presentation below, let us denote by j1 , . . . , jk the indices of those crack tips to be grown, j at which K eq (t0 ) ≈ K IC , j = j1 , . . . , jk . j Now, let the crack tips of indices j1 , . . . , jk grow by length l, in the direction of θC determined by Eq. (12). In   response, we have the cracking vector a0 = ai0 with ai0 = 1 for i = j1 , . . . , jk and ai = 0 for other i’s. Fig. 4 illustrates such a case, where only the second crack tip grows a length of l.

H. Zheng et al. / Comput. Methods Appl. Mech. Engrg. 295 (2015) 150–171

159

After the physical cover corresponding to a0 is obtained, by solving system (24) with λl = 1, we have the stress j field σ |λl =1 , which is used to calculate the equivalent stress intensity factors K eq |λl =1 , j = j1 , . . . , jk . Then, λl is determined by λl =

K IC , K eq |λl =1

(34)

with   j K eq |λl =1 = max K eq |λl =1 , j = j1 , . . . , jk .

(35)

The above evaluation of λl is based on the fact that the equivalent stress intensity factor of a crack tip is proportional to the external load level. As a result, under the action of λl pc on S p and the given cracking vector a0 , the maximum equivalent stress intensity of all the crack tips equals K IC . Subsequently, we will find the cracking vector a under the load factor λl given by Eq. (34). 6.2. Projection-contraction algorithm for cracking vector Observing the fact that for a given λl (or λd ) the cracking vector a uniquely determines the stress field σ and accordingly Hi (a), we design an algorithm called HA for calculating H (a) defined in Eqs. (16) and (17), which is repeatedly invoked from PCA as follows. Algorithm HA Input: the cracking vector a. Output: the vector-valued function H (a). for i = 1 to m //calculate the growth direction of cracks to be growing if ai > 0 then //the crack tip is growing  θCi = 2 tan−1

K Ii −



K Ii

2

 2 +8 K Ii I

4K Ii I

;

//based on stress field at time t0

let the ith crack tip grow a length of lai in direction of θCi ; end if end (for-i) generate the physical cover for the cracking vector a and the cracking direction θCi ; solve system (24) for displacement vector d; calculate the stress field σ ; for i = 1 to m //calculate the vector-valued function H   i = K eq

√  3  2  2 4 2 K Ii I K Ii +3 K Ii +8 K Ii I





K Ii

Hi = 1 − end (for-i) 

2

   2  2 2 +12 K Ii I −K Ii K Ii +8 K Ii I

i K eq K IC ;

3/2 ;

// by Eq. (11)

// by Eq. (17)

Now NCP(H) defined in Eq. (13) can be solved by the following algorithm PCA. The algorithm is designed according to He’s algorithm in [40], which in turn is based on the proposition that vector a solves NCP(H) if and only if a solves the equation a = max [a − βH (a) , 0] . Here and in PCA, β is a positive for scaling H (a) and, the function “max” is vector-valued; for example, suppose the 3-dimensional vector xT = (1.23, 0, −2.16), then   1.23 max (x, 0) =  0  . 0

160

H. Zheng et al. / Comput. Methods Appl. Mech. Engrg. 295 (2015) 150–171

Fig. 4. Let the crack tip with K eq ≈ K IC grow a length of l under load level pc .

Algorithm PCA Input: an initial value of cracking vector a, tolerance ε, with the control parameters ξ , η ∈ (0, 1) and γ ∈ [1, 2). Output: the cracking vector a satisfying the conditions of both equilibrium and fracture toughness. Let β = 1 and k = 0; // counter k do H = HA (a); // call HA to calculate H(a) a˜ = max [a − βH, 0]; // a˜ is the array with the same length as a d (a, a˜ ) = a − a˜ ; // d (a, a˜ ) is the array with the same length as a if ∥d (a, a˜ )∥∞ ≤ ε then break; ˜ = HA (˜a); H ˜ F (a, a˜ ) = H − H; // F (a, a˜ ) is the array with the same length as a β∥F(a,˜a)∥2 r = ∥d(a,˜a)∥ ; 2 while r > ξ   ˜ = HA (˜a); β := 0.7∗ β ∗ min 1, 1 ; a˜ = max [a − βH, 0] ; H r

˜ d (a, a˜ ) = a − a˜ ; r = β∥F(a,˜a)∥2 ; F (a, a˜ ) = H − H; ∥d(a,˜a)∥2 end (while) D (a, a˜ ) = d (a, a˜ ) − βF (a, a˜ ); // array D (a, a˜ ) has the same length as a d(a,˜a)]T D(a,˜a) α=[ ; a := a − γ αd (a, a˜ ); ∥D(a,˜a)∥22

if r ≤ η then β := 1.5∗ β; k := k + 1; while k < K a . // K a = the number of maximum allowable iterations



In the PCA, F (a, a˜ ), D (a, a˜ ) and d (a, a˜ ) represent three temporary arrays determined by arrays a and a˜ . We have to admit that in the presence of a great number of cracks, the PCA will become inefficient. Therefore, a high efficient solver for NCPs is of vital importance for the ability of the proposed procedure to handle a large number of cracks. 6.3. Overall description of the solution of a load step To this point, we can give an overall description of the solution of a load step. At the start of the load step, we have a known stress field and a crack distribution. Our task is to find new cracking increments and the corresponding stress field with a proper load factor λl (or λd ), such that the fracture toughness condition as well as the equilibrium condition is satisfied. The solution of this load step starts by the determination of λl (or λd ) described in Section 6.1, followed by an i < K ; a0 = 1, if K i ≈ K . Then PCA is carried out with initial estimation of cracking vector a0 : ai0 = 0, if K eq IC i IC eq 0 the initial cracking vector a as input.

161

H. Zheng et al. / Comput. Methods Appl. Mech. Engrg. 295 (2015) 150–171

(a) Image of function Wi .

(b) Image of function Wi,x .

Fig. 5. Images of MLS-based weight function and its partial derivative with regard to x.

7. Illustration of numerical tests For all the examples in this study, the mathematical covers are composed of the influence domains of MLS-nodes, with the influence domain of each MLS-node serving as a mathematical patch. The local approximations to the displacement vector are the same as [34]. All the MLS-nodes are equally set parallel to the x- and y-directions, respectively. In order to decrease cutting operations involved in the crack growth, the influence domains of all the MLS-nodes are specified as the squares centered the nodes. In generating the MLS-based shape function wi (x, y) of node (xi , yi ), which also serves as the NMM weight function wi (x, y) of square patch Mi , the MLS-based weight function Wi (x, y) is taken as   Wi (x, y) = W (r x ) W r y , (36) where rx =

|x − xi | , L ix

(37.1)

ry =

|y − yi | , y Li

(37.2)

and

y

with L ix and L i the half lengths of Mi Eq. (36) is defined as  2   − 4r 2 + 4r 3 ,   3 4 W (r ) = 4 − 4r + 4r 2 − r 3 ,   3 3    0,

in the x- and y-directions, respectively; and the univariate function W (r ) in

r ≤ 0.5 (38)

0.5 < r ≤ 1 r > 1. y

Shown in Fig. 5 are the images of function Wi and its partial derivative Wi,x . For all the examples, L ix = L i = L i . Different from the conventional EFGM, the MLS-based NMM allows some nodes to be outside the problem domain. In this way, the integration points close to the domain boundary have the same number of MLS-nodes for interpolation as those far from the domain boundary. The accuracy is accordingly improved. The background grids are the Delaunay triangular meshes generated from the MLS-nodes and then cut by the problem domain. The numerical integration over a triangular cell containing no crack tip is carried out by the Hammer

162

H. Zheng et al. / Comput. Methods Appl. Mech. Engrg. 295 (2015) 150–171

Fig. 6. A plate with a symmetric crack under a uniform tension.

quadrature with six integration points. The numerical integration over a triangular cell containing a crack tip is calculated by the procedure proposed in [41]. For all the examples, the parameters ξ , η and γ , in the algorithm PGA, are taken to the values proposed by He in [40], namely, ξ = 0.9, η = 0.4, and γ = 1.9. The tolerance for the convergence in the algorithm PCA is ε = 0.1%. 7.1. A symmetric central straight crack Shown in Fig. 6 is a central straight crack with the boundary conditions. The example is made absolutely symmetric so as to check the ability of the proposed procedure to treat such extreme cases. An excellent numerical procedure should be able to reproduce the behaviors of symmetric problems as far as possible, no matter whether the entities for discretization, such as the node configuration and the background cell, are symmetrically deployed. The parameters include: L = 2.6 m, W = 1.6 m, the crack length a = 0.3 m, Young’s modulus E = 200 GPa, Poisson’s ratio ν = 0.3, fracture toughness K IC = 1300 MPa mm1/2 . The plane stress condition is assumed, and the maximum crack increment is 0.04 m. A uniformly-deployed MLS-node distribution is set with a node span of 0.02, 0.04, and 0.06 m, respectively, as shown in Fig. 7 for the node configuration with a node span of 0.04 m. Meanwhile, in order to demonstrate the capability of the proposed procedure to reproduce the behaviors of absolutely symmetric problems, the sample is rotated about its center by an angle of 15◦ , 30◦ and 45◦ , respectively, see Fig. 8; and let the maximum crack growth length l be 0.04 m. After rotation, the MLS-node deployment and the background grid become no longer symmetric about the problem domain. However, for all the cases, unless the tolerance ε K for the definition of equivalence of K eq and K IC is below 2%, both the crack tips grow symmetrically, suggesting that the proposed procedure is able to reproduce the behavior exhibited in the symmetric problems. If ε K is set below 1%, only one crack tip grows for the problem domain rotated by 30◦ . But whether the crack grows symmetrically or not, the curves of nominal stress vs. strain almost coincide in all the settings, as shown in Fig. 9. 7.2. A plate with two holes and edge cracks A rectangular plate with two off-center circular holes and two edge cracks is shown in Fig. 10, where the length unit is centimeter, with the material properties as E = 200 GPa, ν = 0.3 and K IC = 1300 MPa mm1/2 . The sample is fixed along its bottom in both horizontal and vertical directions. The top boundary is given an equal vertical displacement as deformation, with no horizontal displacement. So, the displacement-control is conducted and the problem is symmetric about its center.

H. Zheng et al. / Comput. Methods Appl. Mech. Engrg. 295 (2015) 150–171

163

Fig. 7. Uniformly deployed MLS-nodes.

Fig. 8. Rotation of sample and cracking paths.

An equal node span of 2 mm is set, as shown in Fig. 11, but the background grid is the Delaunay mesh automatically generated and not symmetric about the plate center. The maximum cracking increment l is taken as 2, 3 and 4 mm, respectively, to check if the proposed procedure is sensitive to the cracking increment. Let ε K = 10−3 , a very stern tolerance. Fig. 12 displays the cracking paths under three distinct cracking increments, all of them almost coincide, indicating the proposed procedure is insensitive to the cracking increment. As a comparison, the path given by Azadi and Khoei [42] is also drawn in the same figure. Fig. 13 illustrates the variations of K eq of the two crack tips during cracking. It clearly tells us at which steps which crack tips grow and which crack tips keep still—a crack tip is growing within a load step if its value of K IC /K eq equals unity. It also tells us at most steps only one crack tip is growing, which is believed to be consistent with brittle fracture, quite different from the cases of dynamic propagation of multiple cracks. The curve at the top of Fig. 14, called the full process curve, depicts the relationship between the vertical force and displacement on the top boundary under the displacement control. Such a curve has never been reported in the

164

H. Zheng et al. / Comput. Methods Appl. Mech. Engrg. 295 (2015) 150–171

Fig. 9. Full process curve of nominal stress vs. nominal strain (Example 7.1).

Fig. 10. Problem definition of Example 7.2 (length unit: cm).

Fig. 11. Node configuration for Example 7.2.

Fig. 12. Cracking paths under distinct cracking increments (da: cracking increment) given by the proposed method, with the result from Azadi [42].

H. Zheng et al. / Comput. Methods Appl. Mech. Engrg. 295 (2015) 150–171

165

Fig. 13. Variations of K eq ’s of all crack tips during cracking with a cracking increment of 4 mm.

Fig. 14. The full process curve (upper) depicts the relationship between the vertical load and displacement on the top boundary. The digits on both the figures represent the loading steps. A point “×” on the full process curve means that at this step the cracking path intersects with a problem boundary, such as step 39.

literature, and can be obtained only if the fracture toughness condition is obeyed. The digits 0, 2, . . . , on both the full process curve and the cracking path, representing the cracking steps, tell us not only the sequences of cracking, but also the correspondence between the points on the cracking path and the full process curve. Fig. 14 also tells us that

166

H. Zheng et al. / Comput. Methods Appl. Mech. Engrg. 295 (2015) 150–171

Fig. 15. The full process curve corresponding to symmetric cracking, with a relaxed tolerance of ε K = 3%.

the two crack tips grow simultaneously during steps 8–10 because both the left path and the right path are marked by the same step numbers; but beyond these steps, only one crack tip is growing. Interestingly enough, Fig. 14 appears to subvert such a consensus in the fracture mechanics literature: for brittle fracture, once cracking takes place it enters an unstable growth until complete failure of the structure, like the cracking in example 7.1. Nevertheless, Fig. 14 tells us that after cracking takes place, steps 2–4 and steps 6–10 are still stably cracking. If unloading at step 5 and then reloading, the sample is still able to resist deformation stably till step 10 is finished, at which the load level is very close to the collapse load of the plate before cracking. Since the problem is symmetric about its center, it is also possible for the two crack tips to grow symmetrically, although it is an unstable solution. As expected, if tolerance ε K , for defining the equivalence of K eq and K IC , is relaxed to, 3% say, the symmetric cracking path is reached, as shown in Fig. 15, together with the corresponding full process curve. Now we can see that symmetric cracking is able to experience a larger deformation than asymmetric cracking. To compare the two types of cracking, the full process curves are both drawn in Fig. 16, where the latter part of the curve of symmetric cracking is truncated. It can be seen that the curve of symmetric cracking has only two peak points, while that of asymmetric cracking has three peak points. Like the first example, see Fig. 9, whether cracking is symmetric or not, the collapse load is unique, which is believed to be significant in the strength evaluation of structures. For the convenience of comparison, the symmetric cracking paths given by Bouchard et al. [43] are displayed in Fig. 17, by copying and pasting. In reality, symmetric cracking is by no means accessible, because no ideally symmetric samples can be made out and no ideally loading can be realized. Example 7.3. A plate containing 10 random cracks A square plate of 2 m wide is considered containing 10 initial cracks, subjected to a uniform tension traction along the top and bottom boundary, as shown in Fig. 18, with E = 200 GPa, v = 0.3, and K IC = 1300 MPa mm1/2 . The plane stress condition is assumed.

H. Zheng et al. / Comput. Methods Appl. Mech. Engrg. 295 (2015) 150–171

Fig. 16. The full process curves of symmetric and asymmetric cracking.

Fig. 17. Cracking paths given by Bouchard et al. [43].

Fig. 18. Problem definition of Example 7.3.

167

168

H. Zheng et al. / Comput. Methods Appl. Mech. Engrg. 295 (2015) 150–171

Fig. 19. Variations of K eq ’s of all crack tips with loading.

An equal node span of 0.03 m is set, leading to an equal node distribution of 69 × 69 (= 4761). The maximum cracking increment is taken as 0.045 m. Shown in Fig. 19 are the variations of K eq of the 20 crack tips with loading. Again, it demonstrates the fact that in the most loading steps only one crack tip is growing at one moment, but at another moment some other crack tip may be growing. Although there are twenty crack tips present, at last only four of them grew during the entire deformation history, while others kept still. Fig. 20 illustrates the crack paths at different loading steps. Fig. 21 is the full process curve depicting the relationship between the nominal stress and strain. By programming with Matlab and running on a PC with Intel(R) Core(TM) i3-2120 CPU @3.30 GHz 3.30 GHz and a memory of 4.00 GB, the solution of this example takes about 90 min, a tolerable efficiency, we believe. This example was also analyzed by Azadi and Khoei [42], and Budynet al. [44]. Their final cracking paths are shown in Fig. 22, for the convenience of comparison. 8. Discussions and conclusions The static growth of multiple cracks can be reduced to a nonlinear complementarity problem (NCP) in this study, which has never been addressed in the literature. The NCP can be solved only by those algorithms that do not need the Jacobians, such as the projection-contraction algorithm (PCA) as designed in this study. In most of the cracking time, only one crack tip grows at one moment, although at different cracking periods, different crack tips might grow. Even in the presence of multiple cracks, there are only few cracks grown in the whole cracking period. For brittle fracture, there might be more than one critical point on the full curve of nominal stress vs. strain. Therefore, it is worthy to reconsider the conclusion that for brittle fracture once cracking takes place, it will be unstably growing until complete failure of the structure. Even if there exists more than one cracking path for a particular loading way, the collapsed load is unique. The complementarity problem built in this study is well-scaled, which is able to reproduce those behaviors of ideally symmetric problems. Compared with the FE-based NMM, the MLS-based NMM considerably simplifies cutting operations, worthy of further investigations. The proposed procedure is able to simulate growth of multiple cracks in a natural way, which allows the crack tips stop anywhere with no sensitivity to cracking increment or node-span.

H. Zheng et al. / Comput. Methods Appl. Mech. Engrg. 295 (2015) 150–171

Fig. 20. Cracking paths at distinct load steps given by the proposed method.

Fig. 21. The full process curve of nominal stress vs. strain.

169

170

H. Zheng et al. / Comput. Methods Appl. Mech. Engrg. 295 (2015) 150–171

Fig. 22. Final cracking paths given by Azadi and Khoei [42] and Budyn et al. [44].

Acknowledgments This study is supported by the National Basic Research Program of China (973 Program), under the Grant Nos. 2011CB013505 and 2014CB047100; and the National Natural Science Foundation of China, under the Grant No. 11172313. References [1] T. Belytschko, Y.Y. Lu, L. Gu, Element-free Galerkin methods, Internat. J. Numer. Methods Engrg. 37 (1994) 229–256. [2] I. Babuˇska, J.M. Melenk, The partition of unity method, Internat. J. Numer. Methods Engrg. 40 (1997) 727–758. [3] T. Belytschko, T. Black, Elastic crack growth in finite elements with minimal remeshing, Internat. J. Numer. Methods Engrg. 45 (1999) 601–620. [4] T. Strouboulis, I. Babuska, K. Copps, The design and analysis of the generalized finite element method, Comput. Methods Appl. Mech. Engrg. 181 (2000) 43–69. [5] G.H. Shi, Manifold method of material analysis, in: Transactions of the 9th Army Conference on Applied Mathematics and Computing, Report No. 92-1, US Army Research Office, Minneapolis, MN, 1991, pp. 57–76. [6] G.Q. Chen, Y. Ohnishi, T. Ito, Development of high-order manifold method, Internat. J. Numer. Methods Engrg. 43 (1998) 685–712. [7] Q.H. Jiang, C.B. Zhou, D.Q. Li, A three-dimensional numerical manifold method based on tetrahedral meshes, Comput. Struct. 87 (2009) 880–889. [8] X.M. An, L.X. Li, G.W. Ma, H.H. Zhang, Prediction of rank deficiency in partition of unity–based methods with plane triangular or quadrilateral meshes, Comput. Methods Appl. Mech. Engrg. 200 (2011) 665–674. [9] H. Ghasemzadeh, M.A. Ramezanpour, S. Bodaghpour, Dynamic high order numerical manifold method based on weighted residual method, Internat. J. Numer. Methods Engrg. 100 (2014) 596–619. [10] R. Tian, L.F. Wen, Improved XFEM–An extra-dof free, well-conditioning, and interpolating XFEM, Comput. Methods Appl. Mech. Engrg. 285 (2015) 639–658. [11] Y. Tal, Y.H. Hatzor, X.T. Feng, An improved numerical manifold method for simulation of sequential excavation in fractured rocks, Int. J. Rock Mech. Min. Sci. 65 (2014) 116–128. [12] Z. Wu, L.N.Y. Wong, Underground rockfall stability analysis using the numerical manifold method, Adv. Eng. Softw. 76 (2014) 69–85. [13] W.B. Zheng, X.Y. Zhuang, D.D. Tannant, D. Dwayne, Unified continuum/ discontinuum modeling framework for slope stability assessment, Eng. Geol. 179 (2014) 90–101. [14] Z. Wu, L. Fan, The numerical manifold method for elastic wave propagation in rock with time-dependent absorbing boundary conditions, Eng. Anal. Bound. Elem. 46 (2014) 41–50. [15] M. Kurumatani, K. Terada, Finite cover method with multi-cover layers for the analysis of evolving discontinuities in heterogeneous media, Internat. J. Numer. Methods Engrg. 79 (2009) 1–24. [16] D. Kourepinis, C. Pearce, N. Bi´cani´c, Higher-order discontinuous modeling of fracturing in concrete using the numerical manifold method, Int. J. Comput. Methods 7 (2010) 83–116. [17] K. Terada, T. Ishii, T. Kyoya, Y. Kishino, Finite cover method for progressive failure with cohesive zone fracture in heterogeneous solids and structures, Comput. Mech. 39 (2007) 191–210. [18] H.H. Zhang, G.W. Ma, Fracture modeling of isotropic functionally graded materials by the numerical manifold method, Eng. Anal. Bound. Elem. 38 (2014) 61–70. [19] L.Y.N. Wong, Z. Wu, Application of the numerical manifold method to model progressive failure in rock slopes, Eng. Fract. Mech. 119 (2014) 1–20. [20] X.M. An, Y.J. Ning, G.W. Ma, L. He, Modeling progressive failures in rock slopes with non-persistent joints using the numerical manifold method, Int. J. Numer. Anal. Methods Geomech. 38 (2014) 679–701.

H. Zheng et al. / Comput. Methods Appl. Mech. Engrg. 295 (2015) 150–171

171

[21] Y.J. Ning, X.M. An, G.W. Ma, Footwall slope stability analysis with the numerical manifold method, Int. J. Rock Mech. Min. Sci. 48 (2010) 2039–2071. [22] S.K. Yang, G.W. Ma, X.H. Ren, F. Ren, Cover refinement of numerical manifold method for crack propagation simulation, Eng. Anal. Bound. Elem. 43 (2014) 37–49. [23] H.H. Zhang, G.W. Ma, F. Ren, Implementation of the numerical manifold method for thermo-mechanical fracture of planar solids, Eng. Anal. Bound. Elem. 44 (2014) 45–54. [24] H. Zheng, Z.J. Liu, X.R. Ge, Numerical manifold space of Hermitian form and application to Kirchhoff’s thin plate problems, Internat. J. Numer. Methods Engrg. 95 (2013) 721–739. [25] Q.H. Jiang, S.S. Deng, C.B. Zhou, W.B. Lu, Modeling unconfined seepage flow using three-dimensional numerical manifold method, J. Hydrodyn. 22 (2010) 554–561. [26] Y. Wang, M.S. Hu, Q.L. Zhou, J. Rutqvist, Energy-work-based numerical manifold seepage analysis with an efficient scheme to locate the phreatic surface, Int. J. Numer. Anal. Methods Geomech. 38 (2014) 1633–1650. [27] H. Zheng, F. Liu, C.G. Li, Primal mixed solution to unconfined seepage flow in porous media with numerical manifold method, Appl. Math. Model. 39 (2015) 794–808. [28] D.S. Ling, Q.D. Yang, B.N. Cox, An augmented finite element method for modeling arbitrary discontinuities in composite materials, Int. J. Fract. 156 (2009) 53–73. [29] R. de Borst, Numerical aspects of cohesive zone models, Eng. Fract. Mech. 70 (2003) 1743–1757. [30] W. Liu, Q.D. Yang, S. Mohammadizadeh, X.Y. Su, An efficient augmented finite element method for arbitrary cracking and crack interaction in solids, Internat. J. Numer. Methods Engrg. 99 (2014) 438–468. [31] F. Facchinei, J.S. Pang, Finite-Dimensional Variational Inequalities and Complementarity Problems, Springer, New York, 2003. [32] H. Zheng, W. Jiang, Discontinuous deformation analysis based on complementary theory, Sci. China Ser. E 29 (2009) 2547–2554. [33] H. Zheng, X.K. Li, Mixed linear complementarity formulation of discontinuous deformation analysis, Int. J. Rock Mech. Min. Sci 75 (2015) 23–32. [34] H. Zheng, F. Liu, C.G. Li, The MLS-based numerical manifold method with applications to crack analysis, Int. J. Fract. 190 (2014) 47–166. [35] M. Spivak, Calculus on Manifolds, Benjamin, New York, 1965. [36] Y.C. Cai, X.Y. Zhuang, H.H. Zhu, A generalized and efficient method for finite cover generation in the numerical manifold method, Int. J. Comput. Methods 10 (2013) 1350028. [37] T. Belytschko, H. Chen, J. Xu, G. Zi, Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment, Internat. J. Numer. Methods Engrg. 58 (2003) 1873–1905. [38] X. Zhuang, C. Augarde, K.M. Mathisen, Fracture modeling using meshless methods and level sets in 3D: framework and modeling, Internat. J. Numer. Methods Engrg. 92 (2012) 969–998. [39] H. Zheng, J.L. Li, A practical solution for KKT systems, Numer. Algorithms 46 (2007) 105–119. [40] B.S. He, A class of projection and contraction methods for monotone variational inequalities, Appl. Math. Optim. 35 (1997) 69–76. [41] H. Zheng, D.D. Xu, New strategies for some issues of numerical manifold method in simulation of crack propagation, Internat. J. Numer. Methods Engrg. 97 (2014) 986–1010. [42] H. Azadi, A.R. Khoei, Numerical simulation of multiple crack growth in brittle materials with adaptive remeshing, Internat. J. Numer. Methods Engrg. 85 (2011) 1017–1048. [43] P.O. Bouchard, F. Bay, Y. Chastel, Numerical modelling of crack propagation: automatic remeshing and comparison of different criteria, Comput. Methods Appl. Mech. Engrg. 192 (2003) 3887–3908. [44] E. Budyn, G. Zi, N. Mo¨es, T. Belytschko, Numerical modelling of crack propagation: automatic remeshing and comparison of different criteria, Comput. Methods Appl. Mech. Engrg. 192 (2003) 3887–3908.