- Email: [email protected]
A concurrent multiscale method for simulation of crack propagation
Acta Mechanica Solida Sinica, Vol. 28, No. 3, June, 2015 Published by AMSS Press, Wuhan, China
ISSN 0894-9166
A CONCURRENT MULTISCALE METHOD FOR SIMULATION OF CRACK PROPAGATION⋆⋆ Jingkai Wu1,2
Hongwu Zhang1⋆
Yonggang Zheng1
(1 State Key Laboratory of Structural Analysis for Industrial Equipment, Department of Engineering Mechanics, Faculty of Vehicle Engineering and Mechanics, Dalian University of Technology, Dalian 116024, China) (2 Nanjing Research Institute of Electronics Technology, Nanjing 210039, China) Received 19 August 2013, revision received 17 April 2014
ABSTRACT A concurrent multiscale method is developed for simulating quasi-static crack propagation in which the failure processes occur in only a small portion of the structure. For this purpose, a multiscale model is adopted and both scales are discretized with finite-element meshes. The extended finite element method is employed to take into account the propagation of discontinuities on the fine-scale subregions. At the same time, for the other subregions, the coarse-scale mesh is employed and is resolved by using the extended multiscale finite element method. Several representative numerical examples are given to verify the validity of the method.
KEY WORDS concurrent multiscale method, crack propagation, extended finite element method, level set method
I. INTRODUCTION Localization phenomena, such as cracks and shear bands, have been recognized as the main causes of structural failures. In order to predict the life span of materials and analyze the causes of damage, one needs to investigate the behavior of localization phenomena in the structure. For the failure processes, the interesting physics are typically concentrated on a small portion of the structure. Usually, the initial element size of the structure is sufficient to evaluate the global response in the finite element analysis if it is not too coarse for simulation of the localization phenomena induced by the evolving cracks or shear bands. On the other hand, engineers prefer not to change the initial element mesh of the structure during the computational process. In this context, multiscale methods have made a great promise of modeling and simulating these kinds of problems. In recent decades, a large number of multiscale computational methods have been developed. In general, two different kinds of multiscale approaches can be distinguished. The first one is based on the homogenization schemes[1–3] , in which macroscopic equivalent parameters are obtained through the structural responses of representative volume elements (RVEs) that are modeled on the fine scale. This kind of approaches requires that the small-scale length should be much smaller than the coarse-scale length, i.e., the scales must be separated. This restrictive condition can be violated in the vicinity of Corresponding author. E-mail: [email protected] Project supported by the National Natural Science Foundation of China (Nos. 11232003, 11072051, 11272003 and 91315302), the 111 Project (No. B08014), the National Basic Research Program of China (Nos. 2010CB832704 and 2011CB013401), Program for New Century Excellent Talents in University (No. NCET-13-0088) and Ph.D.Programs Foundation of Ministry of Education of China (No. 20130041110050). ⋆
⋆⋆
· 236 ·
ACTA MECHANICA SOLIDA SINICA
2015
critical regions, however, since the failure processes occur on different scales where the shear bands or evolving cracks have to be explicitly modeled on the fine scale. To overcome this difficulty, several global-local methods (also called concurrent multiscale methods) based on hierarchical decomposition technique are developed. The entire region is divided into several subdomains. Critical subdomains, where the shear bands form or where microcracks propagate, are modeled on the fine scale by detailed microstructural computation. A widely used method is the variational multiscale method (VMM), which was initially introduced by Hughes et al.[4, 5] . The VMM has been used to simulate strain localization of both single and multiple dimensions[6–8] and crack propagation problems[9, 10] . Based on the LATIN method, a multiscale XFEM for crack propagation was proposed by Guidault et al.[11] . The method enables one to use a refined mesh only in the vicinity of discontinuities. In Ref.[12], a multiscale projection method was developed to simulate the effects of both macrocracks and microcracks. Belytschko et al.[13] developed a coarse-graining method to handle the multiscale crack propagation with the multiscale aggregating discontinuity method (MAD). There are also some other multiscale methods developed for discontinuous problems, to which the readers are referred[14–18] for more information. Recently, a new multiscale method, which is called the extended multiscale finite element method (EMsFEM), was developed by Zhang et al.[19, 20] for solving solid mechanical problems of multiple scale features. The main idea of the method is to numerically construct the multiscale base functions that can efficiently capture the small-scale properties within each coarse element by locally solving the Dirichlet boundary value problems. In order to consider the coupling effect on different directions in vector field problems, additional coupling terms are added to the multiscale base functions for interpolation of the displacement field. The important small-scale information is then gathered into a coarser scale through the coupling of the equivalent stiffness matrices. The method can be conveniently implemented and adopted to general problems without scale separation and periodicity assumptions. Moreover, the actual microscopic stress and strain within the coarse-scale elements can be obtained simultaneously along with the macroscopic displacement field by means of downscaling computation. The EMsFEM has been successfully used for elasto-plastic[19] and dynamic analyses[20] of heterogeneous materials without the restriction of scale separation and periodic assumptions. It should be noted that the idea of using numerical base functions had earlier been applied to solving elliptic problems arising from heat conductions in composite materials and flows in porous media (See Refs.[21, 22]). Most of these methods are devoted to solving scalar field problems. The simulation of crack propagation with the finite element method (FEM) has been a topic of growing interest in recent years. Various numerical methods have been developed to overcome the meshing difficulty faced when an evolving crack is modeled. Notable among them is the XFEM, which was first introduced by Belytschko et al.[23] . Then, a Heaviside step function in conjunction with twodimensional (2D) asymptotic crack-tip displacement fields are used to account for the crack in Ref.[24]. This allows the crack geometry to be represented independently of the finite element mesh, making remeshing unnecessary for model crack growth. Stolarska et al.[25] presented an algorithm that couples the LSM with the XFEM to solve the crack problem. The LSM was used to represent the crack location, including the location of crack tips, and the XFEM was used to compute the stress and displacement fields necessary for determining the rate of crack growth. The combination of the LSM and the XFEM seems to be quite elegant since the values of enrichment functions can be obtained by the level set functions. The LSM was further used for representing the location of holes and material interfaces by Sukumar et al.[26] . The XFEM has been widely used for modeling strong (displacement) and weak (strain) discontinuities in both 2D and 3D problems. Several surveys on the XFEM can be seen in Refs.[27–29]. Moreover, implementation of the XFEM for crack growth with the finite element codes can be found in Refs.[30–32]. In the present study, a new concurrent multiscale technique is proposed for crack propagation; whichs enables one to use a refined mesh only in the vicinity of discontinuities. For this purpose, the EMsFEM serves as a multiscale framework. During the loading process, the initial macroscale model is transformed into a multiscale model, in which the localized subregions (the vicinity of the discontinuities) are automatically detected and resolved by the XFEM on a fine-scale mesh. Meanwhile, the other subregions that surround the refined zone are resolved by the EMsFEM on a coarse-scale mesh. The displacements on the two different scales are then bridged at the scale-interface by virtue of the multiscale base functions. Level sets are employed here to describe the geometry of the discontinuities
Vol. 28, No. 3 Jingkai Wu et al.: A Concurrent Multiscale Method for Simulation of Crack Propagation · 237 ·
as well as to identify the fine-scale subregions. Note that the crack propagation is usually concentrated in a very small portion of the structure. Hence the dramatic reduction in computational cost.
II. EMsFEM AND ITS IMPLEMENTATION FOR ELASTIC PROBLEMS In this section, we briefly introduce the basic implementation process of the EMsFEM for elastic problems. The readers are referred to Refs.[19,20] for more details. For the EMsFEM, the main task is to numerically construct the multiscale base functions that can capture the small-scale properties within each coarse element by locally solving the Dirichlet boundary value problems. Then the smallscale information is gathered into large-scale through the equivalent stiffness matrices. Consider a 2D heterogeneous material structure occupying a region Ω and having a boundary Γ , as depicted in Fig.1. The equilibrium equation and boundary conditions can be expressed as
Fig. 1. Schematic description of the EMsFEM.
div (D : E (u)) = f σ·n=T u=u
in Ω on Γ σ on Γ u
(1)
where D is the fourth-order stiffness tensor representing material properties, E (u) is the strain tensor given as 1 T E (u) = ∇u + (∇u) (2) 2 and u is the displacement vector; σ is the stress vector which has three independent components, T T σ = σx σy τxy ; f is the body force vector, f = f x f y ; Γσ and Γu are the force and displacement boundary conditions, respectively; n is the unit normal vector. For the elastic problem, there are two computational steps that need to be performed in the EMsFEM. Firstly, the multiscale base functions of coarse elements (unit cells) are constructed numerically; then, upscaling and downscaling computations are implemented by virtue of the multiscale base functions. 2.1. The Construction of Multiscale Base Functions Take one of the coarse elements shown in Fig.1 as an example, which occupies a subregion Ω K , K Ω ⊂ Ω. The multiscale base functions of the element are constructed by solving the equilibrium equation in the subregion Ω K with specified boundary conditions. In refs.[19,20], several different kinds of boundary conditions have been proposed. From Eqs.(1) and (2), the general expression for solving the base functions of a 2D problem can be given as follows: LN ¯i = 0 in Ω k (i = 1, 2, 3, 4) (3) N ¯i (x, y) affined on ∂Ω k
· 238 ·
ACTA MECHANICA SOLIDA SINICA
2015
T where L is the elasticity operator and satisfies Lu = div D : 12 ∇u + (∇u) . Once the base functions are constructed, the microscopic displacement fields within a coarse element can be expressed as 4 4 X X u= N¯ixx u¯i + N¯ixy v¯i (4) i=1
v=
4 X
i=1
N¯iyy v¯i +
i=1
4 X
N¯iyx u¯i
(5)
i=1
where u and v are the microscopic displacement components in the x and y directions, respectively; N¯iyx (or N¯ixy ) is an additional coupling term indicating the displacement field in the y-direction (or x-direction) within the coarse element induced by unit displacement of macroscopic node ¯i in the x-direction (or y-direction). Equations (4) and (5) can be expressed in a unified form: u = Nu ¯E
(6)
where N is the matrix of multiscale base functions, u is the displacement vector of nodes on the fine-scale meshes within the unit cell, and u ¯E is the displacement vector of nodes on the coarse-scale level. They can be written as T u = u1 v1 u2 v2 · · · · · · un vn T T (7) N = R1 RT · · · RT 2 n T u ¯E = u¯1 v¯1 u¯2 v¯2 u¯3 v¯3 u¯4 v¯4
where
N¯1xx (i) N¯1xy (i) Ri = N¯1yx (i) N¯1yy (i) (i = 1, 2, ......, n)
N¯2xx (i) N¯2yx (i)
N¯2xy (i) N¯2yy (i)
N¯3xx (i) N¯3yx (i)
N¯3xy (i) N¯3yy (i)
N¯4xx (i) N¯4yx (i)
N¯4xy (i) N¯4yy (i)
(8)
and n is the total number of nodes on the fine-scale mesh within the unit cell. It should be mentioned that, for materials with periodic microstructures, the multiscale base functions need to be constructed only once on the unit cell and are used for all the coarse elements. 2.2. Upscaling and Downscaling Computations After the multiscale base functions are constructed, we can implement the upscaling and downscaling computations conveniently: the upscaling computation is adopted to obtain the macroscopic equivalent stiffness matrix of the coarse element and the downscaling computation to obtain the microscopic stress and strain information within the coarse element. The equivalent stiffness matrix of a coarse element can be obtained from the point of view of equivalence of strain energy, i.e., KE =
p X
K ′e ,
K ′e = GT e K e Ge
(9)
e=1
where p is the total number of fine-scale elements within the coarse element, and K e is the traditional stiffness matrix of a fine-scale element e (see Fig.1), which can be expressed as Z Ke = BT (10) e D e B e tdΩe Ωe
in which B e and De are the strain-displacement matrix and material property matrix of the element e, respectively; t is the element’s thickness; Ge is a transition matrix that denotes the mapping relation between the displacement vectors of the microscopic scale nodes and the macroscopic scale nodes: T ue = Ge u ¯E , Ge = Re1 Re2 Re3 Re4 (11)
Vol. 28, No. 3 Jingkai Wu et al.: A Concurrent Multiscale Method for Simulation of Crack Propagation · 239 ·
where ue is the displacement vector of the element e. On the other hand, the strain and stress of a fine-scale element (take the eth fine-scale element in the unit cell as an example) can be obtained by εe = B e u e = T e u ¯E ,
T e = B e Ge
(12)
S e = D e B e Ge
(13)
and σ e = D e εe = Se u ¯E ,
III. MODELING OF DISCONTINUITIES ON THE SMALL SCALE In this paper, the discontinuities induced by the material interfaces and crack are modeled on the small scale (fine scale). Consider a body Ω f ⊂ R2 that contains nm circular inclusions and an internal traction-free crack (see Fig.2). We denote the boundary of Ω f by Γ f , which is not only decomposed into displacements Γuf and traction Γσf boundaries, but also into material interfaces Γm = ∪i Γm,i (i = 1, · · · nm ) and crack surface Γc , i.e., Γ f = Γuf ∪ Γσf ∪ Γm ∪ Γc . The domains of inclusions are denoted as Ωin,i , in which Ωin,i ⊂ Ω f . Similar to Eq.(1), the equilibrium equations for the corresponding problem can be expressed as in Ω f div (D : E f (uf )) = f σf · n = T f on Γσf (14) uf = u ¯f on Γuf σf · n = 0 on Γc+ σf · n = 0 on Γc− where uf , σ f , E f are the displacement fields, stress vector, and strain vector on the fine scale, respectively. Both weak discontinuities and strong discontinuities are considered here. The solution to the equation set (14) is worked out using the XFEM. 3.1. BASIC CONCEPTION OF THE XFEM The traditional FEM will encounter meshing difficulty when handling moving discontinuity problems, since the mesh needs to conform with the geometry of the crack. The XFEM provides an efficient technique to overcome this problem. In this section, the XFEM for the mechanical description of the discontinuities is briefly discussed. In the XFEM, the discontinuous enrichment functions are added to the finite element displacement approximation to represent the discontinuities within an element by virtue of the partition of the unity method[33] . In Fig.3, a crack is located on a uniform mesh. The displacement approximations in a classical XFEM[24] can be written as follows: uf (x) =
X
I∈S
NI (x) u ˆI +
X
NJ (x) H(x)bJ +
J∈SH
Fig. 2. Microstructure of a body with inclusions and crack.
X
K∈SC
NK (x)
4 X
Fα (x) cαK
(15)
α=1
Fig. 3. Enriched nodes for 2D crack problems on a uniform mesh. The circled nodes (set of nodes SH ) are enriched by step function and the squared nodes (set of nodes SC ) are enriched by crack tip functions.
· 240 ·
ACTA MECHANICA SOLIDA SINICA
2015
where S is the set of all nodes of the mesh; SH ⊂ S is the set of nodes whose support is bisected by a crack; and SC ⊂ S is the set of nodes whose support contains the crack tips. u ˆI is the classical degrees of freedom at node I; bJ and cαK are the enriched degrees of freedom at the corresponding enriched nodes; H(x) is the generalized Heaviside function whose value is set to +1 above the crack and −1 below the crack; Fα (x) is a function that models the near tip asymptotic field: √ √ √ √ θ θ θ θ {Fα (x)} = (16) r sin , r cos , r sin sin θ, r cos sin θ 2 2 2 2 where r and θ are polar coordinates of the local axes at the tip of crack (xtip , ytip ). The approximation in Eq.(15) will give rise to the emergence of blending elements between the enriched and unenriched elements. In these blending elements, only part of their nodes are enriched, so the approximation does not satisfy the partition of unity. This can lead to a decrease in the rate of convergence and errors in the blending elements. Several modifications have been proposed to circumvent the blending problem (see Refs.[34–36]). For example, the Heaviside enrichment function in the standard XFEM can be replaced with the shifted enrichment function, which can be written as ψJ (x) = NJ (x) (H(x) − H(xJ ))
(17)
Thus, the interpolation property is satisfied since the shifted enrichment function for the discontinuity vanishes at the edges of elements. Also, the enrichment function F can be expressed as X X F (x) = |φI | NI (x) − φI NI (x) (18) I
I
where φI is the level set value at the node I for the material interface level set function (see Ref.[36]). For inclusion problems, the following form can be taken to represent a material interface in the displacement approximation[37]: X X uf (x) = NI (x) u ˆI + aL NL (x) F (x) (19) I∈S
L∈SI
where SI ⊂ S is the set of nodes whose support is split by the interface, and aL are the enriched degrees of freedom at such nodes. 3.2. Geometrical Description of Discontinuities The LSM proposed by Osher and Sethian[38] is used for geometrical description of the discontinuities. Consider a region Ω which is split into two subdomains ΩA and ΩB by the interface Γ . The main idea of the LSM is to represent the interface Γ by means of a scalar function φ(x), which has the following properties: φ(x(t), t) < 0 for x ∈ ΩA φ(x(t), t) > 0 for x ∈ ΩB (20) φ(x(t), t) = 0 for x ∈ Γ Therefore, a moving interface at a given time t can be formulated as Γ (t) = x ∈ R2 |φ(x, t) = 0
(21)
The evolution of the interface can then be governed by solving the Hamilton-Jacobi equation ∂φ(x) + vn (x) |∇φ(x)| = 0 ∂t
(22)
where vn (x) is the velocity of the level set field. In the case of circular inclusions shown in Fig.2, the geometry of the fixed material interfaces Γm = ∪j Γm,j can be described by the following level set function: φm (x) = min {kx − xm,i k − rm,i } | {z } φm,i
for
i = 1, · · · ,nm
(23)
Vol. 28, No. 3 Jingkai Wu et al.: A Concurrent Multiscale Method for Simulation of Crack Propagation · 241 ·
where xm,i is the location of the center of the ith inclusion whose radius is rm,i . In order to track the growth of a 2-D crack , two level set functions φ and ψ are needed, one for the crack’s surface (φ) and the other for the crack tip (ψ). Thus, the crack’s position Γ (t) can be defined by means of the updated values of φ and ψ: Γ (t) := x ∈ R2 |φ(x, t) = 0 and ψ(x, t) ≤ 0 (24)
In general, the scalar values of the level set function are computed only at the finite element nodes, and the level set function can be interpolated by the finite element shape functions X φ(x) = Ni (x)φi (25) i∈I
Thus, the true locations of the discontinuities are approximated by the iso-zero of the level set function φ. Generally, the mesh should be sufficiently fine in the vicinity of discontinuities for their precise location and reduction in numerical errors (see Ref.[37]). As an example, we consider a simple structure with an inclusion shown in Fig.4. A coarse and a fine mesh each are used to model the material interface. Figure 5 shows the level set function for the coarse mesh. We can see that the geometrical representation is poor in this case. At the same time, by use of the fine mesh (see Fig.6), a much better geometrical description for the interface is obtained by the level set function. An adaptive mesh refinement strategy will be introduced in the next section to obtain a finer mesh only in the vicinity of discontinuities.
Fig. 4 Geometry of a body with an inclusion.
Fig. 5. Level set functions in 3D and 2D axes on a coarse-scale mesh.
Fig. 6. Level set functions in 3D and 2D axes on a fine-scale mesh.
· 242 ·
ACTA MECHANICA SOLIDA SINICA
2015
Fig. 7. Illustration of the CMM for crack propagation.
IV. CONCURRENT MULTISCALE METHOD FOR MODELING DISCONTINUITIES In engineering practice, the initial element size of the structure is always sufficient to evaluate the global structural response, but it is not fine enough to represent the discontinuities induced by the fixed material interfaces or the evolving cracks. As pointed out in §3.2, a mesh refinement treatment is required in the vicinity of the discontinuities in order to precisely represent their geometry while accurately resolving the displacement field. On the other hand, engineers prefer not to change the initial finite element mesh of the structure during the evolution of discontinuities since the operation is extremely difficult. In this section, a concurrent multiscale method (CMM) is developed for failure analysis including both weak and strong discontinuities. During the loading process, the entire domain Ω of structure is automatically divided into two different regions, i.e., the coarse-scale discrete region Ω c and the fine-scale discrete region Ω f (see Fig.7). The two regions are solved independently whereby the discontinuities are taken into account on the fine-scale mesh. For this purpose, the EMsFEM is here employed to capture the macroscopic structural response to the coarse-scale region with moderate strains, whereas the XFEM is used to model the material interfaces and crack propagation on the fine-scale region. In this context, the discontinuities are solved by detailed microstructural simulation. The geometry of the material interfaces and the evolving cracks are described on the fine scale by means of the LSM. Then the displacements between the two scales can be bridged at the scale-interface Γ cf using the multiscale base functions. The submodels with two different scales are then combined into one numerical model, the so-called mixed structure model. Note that the discontinuities are usually located in a very small portion of the structure, so that the expensive calculation of detailed fine-scale simulation is limited to a minimum. 4.1. Discretization of Domain An adaptive procedure based on the distance of the element to the discontinuities is used to switch between the two levels of discretization. Take Fig.7 as an example. An evolving crack is located on the initial coarse mesh. Prior to each growth step, the distance of each macroscopic node to the crack tip d1 (x, t) is computed. If the value of d1 (x, t) at any one of the nodes in a coarse element is less than the prescribed critical value R1 , then the corresponding coarse element becomes a fine-scale subregion, and vice versa. The fine-scale subregion can be retransformed into a coarse element when all the values of d1 (x, t) at the four macroscopic nodes exceed R1 . The value of R1 can be set to be approximately 1.5 times the coarse element size. Note that in our method the crack is modeled only on the fine scale; thus the fine-scale subregions have to be adopted everywhere along the crack path. Moreover, in order to ensure that the enriched nodes are not present on the scale-interface Γ cf , a coarse element close to the crack should be replaced by a fine-scale subregion if the distance to the crack surface d2 (x, t) at any one of the macroscopic nodes in the coarse element is less than the prescribed critical value l1 . Here, the value of l1 is approximately 2 times the fine element size. Fortunately, finding that the values of
Vol. 28, No. 3 Jingkai Wu et al.: A Concurrent Multiscale Method for Simulation of Crack Propagation · 243 ·
d1 (x, t) and d2 (x, t) at each coarse-scale node can be obtained by virtue of the level set functions, we have p d1 (x, t) = ψ 2 (x, t) + φ2 (x, t), d2 (x, t) = |φ(x, t)| (26) For the fixed material interface, the same technique can be used to identify the fine-scale subregion. The distances of nodes to the material interface d3 can be obtained by d3 (x) = |φm (x)|
(27)
Fig. 8. Illustration of the CMM for a material interface.
If d3 (x) at any one of the nodes in a coarse element is less than the prescribed critical value l2 , or the element is bisected by the material interface, then the corresponding coarse element becomes a fine-scale subregion. In addition, l2 is equal to l1 in this paper. By taking the structure in Fig.4 as an example, the entire domain Ω is divided into Ω c and Ω f , as shown in Fig.8. In this case, the level set function for the inclusion is plotted in Fig.9. We can see that a good geometrical representation for the interface can be obtained. It should be mentioned that the distance based criterion pro- Fig. 9 Level set functions in 2D axes on a posed above could be replaced by a more sophisticated error es- multiscale mesh. timator in our multiscale scheme; the corresponding work is still going on. 4.2. Bridging between Two Different Scales As mentioned above, two levels of discretization are used in the CMM for the implementation of the EMsFEM and XFEM, respectively. Thus, tying relations should be employed to guarantee compatibility of the displacements at the scale-interface Γ cf between the two regions. From the viewpoint of the EMsFEM we know that the relationship between the microscopic and macroscopic displacement fields within a coarse element can be established by means of the multiscale base functions. Take advantage of this feature and the master-slave constraints can be adopted here to bridge the degrees of freedom between the two different scales. As shown in Fig.10, the microscopic scale nodes at the scale-interface are defined as slave nodes (for example, node i), while the nodes on the adjacent coarse element are defined as master nodes (such as nodes I, J, K, L). The relation of displacements between the salve and master nodes can be expressed as uf (i) = T i u ¯E(I,J,K,L) (28) where uf(i) is the displacement vector of the slave node i, u ¯E(I,J,K,L) is the displacement vector of the corresponding master nodes, T i is the control matrix of the node i and can be deduced from Eq.(11).
· 244 ·
ACTA MECHANICA SOLIDA SINICA
2015
Fig. 10. Master-slave constraints.
4.3. Solution Procedures The direction of crack propagation depends on the chosen physical criteria. In this paper, the maximum circumferential stress[39] is used to evolve a crack. The angle of crack growth can be expressed as s 2 1 K I − sign(KII ) 8 + KI (29) θc = 2 arctan 4 KII KII
where KI and KII are the mixed-mode stress intensity factors and can be calculated by virtue of the interaction integrals[24] . The constant crack growth increment is used to determine the magnitude of the incremental crack growth △a. The solution procedures for modeling quasi-static crack growth in the CMM can be summarized as follows: 1. Preprocessing: Construct the base functions N of the unit cell with the special boundary conditions, and derive the equivalent stiffness matrix K E of the coarse element. 2. According to the current discontinuities’ location, divide the entire domain Ω into the coarse-scale discrete region Ω c and the fine-scale discrete region Ω f . 3. Obtain the stiffness matrix K Ω c of the region Ω c and the stiffness matrix K Ω f of the region Ω f by use of the EMsFEM and XFEM, respectively. Bridge the degrees of freedom between the two different scales at the scale-interface Γ cf . Then, obtain the global stiffness matrix K mix of the mixed structure. 4. Input the displacements and traction boundary conditions. Solution for the macroscopic and microscopic (including the traditional and enriched) degrees of freedom. 5. Calculate KI , KII , θc and △a. Get the new location of crack tip. 6. Repeat the steps 2-5.
V. NUMERICAL EXAMPLES In this section, several representative numerical examples are presented to examine the validity of the proposed method. The various parameters influencing the results of the CMM are studied. Note that plane strain and linear elastic materials are assumed, and all numerical examples are dimensionless. Example 1 In the first numerical example, we investigate the influence of the different parameters (including the resolution of the coarse-scale mesh and the critical value R1 for the mesh refinement of region around the crack tip) on the result of the CMM. As shown in Fig.11, a structure with an edge crack of length 0.5a is loaded by prescribing displacement boundary conditions of mode I along all edges. The material
Vol. 28, No. 3 Jingkai Wu et al.: A Concurrent Multiscale Method for Simulation of Crack Propagation · 245 ·
Fig. 11. Sketch of a structure with edge crack.
parameters are chosen as follows: Young’s modulus E = 1.0 × 106 and the Poisson’s ratio ν = 0.3. The analytical expressions of displacement field (see the work of Griffith) can be written as r r θ KI θ KI r r cos (κ − cos θ), uy (r, θ) = sin (κ − cos θ) (30) ux (r, θ) = 2µ 2π 2 2µ 2π 2 where (r, θ) is the local polar coordinates, µ is the shear modulus and κ is Kolossov’s constant. The stress intensity factor is set to KI = 1.0. The coarse-scale meshes are composed of 9 × 9, 15 × 15, 27 × 27 and 45 × 45 elements, respectively. Moreover, the critical value R1 is set to R1 = 1.5 × lc and R1 = 3.0 × lc , respectively, where lc is the size of coarse-scale mesh. To eliminate the influence of the fine-scale resolution, the size of fine-scale mesh is chosen to be a constant whose value is h = a/405. Figure 12 shows the convergence of the stress intensity factor KI of the crack tip with respect to the number of coarse-scale elements. We can see that the results obtained by both cases can converge toward the theoretical value (1.0) when the number of coarse-scale elements is increasing, especially for the first increasing stage. This means that the number of coarse-scale elements 15 × 15 is sufficient for this test and can yield satisfactory results. It also shows that the accuracy can be obviously improved when a larger value is assigned to the critical value R1 . This is because the displacement gradients vary abruptly in the vicinity of the crack tip, and the fine-scale region should be large enough to simulate this phenomenon. However, more computational degrees of freedom are needed when the value R1 is larger. It should be mentioned that the results shown in Fig.12 are good even for the coarsest mesh resolution (deviation in only the third
Fig. 12. Convergence of stress intensity factor KI .
Fig. 13. The geometry of an edge crack problem.
· 246 ·
ACTA MECHANICA SOLIDA SINICA
2015
Fig. 14. Values of KI with respect to crack growth iteration.
Fig. 15. Geometry of a plate with an inclusion and a void.
digit); thus our multiscale method has a good accuracy. Example 2 Consider a finite plate with an edge crack of length a shown in Fig.13. The width and height of the plate are W and L, respectively. The following values are used, a = 2.5, W = 10, L = 19 and σ = 1. The theoretical Mode I stress intensity factor KI for this problem is given by[24] √ KI = F (λ) σ aπ
(31)
where F (λ) is a finite-geometry correction factor: F (λ) = 1.12 − 0.231λ + 10.55λ2 − 21.72λ3 + 30.39λ4 ,
λ=
a W
(32)
In each iteration, the crack is propagated in constant increment of △a = 0.25, and 15 iterative steps are considered. The number of the coarse-scale elements is kept constant at 10 × 19. In the fine-scale region, each of the coarse-scale elements is discretized by 9 × 9. The material parameters are the same as in Example 1. Figure 14 shows the values of the stress intensity factor KI obtained by the CMM and theoretical expression. It can be seen that the results obtained by the CMM fit fairly well with the exact solutions. Example 3 In this example, we study the problem of crack propagation in the presence of both a hard inclusion and a void. For the finite plate in the preceding example, here a circular inclusion with E = 5.0 × 106 , the Poisson’s ratio ν = 0.3 and radius R = 1.7 plus a circular void with radius R = 1.7 are located in the plate. The geometry of the plate is shown in Fig.15, where L1 = L2 = 5. The crack is propagated for 20 steps. For comparison, a full fine-scale simulation with 90 × 171 elements is accomplished, whose results can be seen as the reference solutions and are denoted as XFEM-F. Note that the size of elements used in the XFEM-F is the same as that of fine-scale elements in the CMM. The crack paths obtained by the two methods are compared in Fig.16, which are in excellent agreement. The crack propagates, as to be expected, close to the void by the influence of the inclusion and void. The mesh is refined only on the subregions where the discontinuities occur. Figure 17 shows the values of stress intensity factors KI and KII , respectively. Moreover, the distributions of the Von-Mises stress calculated by the two methods at three selected steps are illustrated in Fig.18. We can see that the results obtained by the CMM fit well with the reference solutions. For the computing times, the CMM requires about one-fifth the CPU time of the XFEM-F.
Vol. 28, No. 3 Jingkai Wu et al.: A Concurrent Multiscale Method for Simulation of Crack Propagation · 247 ·
Fig. 16. Crack paths at different iterative steps.
Fig. 17. Values of KI and KII with respect to crack growth iterations.
Fig. 18. Distributions of Von-Mises stress at different iterative steps.
Example 4 Let us consider the crack propagation of a three-point bending model, which is loaded by a single force F = 1 imposed at the center of the upper edge. A crack of initial length 0.6 is located on the lower edge. The material parameters are the same as in Example 1. The geometry and boundary conditions of the model are shown in Fig.19. In this study, a total of 20 iterative steps are considered for crack propagation, and the propagation increment is set to △a = 0.3. To test the mesh sensitivity of the multiscale method, two kinds of coarse-scale FE discretization of different element sizes are considered, in which the numbers of the coarse-scale elements are 10×4 and 20×8, respectively. The results obtained by the CMM for these two cases are denoted as CMM1 and CMM2, respectively. In the fine-scale region, square elements with length 1/11 are used. As is known that the stress concentration phenomenon will
· 248 ·
ACTA MECHANICA SOLIDA SINICA
2015
Fig. 19. Geometry and boundary conditions of the model.
emerge at the place where the forces and constraints are. Thus the meshes in the CMM are refined a priori in these subregions to reduce the errors. For the sake of comparison, a full fine-scale simulation with 220 × 88 elements is implemented by the XFEM. Figure 20 shows the crack paths for the XFEM-F, CMM1, and CMM2 at several selected steps. In addition, the deformed meshes for the three methods at the last step are plotted in Fig.21. The results obtained by both the CMM1 and CMM2 are in good agreement with those of the XFEM-F. The crack propagates in a curved path towards the centre of the upper edge. The values of stress intensity factors KI and KII as a function of the crack iterative steps are shown in Fig.22. Moreover, the distributions of Von-Mises stress at several selected steps for the three methods are plotted in Fig.23. The results further illustrate that the CMM proposed have a good accuracy. The method shows insensitivity to the coarse-scale discretization. The mesh refinement in the CMM needs only to be implemented in the vicinity of the crack, thus dramatically reducing the computer memory requirements and CPU time. For this example, the CMM2 requires about one-tenth the CPU time of the XFEM-F.
Fig. 20. Crack paths at different iterative steps.
Vol. 28, No. 3 Jingkai Wu et al.: A Concurrent Multiscale Method for Simulation of Crack Propagation · 249 ·
Fig. 21. Deformed meshes at the last iterative step.
Fig. 22. Values of KI and KII with respect to crack growth iterations.
Fig. 23. Distributions of Von-Mises stress at different iterative steps.
VI. CONCLUSIONS An concurrent multiscale method is developed to model quasi-static crack propagation in elastic material. The EMsFEM serves as a multiscale framework and is extended to take into consideration material failure. For this purpose, the combination of the XFEM and LSM is successfully integrated into the multiscale framework to simulate the crack propagation on a fine scale. Level sets are employed to describe the geometry of discontinuities as well as to automatically identify the fine-scale subregions. Then, the degrees of freedom between the two different scales are bridged at the scale-interfaces by means of multiscale base functions. Since the failure processes are usually concentrated in a very small portion of the structure, the expensive calculation of the fine-scale simulation can be limited to a minimum. Our numerical examples have shown that the proposed method can be successfully used for simulating the evolving discontinuities and for dramatically reducing the computational cost. Thus, the method seems to have great potential for failure analysis of large-scale structures.
References [1] Babuˇska,I., Homogenization approach in engineering. In: Lions,J.L. and Glowinski,R. (Eds.), Computing Methods in Applied Science and Engineering, Lecture Note in Economics and Mathematical Systems. Berlin: Springer, 1976, 134: 137-153. [2] Guedes,J.M. and Kikuchi,N., Preprocessing and postprocessing for materials based on the homogenization method with adaptive finite element methods. Computer Methods in Applied Mechanics and Engineering,
· 250 ·
ACTA MECHANICA SOLIDA SINICA
2015
1990, 83: 143-198. [3] Fish,J., Shek,K., Pandheeradi,M. and Shephard,M., Computational plasticity for composite structures based on mathematical homogenization: theory and practice. Computer Methods in Applied Mechanics and Engineering, 1997, 148: 53-73. [4] Hughes,T.J.R., Multiscale phenomena: Green’s function, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized method. Computer Methods in Applied Mechanics and Engineering, 1995, 127: 387-401. [5] Hughes,T.J.R., Feijoo,G.R., Mazzei,L. and Quincy,J.B., The variational multiscale method-a paradigm for computational mechanics. Computer Methods in Applied Mechanics and Engineering, 1998, 166: 3-24. [6] Garikipati,K. and Hughes,T.J.R., A study of strain localization in a multiple scale framework–the onedimensional problem. Computer Methods in Applied Mechanics and Engineering, 2000, 159: 193-222. [7] Garikipati,K. and Hughes,T.J.R., A variational multiscale approach to strain localization-formulation for multidimensional problems. Computer Methods in Applied Mechanics and Engineering, 2000, 188: 39-60. [8] Yeon,J.H. and Youn,S.K., Meshfree analysis of softening elastoplastic solids using variational multiscale method. International Journal of Solids and Structures, 2005, 42: 4030-4057. [9] Mergheim,J., A variational multiscale method to model crack propagation at finite strains. International Journal for Numerical Methods in Engineering, 2009, 80: 269-289. [10] Hettich,T., Hund,A. and Ramm,E., Modeling of failure in composites by X-FEM and level sets within a multiscale framework. Computer Methods in Applied Mechanics and Engineering, 2008, 197: 414-424. [11] Guidault,P.A., Allix,O., Champaney,L. and Cornuault,C., A multiscale extended finite element method for crack propagation. Computer Methods in Applied Mechanics and Engineering, 2008, 197: 381-399. [12] Loehnert,S. and Belytschko,T., A multiscale projection method for macro/microcrack simulations. International Journal for Numerical Methods in Engineering, 2007, 71: 1466-1482. [13] Belytschko,T. and Song,J.H., Coarse-graining of multiscale crack propagation. International Journal for Numerical Methods in Engineering, 2010, 81: 537-563. [14] Pereira,J.P.A., Kim,D.J. and Duarte,C.A., A two-scale approach for the analysis of propagating threedimensional fractures. Computational Mechanics, 2012, 49: 99-121. [15] Pierres,E., Baietto,M.C. and Gravouil,A., A two-scale extended finite element method for modelling 3D crack growth with interfacial contact. Computer Methods in Applied Mechanics and Engineering, 2010, 199: 1165-1177. [16] Holl,M., Loehnert,S. and Wriggers,P., An adaptive multiscale method for crack propagation and crack coalescence. International Journal for Numerical Methods in Engineering, 2013, 93: 23-51. [17] Kim,D.J., Pereira,J. P. and Duarte,C.A., Analysis of three-dimensional fracture mechanics problems: a two-scale approach using coarse-generalized FEM meshes. International Journal for Numerical Methods in Engineering, 2010, 81: 335-365. [18] Lian,W.D., Legrain,G. and Cartraud,P., Image-based computational homogenization and localization: comparison between X-FEM/levelset and voxel-based approaches. Computational Mechanics, 2013, 51: 279-293. [19] Zhang,H.W., Wu,J.K. and Fu,Z.D., A new multiscale computational method for elasto-plastic analysis of heterogeneous materials. Computational Mechanics, 2012, 49: 149-169. [20] Zhang,H.W., Liu,H. and Wu,J.K., A uniform multiscale method for 2D static and dynamic analyses of heterogeneous materials. International Journal for Numerical Methods in Engineering, 2013, 93: 714-746. [21] Hou,T.Y. and Wu,X.H., A multiscale finite element method for elliptic problems in composite materials and porous media. Journal of Computational Physics, 1997, 134: 169-89. [22] Hou,T.Y., Multiscale modelling and computation of fluid flow. International Journal for Numerical Methods in Fluids, 2005, 47: 707-719. [23] Belytschko,T. and Black,T., Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering, 1999, 45: 601-620. [24] Moes,N., Dolbow,J. and Belytschko,T., A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 1999, 46: 131-150. [25] Stolarska,M., Chopp,D.L., Moes,N. and Belyschko,T., Modelling crack growth by level sets in the extended finite element method. International Journal for Numerical Methods in Engineering, 2001, 51: 943-960. [26] Sukumar,N., Chopp,D., Moes,N. and Belytschko,T., Modeling holes and inclusions by level sets in the extended finite-element method. Computer Methods in Applied Mechanics and Engineering, 2001, 190: 6183-6200. [27] Karihaloo,B.L. and Xiao,Q.Z., Modelling of stationary and growing cracks in FE framework without remeshing: a state-of-the-art review. Computers and Structures, 2003, 81: 119-129. [28] Abdelaziz,Y. and Hamouine,A., A survey of the extended finite element. Computers and structures, 2008, 86: 1141-1151.
Vol. 28, No. 3 Jingkai Wu et al.: A Concurrent Multiscale Method for Simulation of Crack Propagation · 251 · [29] Fries,T.P. and Belytschko,T., The extended/generalized finite element method: An overview of the method and its applications. International Journal for Numerical Methods in Engineering, 2010, 84: 253-304. [30] Pais,M., MATLAB Extended Finite Element (MXFEM) Code, Version 1.2, 2010. (Available from: www.matthewpais.com. [Accessed on April 3, 2012]). [31] Sukumar,N. and Prevost,J.H., Modeling quasi-static crack growth with the extended finite element method part i: computer implementation. International Journal of Solids and Structures, 2003, 40: 7513-7537. [32] Bordas,S., Nguyen,P.V., Dunant,C., Guidoum,A. and Nguyen-Dang,H., An extended finite element library. International Journal for Numerical Methods in Engineering, 2007, 71: 703-732. [33] Melenk,J.M. and Babuska,I., The partition of unity finite element method: basic theory and applications. Computer Methods in Applied Mechanics and Engineering, 1996, 139: 289-314. [34] Belytschko,T., Moes,N., Usui,S. and Parimi,C., Arbitrary discontinuities in finite elements. International Journal for Numerical Methods in Engineering, 2001, 50: 993-1013. [35] Legay,A., Wang,H.W. and Belytschko,T., Strong and weak arbitrary discontinuities in spectral finite elements. International Journal for Numerical Methods in Engineering, 2005, 64: 991-1008. [36] Fries,T.P., A corrected XFEM approximation without problems in blending elements. International Journal for Numerical Methods in Engineering, 2008, 75(5): 503-532. [37] Moes,N., Cloirec,M., Cartraud,P. and Remacle,J., A computational approach to handle complex microstructure geometries. Computer Methods in Applied Mechanics and Engineering, 2003, 192: 3163-3177. [38] Osher,S. and Sethian,J., Fronts propagating with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics, 1988, 79: 12-49. [39] Erdogan,F. and Sih,G.C., On the crack extension in plates under plane loading and transverse shear. Transactions of the ASME. Journal of Basic Engineering, 1963, 85: 519-525.
ISSN 0894-9166
A CONCURRENT MULTISCALE METHOD FOR SIMULATION OF CRACK PROPAGATION⋆⋆ Jingkai Wu1,2
Hongwu Zhang1⋆
Yonggang Zheng1
(1 State Key Laboratory of Structural Analysis for Industrial Equipment, Department of Engineering Mechanics, Faculty of Vehicle Engineering and Mechanics, Dalian University of Technology, Dalian 116024, China) (2 Nanjing Research Institute of Electronics Technology, Nanjing 210039, China) Received 19 August 2013, revision received 17 April 2014
ABSTRACT A concurrent multiscale method is developed for simulating quasi-static crack propagation in which the failure processes occur in only a small portion of the structure. For this purpose, a multiscale model is adopted and both scales are discretized with finite-element meshes. The extended finite element method is employed to take into account the propagation of discontinuities on the fine-scale subregions. At the same time, for the other subregions, the coarse-scale mesh is employed and is resolved by using the extended multiscale finite element method. Several representative numerical examples are given to verify the validity of the method.
KEY WORDS concurrent multiscale method, crack propagation, extended finite element method, level set method
I. INTRODUCTION Localization phenomena, such as cracks and shear bands, have been recognized as the main causes of structural failures. In order to predict the life span of materials and analyze the causes of damage, one needs to investigate the behavior of localization phenomena in the structure. For the failure processes, the interesting physics are typically concentrated on a small portion of the structure. Usually, the initial element size of the structure is sufficient to evaluate the global response in the finite element analysis if it is not too coarse for simulation of the localization phenomena induced by the evolving cracks or shear bands. On the other hand, engineers prefer not to change the initial element mesh of the structure during the computational process. In this context, multiscale methods have made a great promise of modeling and simulating these kinds of problems. In recent decades, a large number of multiscale computational methods have been developed. In general, two different kinds of multiscale approaches can be distinguished. The first one is based on the homogenization schemes[1–3] , in which macroscopic equivalent parameters are obtained through the structural responses of representative volume elements (RVEs) that are modeled on the fine scale. This kind of approaches requires that the small-scale length should be much smaller than the coarse-scale length, i.e., the scales must be separated. This restrictive condition can be violated in the vicinity of Corresponding author. E-mail: [email protected] Project supported by the National Natural Science Foundation of China (Nos. 11232003, 11072051, 11272003 and 91315302), the 111 Project (No. B08014), the National Basic Research Program of China (Nos. 2010CB832704 and 2011CB013401), Program for New Century Excellent Talents in University (No. NCET-13-0088) and Ph.D.Programs Foundation of Ministry of Education of China (No. 20130041110050). ⋆
⋆⋆
· 236 ·
ACTA MECHANICA SOLIDA SINICA
2015
critical regions, however, since the failure processes occur on different scales where the shear bands or evolving cracks have to be explicitly modeled on the fine scale. To overcome this difficulty, several global-local methods (also called concurrent multiscale methods) based on hierarchical decomposition technique are developed. The entire region is divided into several subdomains. Critical subdomains, where the shear bands form or where microcracks propagate, are modeled on the fine scale by detailed microstructural computation. A widely used method is the variational multiscale method (VMM), which was initially introduced by Hughes et al.[4, 5] . The VMM has been used to simulate strain localization of both single and multiple dimensions[6–8] and crack propagation problems[9, 10] . Based on the LATIN method, a multiscale XFEM for crack propagation was proposed by Guidault et al.[11] . The method enables one to use a refined mesh only in the vicinity of discontinuities. In Ref.[12], a multiscale projection method was developed to simulate the effects of both macrocracks and microcracks. Belytschko et al.[13] developed a coarse-graining method to handle the multiscale crack propagation with the multiscale aggregating discontinuity method (MAD). There are also some other multiscale methods developed for discontinuous problems, to which the readers are referred[14–18] for more information. Recently, a new multiscale method, which is called the extended multiscale finite element method (EMsFEM), was developed by Zhang et al.[19, 20] for solving solid mechanical problems of multiple scale features. The main idea of the method is to numerically construct the multiscale base functions that can efficiently capture the small-scale properties within each coarse element by locally solving the Dirichlet boundary value problems. In order to consider the coupling effect on different directions in vector field problems, additional coupling terms are added to the multiscale base functions for interpolation of the displacement field. The important small-scale information is then gathered into a coarser scale through the coupling of the equivalent stiffness matrices. The method can be conveniently implemented and adopted to general problems without scale separation and periodicity assumptions. Moreover, the actual microscopic stress and strain within the coarse-scale elements can be obtained simultaneously along with the macroscopic displacement field by means of downscaling computation. The EMsFEM has been successfully used for elasto-plastic[19] and dynamic analyses[20] of heterogeneous materials without the restriction of scale separation and periodic assumptions. It should be noted that the idea of using numerical base functions had earlier been applied to solving elliptic problems arising from heat conductions in composite materials and flows in porous media (See Refs.[21, 22]). Most of these methods are devoted to solving scalar field problems. The simulation of crack propagation with the finite element method (FEM) has been a topic of growing interest in recent years. Various numerical methods have been developed to overcome the meshing difficulty faced when an evolving crack is modeled. Notable among them is the XFEM, which was first introduced by Belytschko et al.[23] . Then, a Heaviside step function in conjunction with twodimensional (2D) asymptotic crack-tip displacement fields are used to account for the crack in Ref.[24]. This allows the crack geometry to be represented independently of the finite element mesh, making remeshing unnecessary for model crack growth. Stolarska et al.[25] presented an algorithm that couples the LSM with the XFEM to solve the crack problem. The LSM was used to represent the crack location, including the location of crack tips, and the XFEM was used to compute the stress and displacement fields necessary for determining the rate of crack growth. The combination of the LSM and the XFEM seems to be quite elegant since the values of enrichment functions can be obtained by the level set functions. The LSM was further used for representing the location of holes and material interfaces by Sukumar et al.[26] . The XFEM has been widely used for modeling strong (displacement) and weak (strain) discontinuities in both 2D and 3D problems. Several surveys on the XFEM can be seen in Refs.[27–29]. Moreover, implementation of the XFEM for crack growth with the finite element codes can be found in Refs.[30–32]. In the present study, a new concurrent multiscale technique is proposed for crack propagation; whichs enables one to use a refined mesh only in the vicinity of discontinuities. For this purpose, the EMsFEM serves as a multiscale framework. During the loading process, the initial macroscale model is transformed into a multiscale model, in which the localized subregions (the vicinity of the discontinuities) are automatically detected and resolved by the XFEM on a fine-scale mesh. Meanwhile, the other subregions that surround the refined zone are resolved by the EMsFEM on a coarse-scale mesh. The displacements on the two different scales are then bridged at the scale-interface by virtue of the multiscale base functions. Level sets are employed here to describe the geometry of the discontinuities
Vol. 28, No. 3 Jingkai Wu et al.: A Concurrent Multiscale Method for Simulation of Crack Propagation · 237 ·
as well as to identify the fine-scale subregions. Note that the crack propagation is usually concentrated in a very small portion of the structure. Hence the dramatic reduction in computational cost.
II. EMsFEM AND ITS IMPLEMENTATION FOR ELASTIC PROBLEMS In this section, we briefly introduce the basic implementation process of the EMsFEM for elastic problems. The readers are referred to Refs.[19,20] for more details. For the EMsFEM, the main task is to numerically construct the multiscale base functions that can capture the small-scale properties within each coarse element by locally solving the Dirichlet boundary value problems. Then the smallscale information is gathered into large-scale through the equivalent stiffness matrices. Consider a 2D heterogeneous material structure occupying a region Ω and having a boundary Γ , as depicted in Fig.1. The equilibrium equation and boundary conditions can be expressed as
Fig. 1. Schematic description of the EMsFEM.
div (D : E (u)) = f σ·n=T u=u
in Ω on Γ σ on Γ u
(1)
where D is the fourth-order stiffness tensor representing material properties, E (u) is the strain tensor given as 1 T E (u) = ∇u + (∇u) (2) 2 and u is the displacement vector; σ is the stress vector which has three independent components, T T σ = σx σy τxy ; f is the body force vector, f = f x f y ; Γσ and Γu are the force and displacement boundary conditions, respectively; n is the unit normal vector. For the elastic problem, there are two computational steps that need to be performed in the EMsFEM. Firstly, the multiscale base functions of coarse elements (unit cells) are constructed numerically; then, upscaling and downscaling computations are implemented by virtue of the multiscale base functions. 2.1. The Construction of Multiscale Base Functions Take one of the coarse elements shown in Fig.1 as an example, which occupies a subregion Ω K , K Ω ⊂ Ω. The multiscale base functions of the element are constructed by solving the equilibrium equation in the subregion Ω K with specified boundary conditions. In refs.[19,20], several different kinds of boundary conditions have been proposed. From Eqs.(1) and (2), the general expression for solving the base functions of a 2D problem can be given as follows: LN ¯i = 0 in Ω k (i = 1, 2, 3, 4) (3) N ¯i (x, y) affined on ∂Ω k
· 238 ·
ACTA MECHANICA SOLIDA SINICA
2015
T where L is the elasticity operator and satisfies Lu = div D : 12 ∇u + (∇u) . Once the base functions are constructed, the microscopic displacement fields within a coarse element can be expressed as 4 4 X X u= N¯ixx u¯i + N¯ixy v¯i (4) i=1
v=
4 X
i=1
N¯iyy v¯i +
i=1
4 X
N¯iyx u¯i
(5)
i=1
where u and v are the microscopic displacement components in the x and y directions, respectively; N¯iyx (or N¯ixy ) is an additional coupling term indicating the displacement field in the y-direction (or x-direction) within the coarse element induced by unit displacement of macroscopic node ¯i in the x-direction (or y-direction). Equations (4) and (5) can be expressed in a unified form: u = Nu ¯E
(6)
where N is the matrix of multiscale base functions, u is the displacement vector of nodes on the fine-scale meshes within the unit cell, and u ¯E is the displacement vector of nodes on the coarse-scale level. They can be written as T u = u1 v1 u2 v2 · · · · · · un vn T T (7) N = R1 RT · · · RT 2 n T u ¯E = u¯1 v¯1 u¯2 v¯2 u¯3 v¯3 u¯4 v¯4
where
N¯1xx (i) N¯1xy (i) Ri = N¯1yx (i) N¯1yy (i) (i = 1, 2, ......, n)
N¯2xx (i) N¯2yx (i)
N¯2xy (i) N¯2yy (i)
N¯3xx (i) N¯3yx (i)
N¯3xy (i) N¯3yy (i)
N¯4xx (i) N¯4yx (i)
N¯4xy (i) N¯4yy (i)
(8)
and n is the total number of nodes on the fine-scale mesh within the unit cell. It should be mentioned that, for materials with periodic microstructures, the multiscale base functions need to be constructed only once on the unit cell and are used for all the coarse elements. 2.2. Upscaling and Downscaling Computations After the multiscale base functions are constructed, we can implement the upscaling and downscaling computations conveniently: the upscaling computation is adopted to obtain the macroscopic equivalent stiffness matrix of the coarse element and the downscaling computation to obtain the microscopic stress and strain information within the coarse element. The equivalent stiffness matrix of a coarse element can be obtained from the point of view of equivalence of strain energy, i.e., KE =
p X
K ′e ,
K ′e = GT e K e Ge
(9)
e=1
where p is the total number of fine-scale elements within the coarse element, and K e is the traditional stiffness matrix of a fine-scale element e (see Fig.1), which can be expressed as Z Ke = BT (10) e D e B e tdΩe Ωe
in which B e and De are the strain-displacement matrix and material property matrix of the element e, respectively; t is the element’s thickness; Ge is a transition matrix that denotes the mapping relation between the displacement vectors of the microscopic scale nodes and the macroscopic scale nodes: T ue = Ge u ¯E , Ge = Re1 Re2 Re3 Re4 (11)
Vol. 28, No. 3 Jingkai Wu et al.: A Concurrent Multiscale Method for Simulation of Crack Propagation · 239 ·
where ue is the displacement vector of the element e. On the other hand, the strain and stress of a fine-scale element (take the eth fine-scale element in the unit cell as an example) can be obtained by εe = B e u e = T e u ¯E ,
T e = B e Ge
(12)
S e = D e B e Ge
(13)
and σ e = D e εe = Se u ¯E ,
III. MODELING OF DISCONTINUITIES ON THE SMALL SCALE In this paper, the discontinuities induced by the material interfaces and crack are modeled on the small scale (fine scale). Consider a body Ω f ⊂ R2 that contains nm circular inclusions and an internal traction-free crack (see Fig.2). We denote the boundary of Ω f by Γ f , which is not only decomposed into displacements Γuf and traction Γσf boundaries, but also into material interfaces Γm = ∪i Γm,i (i = 1, · · · nm ) and crack surface Γc , i.e., Γ f = Γuf ∪ Γσf ∪ Γm ∪ Γc . The domains of inclusions are denoted as Ωin,i , in which Ωin,i ⊂ Ω f . Similar to Eq.(1), the equilibrium equations for the corresponding problem can be expressed as in Ω f div (D : E f (uf )) = f σf · n = T f on Γσf (14) uf = u ¯f on Γuf σf · n = 0 on Γc+ σf · n = 0 on Γc− where uf , σ f , E f are the displacement fields, stress vector, and strain vector on the fine scale, respectively. Both weak discontinuities and strong discontinuities are considered here. The solution to the equation set (14) is worked out using the XFEM. 3.1. BASIC CONCEPTION OF THE XFEM The traditional FEM will encounter meshing difficulty when handling moving discontinuity problems, since the mesh needs to conform with the geometry of the crack. The XFEM provides an efficient technique to overcome this problem. In this section, the XFEM for the mechanical description of the discontinuities is briefly discussed. In the XFEM, the discontinuous enrichment functions are added to the finite element displacement approximation to represent the discontinuities within an element by virtue of the partition of the unity method[33] . In Fig.3, a crack is located on a uniform mesh. The displacement approximations in a classical XFEM[24] can be written as follows: uf (x) =
X
I∈S
NI (x) u ˆI +
X
NJ (x) H(x)bJ +
J∈SH
Fig. 2. Microstructure of a body with inclusions and crack.
X
K∈SC
NK (x)
4 X
Fα (x) cαK
(15)
α=1
Fig. 3. Enriched nodes for 2D crack problems on a uniform mesh. The circled nodes (set of nodes SH ) are enriched by step function and the squared nodes (set of nodes SC ) are enriched by crack tip functions.
· 240 ·
ACTA MECHANICA SOLIDA SINICA
2015
where S is the set of all nodes of the mesh; SH ⊂ S is the set of nodes whose support is bisected by a crack; and SC ⊂ S is the set of nodes whose support contains the crack tips. u ˆI is the classical degrees of freedom at node I; bJ and cαK are the enriched degrees of freedom at the corresponding enriched nodes; H(x) is the generalized Heaviside function whose value is set to +1 above the crack and −1 below the crack; Fα (x) is a function that models the near tip asymptotic field: √ √ √ √ θ θ θ θ {Fα (x)} = (16) r sin , r cos , r sin sin θ, r cos sin θ 2 2 2 2 where r and θ are polar coordinates of the local axes at the tip of crack (xtip , ytip ). The approximation in Eq.(15) will give rise to the emergence of blending elements between the enriched and unenriched elements. In these blending elements, only part of their nodes are enriched, so the approximation does not satisfy the partition of unity. This can lead to a decrease in the rate of convergence and errors in the blending elements. Several modifications have been proposed to circumvent the blending problem (see Refs.[34–36]). For example, the Heaviside enrichment function in the standard XFEM can be replaced with the shifted enrichment function, which can be written as ψJ (x) = NJ (x) (H(x) − H(xJ ))
(17)
Thus, the interpolation property is satisfied since the shifted enrichment function for the discontinuity vanishes at the edges of elements. Also, the enrichment function F can be expressed as X X F (x) = |φI | NI (x) − φI NI (x) (18) I
I
where φI is the level set value at the node I for the material interface level set function (see Ref.[36]). For inclusion problems, the following form can be taken to represent a material interface in the displacement approximation[37]: X X uf (x) = NI (x) u ˆI + aL NL (x) F (x) (19) I∈S
L∈SI
where SI ⊂ S is the set of nodes whose support is split by the interface, and aL are the enriched degrees of freedom at such nodes. 3.2. Geometrical Description of Discontinuities The LSM proposed by Osher and Sethian[38] is used for geometrical description of the discontinuities. Consider a region Ω which is split into two subdomains ΩA and ΩB by the interface Γ . The main idea of the LSM is to represent the interface Γ by means of a scalar function φ(x), which has the following properties: φ(x(t), t) < 0 for x ∈ ΩA φ(x(t), t) > 0 for x ∈ ΩB (20) φ(x(t), t) = 0 for x ∈ Γ Therefore, a moving interface at a given time t can be formulated as Γ (t) = x ∈ R2 |φ(x, t) = 0
(21)
The evolution of the interface can then be governed by solving the Hamilton-Jacobi equation ∂φ(x) + vn (x) |∇φ(x)| = 0 ∂t
(22)
where vn (x) is the velocity of the level set field. In the case of circular inclusions shown in Fig.2, the geometry of the fixed material interfaces Γm = ∪j Γm,j can be described by the following level set function: φm (x) = min {kx − xm,i k − rm,i } | {z } φm,i
for
i = 1, · · · ,nm
(23)
Vol. 28, No. 3 Jingkai Wu et al.: A Concurrent Multiscale Method for Simulation of Crack Propagation · 241 ·
where xm,i is the location of the center of the ith inclusion whose radius is rm,i . In order to track the growth of a 2-D crack , two level set functions φ and ψ are needed, one for the crack’s surface (φ) and the other for the crack tip (ψ). Thus, the crack’s position Γ (t) can be defined by means of the updated values of φ and ψ: Γ (t) := x ∈ R2 |φ(x, t) = 0 and ψ(x, t) ≤ 0 (24)
In general, the scalar values of the level set function are computed only at the finite element nodes, and the level set function can be interpolated by the finite element shape functions X φ(x) = Ni (x)φi (25) i∈I
Thus, the true locations of the discontinuities are approximated by the iso-zero of the level set function φ. Generally, the mesh should be sufficiently fine in the vicinity of discontinuities for their precise location and reduction in numerical errors (see Ref.[37]). As an example, we consider a simple structure with an inclusion shown in Fig.4. A coarse and a fine mesh each are used to model the material interface. Figure 5 shows the level set function for the coarse mesh. We can see that the geometrical representation is poor in this case. At the same time, by use of the fine mesh (see Fig.6), a much better geometrical description for the interface is obtained by the level set function. An adaptive mesh refinement strategy will be introduced in the next section to obtain a finer mesh only in the vicinity of discontinuities.
Fig. 4 Geometry of a body with an inclusion.
Fig. 5. Level set functions in 3D and 2D axes on a coarse-scale mesh.
Fig. 6. Level set functions in 3D and 2D axes on a fine-scale mesh.
· 242 ·
ACTA MECHANICA SOLIDA SINICA
2015
Fig. 7. Illustration of the CMM for crack propagation.
IV. CONCURRENT MULTISCALE METHOD FOR MODELING DISCONTINUITIES In engineering practice, the initial element size of the structure is always sufficient to evaluate the global structural response, but it is not fine enough to represent the discontinuities induced by the fixed material interfaces or the evolving cracks. As pointed out in §3.2, a mesh refinement treatment is required in the vicinity of the discontinuities in order to precisely represent their geometry while accurately resolving the displacement field. On the other hand, engineers prefer not to change the initial finite element mesh of the structure during the evolution of discontinuities since the operation is extremely difficult. In this section, a concurrent multiscale method (CMM) is developed for failure analysis including both weak and strong discontinuities. During the loading process, the entire domain Ω of structure is automatically divided into two different regions, i.e., the coarse-scale discrete region Ω c and the fine-scale discrete region Ω f (see Fig.7). The two regions are solved independently whereby the discontinuities are taken into account on the fine-scale mesh. For this purpose, the EMsFEM is here employed to capture the macroscopic structural response to the coarse-scale region with moderate strains, whereas the XFEM is used to model the material interfaces and crack propagation on the fine-scale region. In this context, the discontinuities are solved by detailed microstructural simulation. The geometry of the material interfaces and the evolving cracks are described on the fine scale by means of the LSM. Then the displacements between the two scales can be bridged at the scale-interface Γ cf using the multiscale base functions. The submodels with two different scales are then combined into one numerical model, the so-called mixed structure model. Note that the discontinuities are usually located in a very small portion of the structure, so that the expensive calculation of detailed fine-scale simulation is limited to a minimum. 4.1. Discretization of Domain An adaptive procedure based on the distance of the element to the discontinuities is used to switch between the two levels of discretization. Take Fig.7 as an example. An evolving crack is located on the initial coarse mesh. Prior to each growth step, the distance of each macroscopic node to the crack tip d1 (x, t) is computed. If the value of d1 (x, t) at any one of the nodes in a coarse element is less than the prescribed critical value R1 , then the corresponding coarse element becomes a fine-scale subregion, and vice versa. The fine-scale subregion can be retransformed into a coarse element when all the values of d1 (x, t) at the four macroscopic nodes exceed R1 . The value of R1 can be set to be approximately 1.5 times the coarse element size. Note that in our method the crack is modeled only on the fine scale; thus the fine-scale subregions have to be adopted everywhere along the crack path. Moreover, in order to ensure that the enriched nodes are not present on the scale-interface Γ cf , a coarse element close to the crack should be replaced by a fine-scale subregion if the distance to the crack surface d2 (x, t) at any one of the macroscopic nodes in the coarse element is less than the prescribed critical value l1 . Here, the value of l1 is approximately 2 times the fine element size. Fortunately, finding that the values of
Vol. 28, No. 3 Jingkai Wu et al.: A Concurrent Multiscale Method for Simulation of Crack Propagation · 243 ·
d1 (x, t) and d2 (x, t) at each coarse-scale node can be obtained by virtue of the level set functions, we have p d1 (x, t) = ψ 2 (x, t) + φ2 (x, t), d2 (x, t) = |φ(x, t)| (26) For the fixed material interface, the same technique can be used to identify the fine-scale subregion. The distances of nodes to the material interface d3 can be obtained by d3 (x) = |φm (x)|
(27)
Fig. 8. Illustration of the CMM for a material interface.
If d3 (x) at any one of the nodes in a coarse element is less than the prescribed critical value l2 , or the element is bisected by the material interface, then the corresponding coarse element becomes a fine-scale subregion. In addition, l2 is equal to l1 in this paper. By taking the structure in Fig.4 as an example, the entire domain Ω is divided into Ω c and Ω f , as shown in Fig.8. In this case, the level set function for the inclusion is plotted in Fig.9. We can see that a good geometrical representation for the interface can be obtained. It should be mentioned that the distance based criterion pro- Fig. 9 Level set functions in 2D axes on a posed above could be replaced by a more sophisticated error es- multiscale mesh. timator in our multiscale scheme; the corresponding work is still going on. 4.2. Bridging between Two Different Scales As mentioned above, two levels of discretization are used in the CMM for the implementation of the EMsFEM and XFEM, respectively. Thus, tying relations should be employed to guarantee compatibility of the displacements at the scale-interface Γ cf between the two regions. From the viewpoint of the EMsFEM we know that the relationship between the microscopic and macroscopic displacement fields within a coarse element can be established by means of the multiscale base functions. Take advantage of this feature and the master-slave constraints can be adopted here to bridge the degrees of freedom between the two different scales. As shown in Fig.10, the microscopic scale nodes at the scale-interface are defined as slave nodes (for example, node i), while the nodes on the adjacent coarse element are defined as master nodes (such as nodes I, J, K, L). The relation of displacements between the salve and master nodes can be expressed as uf (i) = T i u ¯E(I,J,K,L) (28) where uf(i) is the displacement vector of the slave node i, u ¯E(I,J,K,L) is the displacement vector of the corresponding master nodes, T i is the control matrix of the node i and can be deduced from Eq.(11).
· 244 ·
ACTA MECHANICA SOLIDA SINICA
2015
Fig. 10. Master-slave constraints.
4.3. Solution Procedures The direction of crack propagation depends on the chosen physical criteria. In this paper, the maximum circumferential stress[39] is used to evolve a crack. The angle of crack growth can be expressed as s 2 1 K I − sign(KII ) 8 + KI (29) θc = 2 arctan 4 KII KII
where KI and KII are the mixed-mode stress intensity factors and can be calculated by virtue of the interaction integrals[24] . The constant crack growth increment is used to determine the magnitude of the incremental crack growth △a. The solution procedures for modeling quasi-static crack growth in the CMM can be summarized as follows: 1. Preprocessing: Construct the base functions N of the unit cell with the special boundary conditions, and derive the equivalent stiffness matrix K E of the coarse element. 2. According to the current discontinuities’ location, divide the entire domain Ω into the coarse-scale discrete region Ω c and the fine-scale discrete region Ω f . 3. Obtain the stiffness matrix K Ω c of the region Ω c and the stiffness matrix K Ω f of the region Ω f by use of the EMsFEM and XFEM, respectively. Bridge the degrees of freedom between the two different scales at the scale-interface Γ cf . Then, obtain the global stiffness matrix K mix of the mixed structure. 4. Input the displacements and traction boundary conditions. Solution for the macroscopic and microscopic (including the traditional and enriched) degrees of freedom. 5. Calculate KI , KII , θc and △a. Get the new location of crack tip. 6. Repeat the steps 2-5.
V. NUMERICAL EXAMPLES In this section, several representative numerical examples are presented to examine the validity of the proposed method. The various parameters influencing the results of the CMM are studied. Note that plane strain and linear elastic materials are assumed, and all numerical examples are dimensionless. Example 1 In the first numerical example, we investigate the influence of the different parameters (including the resolution of the coarse-scale mesh and the critical value R1 for the mesh refinement of region around the crack tip) on the result of the CMM. As shown in Fig.11, a structure with an edge crack of length 0.5a is loaded by prescribing displacement boundary conditions of mode I along all edges. The material
Vol. 28, No. 3 Jingkai Wu et al.: A Concurrent Multiscale Method for Simulation of Crack Propagation · 245 ·
Fig. 11. Sketch of a structure with edge crack.
parameters are chosen as follows: Young’s modulus E = 1.0 × 106 and the Poisson’s ratio ν = 0.3. The analytical expressions of displacement field (see the work of Griffith) can be written as r r θ KI θ KI r r cos (κ − cos θ), uy (r, θ) = sin (κ − cos θ) (30) ux (r, θ) = 2µ 2π 2 2µ 2π 2 where (r, θ) is the local polar coordinates, µ is the shear modulus and κ is Kolossov’s constant. The stress intensity factor is set to KI = 1.0. The coarse-scale meshes are composed of 9 × 9, 15 × 15, 27 × 27 and 45 × 45 elements, respectively. Moreover, the critical value R1 is set to R1 = 1.5 × lc and R1 = 3.0 × lc , respectively, where lc is the size of coarse-scale mesh. To eliminate the influence of the fine-scale resolution, the size of fine-scale mesh is chosen to be a constant whose value is h = a/405. Figure 12 shows the convergence of the stress intensity factor KI of the crack tip with respect to the number of coarse-scale elements. We can see that the results obtained by both cases can converge toward the theoretical value (1.0) when the number of coarse-scale elements is increasing, especially for the first increasing stage. This means that the number of coarse-scale elements 15 × 15 is sufficient for this test and can yield satisfactory results. It also shows that the accuracy can be obviously improved when a larger value is assigned to the critical value R1 . This is because the displacement gradients vary abruptly in the vicinity of the crack tip, and the fine-scale region should be large enough to simulate this phenomenon. However, more computational degrees of freedom are needed when the value R1 is larger. It should be mentioned that the results shown in Fig.12 are good even for the coarsest mesh resolution (deviation in only the third
Fig. 12. Convergence of stress intensity factor KI .
Fig. 13. The geometry of an edge crack problem.
· 246 ·
ACTA MECHANICA SOLIDA SINICA
2015
Fig. 14. Values of KI with respect to crack growth iteration.
Fig. 15. Geometry of a plate with an inclusion and a void.
digit); thus our multiscale method has a good accuracy. Example 2 Consider a finite plate with an edge crack of length a shown in Fig.13. The width and height of the plate are W and L, respectively. The following values are used, a = 2.5, W = 10, L = 19 and σ = 1. The theoretical Mode I stress intensity factor KI for this problem is given by[24] √ KI = F (λ) σ aπ
(31)
where F (λ) is a finite-geometry correction factor: F (λ) = 1.12 − 0.231λ + 10.55λ2 − 21.72λ3 + 30.39λ4 ,
λ=
a W
(32)
In each iteration, the crack is propagated in constant increment of △a = 0.25, and 15 iterative steps are considered. The number of the coarse-scale elements is kept constant at 10 × 19. In the fine-scale region, each of the coarse-scale elements is discretized by 9 × 9. The material parameters are the same as in Example 1. Figure 14 shows the values of the stress intensity factor KI obtained by the CMM and theoretical expression. It can be seen that the results obtained by the CMM fit fairly well with the exact solutions. Example 3 In this example, we study the problem of crack propagation in the presence of both a hard inclusion and a void. For the finite plate in the preceding example, here a circular inclusion with E = 5.0 × 106 , the Poisson’s ratio ν = 0.3 and radius R = 1.7 plus a circular void with radius R = 1.7 are located in the plate. The geometry of the plate is shown in Fig.15, where L1 = L2 = 5. The crack is propagated for 20 steps. For comparison, a full fine-scale simulation with 90 × 171 elements is accomplished, whose results can be seen as the reference solutions and are denoted as XFEM-F. Note that the size of elements used in the XFEM-F is the same as that of fine-scale elements in the CMM. The crack paths obtained by the two methods are compared in Fig.16, which are in excellent agreement. The crack propagates, as to be expected, close to the void by the influence of the inclusion and void. The mesh is refined only on the subregions where the discontinuities occur. Figure 17 shows the values of stress intensity factors KI and KII , respectively. Moreover, the distributions of the Von-Mises stress calculated by the two methods at three selected steps are illustrated in Fig.18. We can see that the results obtained by the CMM fit well with the reference solutions. For the computing times, the CMM requires about one-fifth the CPU time of the XFEM-F.
Vol. 28, No. 3 Jingkai Wu et al.: A Concurrent Multiscale Method for Simulation of Crack Propagation · 247 ·
Fig. 16. Crack paths at different iterative steps.
Fig. 17. Values of KI and KII with respect to crack growth iterations.
Fig. 18. Distributions of Von-Mises stress at different iterative steps.
Example 4 Let us consider the crack propagation of a three-point bending model, which is loaded by a single force F = 1 imposed at the center of the upper edge. A crack of initial length 0.6 is located on the lower edge. The material parameters are the same as in Example 1. The geometry and boundary conditions of the model are shown in Fig.19. In this study, a total of 20 iterative steps are considered for crack propagation, and the propagation increment is set to △a = 0.3. To test the mesh sensitivity of the multiscale method, two kinds of coarse-scale FE discretization of different element sizes are considered, in which the numbers of the coarse-scale elements are 10×4 and 20×8, respectively. The results obtained by the CMM for these two cases are denoted as CMM1 and CMM2, respectively. In the fine-scale region, square elements with length 1/11 are used. As is known that the stress concentration phenomenon will
· 248 ·
ACTA MECHANICA SOLIDA SINICA
2015
Fig. 19. Geometry and boundary conditions of the model.
emerge at the place where the forces and constraints are. Thus the meshes in the CMM are refined a priori in these subregions to reduce the errors. For the sake of comparison, a full fine-scale simulation with 220 × 88 elements is implemented by the XFEM. Figure 20 shows the crack paths for the XFEM-F, CMM1, and CMM2 at several selected steps. In addition, the deformed meshes for the three methods at the last step are plotted in Fig.21. The results obtained by both the CMM1 and CMM2 are in good agreement with those of the XFEM-F. The crack propagates in a curved path towards the centre of the upper edge. The values of stress intensity factors KI and KII as a function of the crack iterative steps are shown in Fig.22. Moreover, the distributions of Von-Mises stress at several selected steps for the three methods are plotted in Fig.23. The results further illustrate that the CMM proposed have a good accuracy. The method shows insensitivity to the coarse-scale discretization. The mesh refinement in the CMM needs only to be implemented in the vicinity of the crack, thus dramatically reducing the computer memory requirements and CPU time. For this example, the CMM2 requires about one-tenth the CPU time of the XFEM-F.
Fig. 20. Crack paths at different iterative steps.
Vol. 28, No. 3 Jingkai Wu et al.: A Concurrent Multiscale Method for Simulation of Crack Propagation · 249 ·
Fig. 21. Deformed meshes at the last iterative step.
Fig. 22. Values of KI and KII with respect to crack growth iterations.
Fig. 23. Distributions of Von-Mises stress at different iterative steps.
VI. CONCLUSIONS An concurrent multiscale method is developed to model quasi-static crack propagation in elastic material. The EMsFEM serves as a multiscale framework and is extended to take into consideration material failure. For this purpose, the combination of the XFEM and LSM is successfully integrated into the multiscale framework to simulate the crack propagation on a fine scale. Level sets are employed to describe the geometry of discontinuities as well as to automatically identify the fine-scale subregions. Then, the degrees of freedom between the two different scales are bridged at the scale-interfaces by means of multiscale base functions. Since the failure processes are usually concentrated in a very small portion of the structure, the expensive calculation of the fine-scale simulation can be limited to a minimum. Our numerical examples have shown that the proposed method can be successfully used for simulating the evolving discontinuities and for dramatically reducing the computational cost. Thus, the method seems to have great potential for failure analysis of large-scale structures.
References [1] Babuˇska,I., Homogenization approach in engineering. In: Lions,J.L. and Glowinski,R. (Eds.), Computing Methods in Applied Science and Engineering, Lecture Note in Economics and Mathematical Systems. Berlin: Springer, 1976, 134: 137-153. [2] Guedes,J.M. and Kikuchi,N., Preprocessing and postprocessing for materials based on the homogenization method with adaptive finite element methods. Computer Methods in Applied Mechanics and Engineering,
· 250 ·
ACTA MECHANICA SOLIDA SINICA
2015
1990, 83: 143-198. [3] Fish,J., Shek,K., Pandheeradi,M. and Shephard,M., Computational plasticity for composite structures based on mathematical homogenization: theory and practice. Computer Methods in Applied Mechanics and Engineering, 1997, 148: 53-73. [4] Hughes,T.J.R., Multiscale phenomena: Green’s function, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized method. Computer Methods in Applied Mechanics and Engineering, 1995, 127: 387-401. [5] Hughes,T.J.R., Feijoo,G.R., Mazzei,L. and Quincy,J.B., The variational multiscale method-a paradigm for computational mechanics. Computer Methods in Applied Mechanics and Engineering, 1998, 166: 3-24. [6] Garikipati,K. and Hughes,T.J.R., A study of strain localization in a multiple scale framework–the onedimensional problem. Computer Methods in Applied Mechanics and Engineering, 2000, 159: 193-222. [7] Garikipati,K. and Hughes,T.J.R., A variational multiscale approach to strain localization-formulation for multidimensional problems. Computer Methods in Applied Mechanics and Engineering, 2000, 188: 39-60. [8] Yeon,J.H. and Youn,S.K., Meshfree analysis of softening elastoplastic solids using variational multiscale method. International Journal of Solids and Structures, 2005, 42: 4030-4057. [9] Mergheim,J., A variational multiscale method to model crack propagation at finite strains. International Journal for Numerical Methods in Engineering, 2009, 80: 269-289. [10] Hettich,T., Hund,A. and Ramm,E., Modeling of failure in composites by X-FEM and level sets within a multiscale framework. Computer Methods in Applied Mechanics and Engineering, 2008, 197: 414-424. [11] Guidault,P.A., Allix,O., Champaney,L. and Cornuault,C., A multiscale extended finite element method for crack propagation. Computer Methods in Applied Mechanics and Engineering, 2008, 197: 381-399. [12] Loehnert,S. and Belytschko,T., A multiscale projection method for macro/microcrack simulations. International Journal for Numerical Methods in Engineering, 2007, 71: 1466-1482. [13] Belytschko,T. and Song,J.H., Coarse-graining of multiscale crack propagation. International Journal for Numerical Methods in Engineering, 2010, 81: 537-563. [14] Pereira,J.P.A., Kim,D.J. and Duarte,C.A., A two-scale approach for the analysis of propagating threedimensional fractures. Computational Mechanics, 2012, 49: 99-121. [15] Pierres,E., Baietto,M.C. and Gravouil,A., A two-scale extended finite element method for modelling 3D crack growth with interfacial contact. Computer Methods in Applied Mechanics and Engineering, 2010, 199: 1165-1177. [16] Holl,M., Loehnert,S. and Wriggers,P., An adaptive multiscale method for crack propagation and crack coalescence. International Journal for Numerical Methods in Engineering, 2013, 93: 23-51. [17] Kim,D.J., Pereira,J. P. and Duarte,C.A., Analysis of three-dimensional fracture mechanics problems: a two-scale approach using coarse-generalized FEM meshes. International Journal for Numerical Methods in Engineering, 2010, 81: 335-365. [18] Lian,W.D., Legrain,G. and Cartraud,P., Image-based computational homogenization and localization: comparison between X-FEM/levelset and voxel-based approaches. Computational Mechanics, 2013, 51: 279-293. [19] Zhang,H.W., Wu,J.K. and Fu,Z.D., A new multiscale computational method for elasto-plastic analysis of heterogeneous materials. Computational Mechanics, 2012, 49: 149-169. [20] Zhang,H.W., Liu,H. and Wu,J.K., A uniform multiscale method for 2D static and dynamic analyses of heterogeneous materials. International Journal for Numerical Methods in Engineering, 2013, 93: 714-746. [21] Hou,T.Y. and Wu,X.H., A multiscale finite element method for elliptic problems in composite materials and porous media. Journal of Computational Physics, 1997, 134: 169-89. [22] Hou,T.Y., Multiscale modelling and computation of fluid flow. International Journal for Numerical Methods in Fluids, 2005, 47: 707-719. [23] Belytschko,T. and Black,T., Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering, 1999, 45: 601-620. [24] Moes,N., Dolbow,J. and Belytschko,T., A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 1999, 46: 131-150. [25] Stolarska,M., Chopp,D.L., Moes,N. and Belyschko,T., Modelling crack growth by level sets in the extended finite element method. International Journal for Numerical Methods in Engineering, 2001, 51: 943-960. [26] Sukumar,N., Chopp,D., Moes,N. and Belytschko,T., Modeling holes and inclusions by level sets in the extended finite-element method. Computer Methods in Applied Mechanics and Engineering, 2001, 190: 6183-6200. [27] Karihaloo,B.L. and Xiao,Q.Z., Modelling of stationary and growing cracks in FE framework without remeshing: a state-of-the-art review. Computers and Structures, 2003, 81: 119-129. [28] Abdelaziz,Y. and Hamouine,A., A survey of the extended finite element. Computers and structures, 2008, 86: 1141-1151.
Vol. 28, No. 3 Jingkai Wu et al.: A Concurrent Multiscale Method for Simulation of Crack Propagation · 251 · [29] Fries,T.P. and Belytschko,T., The extended/generalized finite element method: An overview of the method and its applications. International Journal for Numerical Methods in Engineering, 2010, 84: 253-304. [30] Pais,M., MATLAB Extended Finite Element (MXFEM) Code, Version 1.2, 2010. (Available from: www.matthewpais.com. [Accessed on April 3, 2012]). [31] Sukumar,N. and Prevost,J.H., Modeling quasi-static crack growth with the extended finite element method part i: computer implementation. International Journal of Solids and Structures, 2003, 40: 7513-7537. [32] Bordas,S., Nguyen,P.V., Dunant,C., Guidoum,A. and Nguyen-Dang,H., An extended finite element library. International Journal for Numerical Methods in Engineering, 2007, 71: 703-732. [33] Melenk,J.M. and Babuska,I., The partition of unity finite element method: basic theory and applications. Computer Methods in Applied Mechanics and Engineering, 1996, 139: 289-314. [34] Belytschko,T., Moes,N., Usui,S. and Parimi,C., Arbitrary discontinuities in finite elements. International Journal for Numerical Methods in Engineering, 2001, 50: 993-1013. [35] Legay,A., Wang,H.W. and Belytschko,T., Strong and weak arbitrary discontinuities in spectral finite elements. International Journal for Numerical Methods in Engineering, 2005, 64: 991-1008. [36] Fries,T.P., A corrected XFEM approximation without problems in blending elements. International Journal for Numerical Methods in Engineering, 2008, 75(5): 503-532. [37] Moes,N., Cloirec,M., Cartraud,P. and Remacle,J., A computational approach to handle complex microstructure geometries. Computer Methods in Applied Mechanics and Engineering, 2003, 192: 3163-3177. [38] Osher,S. and Sethian,J., Fronts propagating with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics, 1988, 79: 12-49. [39] Erdogan,F. and Sih,G.C., On the crack extension in plates under plane loading and transverse shear. Transactions of the ASME. Journal of Basic Engineering, 1963, 85: 519-525.