- Email: [email protected]
Multiscale analysis of interaction between macro crack and microdefects by using multiscale projection method
Accepted Manuscript Multiscale analysis of interaction between macro crack and microdefects by using multiscale projection method Guangzhong Liu, Dai Zhou, Yan Bao, Jin Ma, Zhaolong. Han PII: DOI: Reference:
S0167-8442(16)30382-2 http://dx.doi.org/10.1016/j.tafmec.2017.03.002 TAFMEC 1818
To appear in:
Theoretical and Applied Fracture Mechanics
Received Date: Revised Date: Accepted Date:
18 November 2016 22 February 2017 2 March 2017
Please cite this article as: G. Liu, D. Zhou, Y. Bao, J. Ma, Zhaolong. Han, Multiscale analysis of interaction between macro crack and microdefects by using multiscale projection method, Theoretical and Applied Fracture Mechanics (2017), doi: http://dx.doi.org/10.1016/j.tafmec.2017.03.002
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Multiscale analysis of interaction between macro crack and microdefects by using multiscale projection method a, b, c 1 Guangzhong Liu , Dai Zhou * , Yan Baoa, Jin Maa, Zhaolong. Hana a
a. Department of Civil Engineering, School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, No. 800, Dongchuan Road, Shanghai, 200240, China b. State Key Laboratory of Ocean Engineering, Shanghai Jiao Tong University, No. 800, Dongchuan Road, Shanghai, 200240, China c. Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration, No. 800, Dongchuan Road, Shanghai, 200240, China
Abstract: The presence of microdefects near a macro crack tip can induce either amplification or shielding of macro crack propagation. In the present work, microdefects (micro cracks, inclusions, voids) in vicinity of a solitary macro crack tip are numerical simulated by using multiscale projection method. In this method, macro crack and microdefects are efficiently simulated by two scale decomposition and parallelization. The discontinuities are taken into account in the framework of XFEM, and the influence of the microdefects is shown by the stress intensity factor (SIF) calculated at the macro crack tip. The domain contains a macro edge crack along with microdefects is modeled under mode I condition as well as mode II condition. The location of microdefects is varied to study the relation between the location and their effects on the SIF of macro crack. The influence of micro voids on the main crack tip is investigated thoroughly for the first time. The results obtained from these simulations are useful for accurately evaluating the residual strength, and deliberately suppress the main crack propagation. Keyword: multiscale projection method; XFEM; micro cracks; inclusions; voids 1. Introduction The evolution of localized phenomena such as cracks, voids, inclusions, shear bands have been recognized as main causes of material failure. Even in ductile materials, accumulation and development of such microdefects in the vicinity of a major crack tip can be observed [1]. The presence of microdefects near a solitary macro crack tip can greatly influence the propagation of the main crack. Each microdefect has either amplification or shielding effect on the macro crack, depending on the microdefect’s location, and geometry. Thus, the mechanism of main crack interacting with microdefects is an important problem in material strength analysis and fracture mechanics. Numerous analytical and numerical studies have been conducted to investigate the effect the micro cracks have on the macro crack. For analytical solutions, several studies using integral equations have been reported concerning interaction between a single macro crack and arrays of micro cracks.
*Correspondence to: Department of Civil Engineering, School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, No. 800, Dongchuan Road, Shanghai, 200240, China. Tel.:02134206195 E-mail address: [email protected] (Dai Zhou)
Due to the complexity of the integral equations, only approximate solutions of integral equations can be obtained for the special few cases under restrictive assumptions. Till today, analytical investigations on the influence of other kinds of micro defects (such as inclusions and voids) are quite limited and still under further investigation. A review of studies devoted to micro/macrocrack interaction problems is given in Tamuzs and Petrova [2]. For numerical analysis, two main different kinds of approaches can be distinguished. The first one is continuum approach, where the region containing microdefects were replaced by a homogeneous material with reduced stiffness. This approach is computational cost saving, but may lead to accuracy losing, as discussed in detail by Kachanov et al. [3], in that the presence of microdefects effect the macro crack not only by material softening, but also by local stress concentration. The second approach is discrete approach, in which the geometries of the microdefects are explicitly modelled, although the computational cost could be severely high due to the refined mesh. For example, the well-known FEM simulations of a semi-infinite crack interacting with microdefects obtained by Gong [4] and Meguid [5] were in good agreement with results calculated form the asymptotic formulas, in terms of the main crack SIF. Brencich and Carpinteri [6] numerically simulated the macro crack interacting with arrays of micro cracks by using displacement discontinuity boundary elements. It was found that, those micro cracks with tip inside the process zone ahead of the macro crack tip would amplify the SIF value of the macro crack. In recent years, various multiscale computational methods have been developed which offer great promise in modelling microdefects problem, in the sense that, considerably computational costs can be saved. In general, the multiscale method enables macro-micro scale decomposition and coupling, where the localized phenomena and the whole bulk are modelled by global and local meshes. When it comes to the interaction between different scales, the multiscale methods based on homogenization schemes [7-9] are only useful for the entire bulk of structure, but usually fail to predict the accurate stress field in the vicinity of microdefects. To overcome this difficulty, several concurrent multiscale methods based on decomposition technique were developed, such as Dirichlet projection method [10], the variational multiscale method (VMM) [11], domain decomposition methods [12], multiscale aggregating method [13], multiscale projection method [14]. Among them, the multiscale projection method proposed by S. Loehnert and T. Belytschko enables the parallel computations of the microscale and macroscale problems. In this method, the different scale decomposition is achieved by transition of field variables between different scales, where the effect of microfield in the macroscale formulation are emphasized. Therefore, this multiscale method is efficient and ideal for investigating interaction between macro crack and microdefects. The simulation of fracture problem with the finite element method (FEM) requires that, the mesh is conformal to each discontinuity and extremely refined, and special singular elements are created to capture the stress singularity. Recent years, a novel numerical method has been developed, known as extended finite element method (XFEM) [15], in which the discontinuities are modelled by adding enrichment functions to standard FEM approximate formulation, so that the mesh can be completely independent of discontinuity geometry. So far, the XFEM has been applied to a variety of materials and
discontinuities such as cracks [16], voids [17], inclusions [18], in both static and dynamic problems[19-21]. Alternative method is novel extended 4-node quadrilateral finite element method (XCQ4) proposed by Bui [22,23], which made full use of both consecutive-interpolation procedure and XFEM. The method has been successfully applied in dynamic and 3D problems [24,25]. Also, extended isogeometric analysis (XIGA) [19,20] is developed, which taking the advantages of high order NURBS basis functions and local enrichment methods. Till today, most of the research works mainly have been focused on the influence of micro cracks on a main crack, whereas the interaction between main crack and microscale inclusions, voids have not been given much attention. In our previous work [26], we modified the multiscale XFEM by employing corrected XFEM on the macroscale level to further save computational cost, the improvement of the method was demonstrated by microcracking examples. In present work, the interaction between microdefects and a solitary main crack was numerically investigated by using multiscale projection technique in conjugation with XFEM. The shielding or amplification effect of microdefects is evaluated in terms of the SIF of main crack tip, which is calculated by interaction integral method. A series of microdefects problems with different configurations are numerically simulated, the locations of the microdefects are altered to investigate the relation between the location and their effects. 2. The multiscale strategy and procedure In this section, we briefly introduce the basic implementation process of the multiscale projection method. For more details, the readers are referred to S. Loehnert [14]. For simplicity, only two scales problems are considered in the present paper, although the method can be applied to multiple scales. Consider a 2D material structure which contains a major crack, with multiple microdefects near the major crack tip. It is assumed that, the sizes of microdefects are much smaller than the size of macro crack. 2.1. Multiscale test functions Let the whole computational domain under consideration be Ω0, with a boundary ∂Ω0. A subdomain containing the macro crack tip and multiple microdefects is taken out for microscale analysis. The subdomain is denoted as Ω1, with a boundary ∂Ω1. The scale of macroscale crack is denoted as l0, and the scale of the microscale defects is denoted as l1. We assume that the geometry and response of the microscale defects is significantly smaller than that of macroscale crack, in that l0>> l1. It can be postulated that, along ∂Ω1, the fluctuations of the variables contributing to the existence of microstructure is negligible. According to the study made by Loehnert, the radius of Ω1 should be at least third times as big as the microdefects to ensure high accuracy.
Fig. 1. Multiscale domain with a main crack in presence of multiple microdefects In the subdomain Ω1, the displacement field is superposed as u1 u0 u 1 , where, the displacement u 0 and u 1 represent macro and micro displacement in the subdomain. The micro displacement field u 1 is the contribution of microdefects on the domain. When the size of microdefects is significantly smaller compared with the macro crack, the influence of microdefects on the domain Ω0\Ω1 is negligible, therefore we specify that u 1 0 in the domain Ω0\Ω1. For the whole domain Ω0, the weak form of governing equation can be obtained by multiplying equilibrium equation with macroscale test function η0,
div (u 0
0
u 1 ) 0 d
0
t 0d0 0 b 0d+ 0 t c D0 d
D
(1)
where, 0D denotes the macroscale discontinuity, tc is cohesive tension, the macroscale test function η0 consists of continuous part and discontinuous part:
0 ( x) C0 ( x) D0 ( x)
(2)
The macroscale test function η0 is continuous across all those micro defects within the subdomain. Thus, by the use of macroscale test function, the role of microscale discontinuities is eliminated from the coarse scale equilibrium equation. However, their effect is still included, because u 1 depends on the presence of microdefects. In the subdomain Ω1, the following equation can be obtained by multiplying equilibrium equation of subdomain with microscale test function, 1 1 1 div (u ) d b d 1
1
(3)
where, η1 is the microscale test function, which also consists of continuous part and discontinuous part. In subdomain, η1 accounts for both micro defects and main crack.
1 ( x) C1 ( x) 1D ( x)
(4)
The test function of both scale should be decomposable into continuous part and discontinuous part, in order to eliminate the role of microdefects in the right side of macroscale equilibrium equation. The serve this purpose, we employ XFEM to solve both scale problems. Because, XFEM models
discontinuities by adding enrichment functions into continuous finite element approximation, the enrichment functions serve as the discontinuous part of the test functions ηD(x), which will be elaborated in Section 3. 2.2. Multiscale projection technique In the multiscale projection method, the macroscale and microscale problem is solved with coarse scale mesh and fine scale mesh separately, as depicted in Fig. 2. The domain Ω0 is subdivided into macroscale elements, of which, those macroscale elements that constitutes the subdomain are called multiscale elements. Each multiscale element is subdivided into a mesh of microscale elements congruent to the coarse scale mesh.
Fig. 2. sketch of coarse scale mesh and fine scale mesh For the fine scale domain, the boundary conditions are u1 u0 on ∂Ω1
(5)
The main task of the multiscale projection is to enable transition of field variables between coarse and fine meshes, and meet the requirement of compatibility condition. In many multiscale methods in conjugation with FEM, the boundary condition is enforced by a curvilinear integral along the boundary ∂Ω1. However, in XFEM, since the enriched scheme at the coarse and fine scale meshes differ, the boundary condition above is enforced by a least-square projection on a specified domain 1 as:
( Nˆ uˆ
1 1
1
Nˆ 0uˆ 0 ) ( Nˆ 1 uˆ1 )dΩ 0
(6)
where, δu is the variational displacement. The domain 1 is chosen to be a strip of a width of one microscale element along the boundary ∂Ω1, as marked yellow in Fig. 2. N and u denote the enriched shape function and degree of freedom in the XFEM. The discretized form of Eq. (6) can be written as
P uˆ dˆ
(7)
Pij 1 Nˆ i ( Nˆ j )T d
(8)
dˆi 1 Nˆ i1 ( Nˆ 0uˆ 0 )d
(9)
The approximation of the coarse scale mesh may be relatively coarse, as it can not capture the multiple microscale discontinuities and a much coarser mesh. To avoid a stiff boundary condition on the microscale problem, the domain size of Ω1 should be adequately large, so that fluctuations of variable field due to the microstructure are relatively small on the boundary ∂Ω1.
2.3. Solution procedure The solution procedure of the multiscale method consists of a few sequential iterations, in which coarse scale and fine scale equation are solved by XFEM collaterally, as listed bellow. In the list, B denotes the enriched strain matrix in XFEM. The solution of coarse scale problem is corrected based on 0 the fine scale solution, by means of the internal force fint,k , in the following step 2. The present
multiscale method has excellent convergence behavior, such that merely two or three sequential iterative steps are needed. We restricted the present work to linear elastic material, even though the method can easily be extended to simulate plastic behavior by replacing the solution step 4 with a Newton iteration step. As for elastic-plastic XFEM, please refer to the literature [27]. 1.
Initialize variables
k 0 , uˆk0 0 , uˆ1k 0 and k0 0 , 1k 0 2.
Solve the macroscale problem 0 0 ˆk0 , with K 0 uˆk0 fext fint, k for u
Kij0 0 ( Bˆi0 )T D Bˆ 0j d , fi0,ext
0
0 0 T 0 1 Ni0td and fint, k 0 ( B ) (uk uk )d
Set uˆk01 uˆk0 uˆk0 3.
Based on uˆk01 , project the boundary conditions onto ∂Ω1 P uˆ1k 1 dˆk 1 for uˆ1k 1 on the boundary ∂Ω1 with Pij 1 Nˆ i1 ( Nˆ 1j )T d , dˆi , k 1 1 Nˆ i1 ( Nˆ 0 uˆ k01 )d .
4.
Solve the microscale problem 1 1 1 K 1 uˆ1k 1 fext, k 1 for uˆk 1 on entire microscale domain Ω with
1 1 1 ˆL,k 1 , and Kij1 0 ( Bˆi1 )T D Bˆ1j d , where L represents all the nodes on ∂Ω fext, k 1 K u
1
5.
If uˆ1k 1 uˆ1k tol1 or uˆk0 tol0 , then set k=k+1, and return to step 1, until the tolerances are reached.
3. Formulation of XFEM Based on partition unity conception, certain enriched test functions and enriched degree of freedom are added into classical finite element approximation, in order to model the discontinuities, such as cracks, inclusions, and voids. The generalized form of XFEM approximation for a 2D problem can be presented as:
uˆ( x ) N j ( x )u j N j ( x )g ( x )e j jA
(9)
jE
where, A is the set of all the nodes in the domain, Nj is the classical finite element shape function, and uj is the nodal displacement vector of all nodes. E is the set of enriched nodes which belonging to elements containing discontinuities or interfaces. g(x) is the enriched test functions introduced to capture the discontinuities, ej is the additional degree of freedom associated with E set of nodes. Different enriched test functions have been developed for cracks, inclusions, which will be detailed below. For simplicity of the notation, we define
N j N j ( x ) N ( x ) j g ( x ) u j u j
e j
T
(10)
T
(11)
3.1. XFEM approximation for cracks For modeling of 2D domain problem with cracks, the approximation takes the form as: uˆ( x ) N j ( x )u j jA
N j ( x)[ H ( x) H ( x j )]a j
jM
jI J
4
N j ( x ) R( x ) [ F ( x ) F ( x j )]bj
(12)
1
where, M is the set of nodes belonging to those elements which intersects with the crack; I is the set of nodes belonging to the elements containing the crack tip; J is the second layer nodes around the crack tip. All the nodal subsets is marked in Fig. 3. H(x) is the heaviside function used to model the discontinuity in displacement, which takes +1 on one side of the crack surface and -1 on the other side. aj are the nodal unknowns enriched on the M set of nodes. Fα(x) (α=1~4) are four crack tip enrichment test functions, used to incorporate the crack tip asymptotic displacement field. bj (α=1~4) are the nodal unknowns enriched on the I and J set of nodes.
Fig. 3. Enriched nodal subsets in XFEM for a crack Considering the local polar coordinates r and θ centered at the crack tip, these four branch functions can be written as: F ( x ), 1~4 r [cos
2
sin
2
cos sin 2
sin sin ] 2
(13)
R(x) is a ramp function introduced into the enrichment functions to cope with the blending element
problem, more details see the work of Fries [28]. R( x ) N j ( x )
(14)
jI
3.2. XFEM approximation for inclusions Consider a domain containing a circular inclusion, as marked with a red circle in Fig. 5. The approximation for 2D domain containing inclusions in XFEM can be written as uˆ( x ) N j ( x )u j N j ( x ) ( x )c j jA
(12)
jL
where, cj is the nodal unknowns enriched on the L set of nodes, which is marked in Fig. 4. γ(x) is the modified enrichment test function [29], defined as
( x) N j ( x) j ( x) jL
N j ( x) j ( x)
(13)
jL
where, φ(x) is a distance function defined by ( x ) x2 y 2 r , are the coordinates originated at the center of the circular inclusion.
Fig. 4. Enriched nodal subset in XFEM for void 3.3. XFEM approximation for voids Void enrichment function was first developed by Sukumar et al, by using a signed distance function φ(x). Unlike the use of enrichment functions presented before, the void enrichment function V(x) here does not require additional unknowns [30]. The approximation for voids can be written as uˆ( x ) V ( x ) N j ( x )u j
(14)
jA
The void enrichment function V(x) is a signed distance function, which takes +1 outside of the void, and takes 0 inside the void. In practice, integration is simply skipped where V(x)=0. Additionally, nodes whose supports are completely within the voids are considered as fixed DOF. 3.4. Integration schemes at discontinuities Due to the singularity of the enrichment functions, the ordinary gauss quadrature technique cannot
accurately produce integration result. Those elements containing interfaces are partitioned into several sub-quads or sub-triangles, as shown in Fig. 5. Within each sub-quads or sub-triangles, the integrands keep continuous, and standard gauss quadrature is used. Particularly, in the element containing crack tip, a singular mapping technique is used [31], where the gauss points of the each sub-triangles are mapped by those of a reference quadrangular cell, as the case of Fig. 5c.
b
a
c
Fig. 5. Division of enriched elements and gauss points distribution (a) Split element into quadrilaterals (b) Split element into triangles and (c) Tip element 3.5. Computation of stress Intensity Factors In this paper, an interaction integral approach is employed to compute the SIF at the tip of the main crack. In this method, J-integral is obtained by calculation of a domain based interaction integral. The interaction integral is given as q u (2) u (1) M (1,2) * ij(1) i ij(2) i W (1,2)1 j s dA, A x x 1 1 x j
(15)
where A* is the integral domain which contains the crack tip; qs is a smoothing function; W(1,2) is the interaction strain energy density. The superscripts 1 and 2 represent actual and auxiliary state respectively. Then, by intentionally selecting the auxiliary state as K(2) 1 , K
(2)
0 , and then K(2) 0 ,
K (2) 1 , the mixed SIFs can be obtained through the following equation: M (1,2)
2 ( K(1) K (2) K (1) K (2) ) E
(16)
4. Problem descriptions Loehnert and Belytschko[32] investigated the effect of large array of micro cracks by using single scale XFEM, and found out that, even large arrays of micro cracks have little additional effect beyond those which are in the immediate proximity of the crack. Therefore, the problems considered in the present paper are only closely located microdefects in pair, which make pure mode I or II problems. Four series of cases are considered, i.e. a solitary major edge crack along with symmetrically
located microscale cracks, circular inclusions under mode I condition, circular inclusions under mode Ⅱ condition and circular voids. Their locations are altered to explore the amplification or shielding effects on the SIF of major crack. A rectangular domain with the size 100 mm×200 mm of 2024-T4 aluminum alloy with E=73.1GPa and ν=0.33 is considered in the present study. The major crack located at the edge of the domain is with the length a=50mm. The material properties of the inclusions are taken as E=20GPa and ν=0.33. Uniform distributed traction is applied at the upper and bottom boundary of the domain along y-direction. Plane stress condition and linear elastic material are assumed. SIF is calculated at the tip of the major crack by interaction integral approach. 5. Results and discussions For demonstration of the present multiscale method, please refer to our previous work [26], in which, the present multiscale method was compared with one single XFEM, as well as analytical approach. In real metal material, during the fatigue life, multiple micro-cracks and micro-voids can be found randomly distributed around the crack tip. However, as superposition principle [2] stipulates, the SIF of the main crack can be evaluated as a sum of contributions of each micro-defects. (which is valid for distances between micro-defects are larger than half of their length). Therefore, in this section, only two symmetrically located microdefects are considered. The locations of the microdefects are altered to study the relation between the position of micro defects and their effects. 5.1. Case I: A major crack with symmetrically located micro cracks Two symmetrically located micro cracks are incorporated near a major edge crack tip, the dimension and boundary conditions are depicted at in Fig. 6. The geometric parameters that determine the effect of micro cracks, are micro crack orientation ψ, and the location of the micro cracks (r, θ). The distance from center of micro crack to macro crack tip r and the micro crack length 2c only effect the magnitude of the amplification or shielding. Therefore, we set r and 2c constantly at small values as r=4mm and 2c=2mm, and increase θ and ψ with increment of 15°.
Fig. 6. Geometric and boundary conditions for micro crack problem (unit: mm) As an example, the deformed geometry for the case of θ=60° and ψ=60° is shown in Fig. 7, by stepwise zooming. The Mises stress contour plot for the subdomain of case θ=30° and θ=120° is shown
in Fig. 8. By the comparison made in Fig. 8, it can be seen that, when the micro cracks is located in front of the major crack tip, the presence of micro cracks can increase the stress gradient, which is the case of θ=30°. On the contrary, the micro cracks behind the major crack tip can reduce the stress gradient.
a
b
c
Fig. 7. Deformed geometry of example case for micro crack problem
a
b
Fig. 8. Mises stress contour plot for case (a) θ=30° and (b) θ=120° (unit: MPa) The mode-I SIF of the major crack with micro flaws K MF was obtained for various cases of different θ and ψ. KI represents the SIF of the major crack tip in absence of micro crack. The influence factor KMF / K was plotted as contour θ and ψ shown in Fig. 9. In Fig. 9, the blue represent the
shielding zones where 1 , the red represent the amplification zones where 1 . The amplification-shielding zones are also influenced by the distance ρ=r/a. The contour result obtained by present numerical simulation reveals the same trend, as compared with semi-finite crack solution, which obtained by approximation solution of analytical formula. 150
1.1
120
θ
1.05 90 1 60 0.95 30 0
30
60
ψ
90
120
150
Fig. 9. Amplification (red) and shielding (blue) zones From Fig. 9, it can be observed that, the influence factor λ changes dramatically with variation of θ, when ψ=0°. When ψ approaches 90°, the presence of micro cracks has insignificant influence on the SIF of macro crack tip. The neutral location of the micro crack, which does not affect the SIF of macro crack is called neutral angle (θn) in the literature. Many researches predicted the neutral angle by approximation solution of analytical formula, and found out θn≈70°, for ψ=0°, and the value fluctuates depending on ρ=r/a. The result obtained by present study is θn≈76°, for ψ=0°, which is in general agreement with results by exact analytical methods [33]. 5.2. Case II: symmetrical micro inclusions under mode-I condition In this series of cases, two symmetrical micro inclusions near a major crack are considered, under mode-I condition. The boundary conditions and dimension sizes are same with those in case I, which are depicted in Fig. 10. The radii of two circular inclusions are kept constant as small as 1mm. The distance between major crack and micro inclusions is increased in step from 4mm to 6mm, and θ is increased from 0° to 150°, with increment of 30°.
Fig. 10. A major edge crack with micro inclusions under mode-I condition (unit: mm) The mode-I SIF is calculated at the tip of the major edge crack. In absence of micro flaws, the SIF is obtained as KI 338.2MPa mm . Fig. 11 shows the plots of the influence factor KMF / K for various θ and distances d. It can be observed from Fig. 11 that, the amplification effect reaches maximum at 30°, and the shielding effect reaches maximum at 120°. The neutral angle is estimated as θn≈80°, which is quite similar with the case of major crack along with two symmetrical micro cracks of ψ=0°. As expected, the distance only influences the magnitude of the amplification or shielding effect. d=4mm d=5mm d=6mm
1.08
λ
1.04 1.00 0.96 0
30
60
90
120
150
θ(°) Fig. 11. Mode-I SIF variation of a major edge crack with micro inclusions For examples, the Von Mises stress contour plots for the cases of θ=30° and θ=120° are presented in Fig. 12. Similar to the cases of micro cracks, when the micro inclusions are located in front of the major crack tip, the presence of micro inclusions can increase the stress gradient, and have amplification effect. It also can be observed that, the high stress is released within the micro inclusions, and stress concentration appears around the inclusions as expected.
b
a
Fig. 12. Von Mises stress plot for the subdomain of the case with (a) θ=30° and (b) θ=120° (unit: MPa) 5.3 Case III: symmetrical micro inclusions under mode-II condition In this case, two symmetrical micro inclusions are incorporated in a major edge crack mode-II problem. Fig. 13 depicts the geometric dimensions and boundary conditions for this series of problems. The material properties of the main plate and inclusions are the same as in case II. We considered the influence of distance, which are varied from 4mm to 6mm, and inclusion angle θ which increases from 0° to 150°.
Fig. 13. A major edge crack with micro inclusions under mode-II condition (unit: mm) The mode-II SIF is calculated at the tip of the major edge crack. In absence of the micro inclusions, the SIF of the major crack tip is obtained as KII 186.2MPa mm ． Fig. 14 shows the plots of the influence factor λ for various θ and distance d. When θ=60°, the influence factor λ reaches a polar point as, λ≈1. When θ ranges from 60° to 120°, the micro inclusions have amplification effect, which can be explained by the high stress zone shape under mode-II tearing condition. d=4mm d=5mm d=6mm
1.08
λ
1.04 1.00 0.96 0.92 0
30
60
90
θ(°)
120 150
Fig. 14. Mode-II SIF variation of a major edge crack with micro inclusions Fig. 15 provides the Mises stress contour plots for the cases of θ=30°, θ=60° and θ=150°. When θ=60°, the presence of the micro inclusions have minimal effect on the stress concentration of the major crack tip. When θ=150°, stress concentration appears around the circular inclusions, and the stress concentration of the major crack is relieved. When θ=30°, the stress concentration is enhanced.
a
b
c
Fig. 15. Von Mises stress plot for the subdomain of the case with (a) θ=30°, (b) θ=60° and (c) θ=150° (unit: MPa) 5.4 Case IV: edge crack in presence of symmetrical micro voids The effect of micro voids is quite like that of micro inclusions. The voids can be considered as inclusions with zero elastic modulus. The effects of the micro voids on the major crack tip have the same trend with those of micro inclusions, only with enhanced magnitude. Therefore, in this subsection is only provided the comparison of the influence factors. Both mode-I and mode-II problems are considered, the configuration is the same as that shown in Fig. 10 and Fig. 13. The deformed geometry is presented by stepwise zooming in Fig. 16 as an example.
b
a
c
Fig. 16. Deformed geometry of example case for micro voids problem In Fig. 17 is shown the comparison between the influence factors of micro inclusions and that of micro voids, with the same distance of d=4mm. As expected, the presences of micro voids have the
same kind of effect as that of micro inclusion, but with enhanced magnitude. The magnitude of the influence factor λ depends on elastic modulus of the inclusion and the distance d.
a inclusion void
1.1
λ
λ
1.1
b
1.0
1.0 0.9
inclusion void
0.9 0
30
60
90
θ(°)
120 150
0.8
0
30
60
90
θ(°)
120 150
Fig. 17. Comparison of influence factors under (a) mode-I configuration and (b) mode- II configuration In order to enhance the fatigue resistance of structures working with a main crack, symmetrical minor-voids of can be inserted in particular location to reduce the SIF of the main crack tip. Micro voids introduce smaller stress concentration in the vicinity of the defects than micro-cracks. According to the results present in this sub-section, symmetrical micro voids should be inserted at location as θ=60° for mode-I condition and θ=150° for mode- II condition. 6. Conclusions In the present work, the multiscale analysis of micro defects (minor cracks, inclusions, voids) in vicinity of a major crack tip is carried out by using multiscale XFEM method. The use of multiscale projection method makes it efficient to numerically simulate micro discontinuities with significantly small geometric dimension. The influences of the micro defects are evaluated by the SIF of major crack tip. Various configurations of micro defects are considered. On the basis of these simulations, the following conclusions can be drawn: (1) The present multiscale method has been found ideal for simulation of interaction between macrocrack and multiple micro defects in immediate proximity. (2) The SIF of the main crack tip is significantly affected due to the presence of micro defects, the magnitude of the effect the micro defects depend on their positions. (3) In general, those micro defects in front of the major crack tip have amplification effect, and behind the major crack tip have shielding effect. (4) For mode-I problem, the neutral angle that differentiate the amplification and shielding effect is found out to be θn≈70°~80°, which is in good agreement with that obtained by exact analytical method. (5) The relation between influence factor λ and the location of the micro defects are quite different for mode-Ⅰ problem and mode-Ⅱ problem. The failure of component is mainly decided on the SIF of main crack tip, which influenced by close micro defects. According to superposition principle, . The present study help to accurately predict the failure of crack propagation. Micro voids can be deliberately inserted in particular location
according to the results of present study, to reduce the SIF of the main crack tip, consequently enhance the fracture resistance of structures working with main crack. By the virtue of XFEM, the main crack propagation in presence of micro defects can also be easily simulated by the present method without any difficulty. The present method can be further extended to 3D cracking problem. Acknowledgements Support from National Natural Science Foundation of China (Nos. 51278297 and 51679139), the Major Program of the National Natural Science Foundation of China (No. 51490674), Research Program of Shanghai Leader Talent（No.20）and Doctoral Disciplinary Special Research Project of Chinese Ministry of Education (No.20130073110096) are acknowledged.
References: [1] A.P. Jivkov, Structure of micro-crack population and damage evolution in quasi-brittle media, THEOR APPL FRACT MEC. 70 (2014) 1-9. [2] V.P. Tamuzs, V.E. Petrova, On macrocrack-microdefect interaction, International Applied Mechanics. 38 (2002) 1157-1177. [3] M. Kachanov, J. Laures, Three-dimensional problems of strongly interacting arbitrarily located penny-shaped cracks, INT J FRACTURE. 41 (1989) 289-313. [4] S.X. Gong, On the formation of near-tip microcracking and associated toughening effects, ENG FRACT MECH. 50 (1995) 29-39. [5] S.A. Meguid, P.E. Gaultier, S.X. Gong, A comparison between analytical and finite element analysis of main crack-microcrack interaction, ENG FRACT MECH. 38 (1991) 451-465. [6] A. Brencich, A. Carpinteri, Interaction of a main crack with ordered distributions of microcracks: a numerical technique by displacement discontinuity boundary elements, INT J FRACTURE. 76 (1995) 373-389. [7] E. Svenning, M. Fagerstrom, F. Larsson, Computational homogenization of microfractured continua using weakly periodic boundary conditions, COMPUT METHOD APPL M. 299 (2016) 1-21. [8] G. Bhardwaj, I.V. Singh, Fatigue crack growth analysis of a homogeneous plate in the presence of multiple defects using extended isogeometric analysis, J BRAZ SOC MECH SCI. 37 (2015) 1065-1082. [9] S. Kumar, I.V. Singh, B.K. Mishra, A homogenized XFEM approach to simulate fatigue crack growth problems, COMPUT STRUCT. 150 (2015) 1-22. [10] J. Ortiz, A.P. Cisilino, J.L. Otegui, Effect of microcracking on the micromechanics of fatigue crack growth in austempered ductile iron, FATIGUE FRACT ENG M. 24 (2001) 591-605. [11] J. Mergheim, A variational multiscale method to model crack propagation at finite strains, INT J NUMER METH ENG. 80 (2009) 269-289. [12] X. Bian, Z. Li, G.E. Karniadakis, Multi-resolution flow simulations by smoothed particle hydrodynamics via domain decomposition, J COMPUT PHYS. 297 (2015) 132-155. [13] J. Song, Y. Yoon, Multiscale failure analysis with coarse-grained micro cracks and damage, THEOR APPL FRACT MEC. 72 (2014) 100-109.
[14] S. Loehnert, T. Belytschko, A multiscale projection method for macro/microcrack simulations, INT J NUMER METH ENG. 71 (2007) 1466-1482. [15] C. Ye, J. Shi, G.J. Cheng, An eXtended Finite Element Method (XFEM) study on the effect of reinforcing particles on the crack propagation behavior in a metal–matrix composite, INT J FATIGUE. 44 (2012) 151-156. [16] I.V. Singh, B.K. Mishra, S. Bhattacharya, R.U. Patil, The numerical simulation of fatigue crack growth using extended finite element method, INT J FATIGUE. 36 (2012) 109-119. [17] I.V. Singh, B.K. Mishra, S. Bhattacharya, XFEM simulation of cracks, holes and inclusions in functionally graded materials, INT J MECH MATER DES. 7 (2011) 199-218. [18] N. Sukumar, D.L. Chopp, N. Mo S, T. Belytschko, Modeling holes and inclusions by level sets in the extended finite-element method, COMPUT METHOD APPL M. 190 (2001) 6183-6200. [19] Q.B. Tinh, S. Hirose, C. Zhang, T. Rabczuk, C. Wu, T. Saitoh, J. Lei, Extended isogeometric analysis for dynamic fracture in multiphase piezoelectric/piezomagnetic composites, MECH MATER. 97 (2016) 135-163. [20] Q.B. Tinh, Extended isogeometric dynamic and static fracture analysis for cracks in piezoelectric materials using NURBS, COMPUT METHOD APPL M. 295 (2015) 470-509. [21] Q.B. Tinh, C. Zhang, Extended finite element simulation of stationary dynamic cracks in piezoelectric solids under impact loading, COMP MATER SCI. 62 (2012) 243-257. [22] Z. Kang, Q.B. Tinh, N. Du Dinh, T. Saitoh, S. Hirose, An extended consecutive-interpolation quadrilateral element (XCQ4) applied to linear elastic fracture mechanics, ACTA MECH. 226 (2015) 3991-4015. [23] Z. Kang, T.Q. Bui, T. Saitoh, S. Hirose, Quasi-static crack propagation simulation by an enhanced nodal gradient finite element with different enrichments, THEOR APPL FRACT MEC. (2017) 61-77. [24] Z. Wang, T. Yu, Q.B. Tinh, S. Tanaka, C. Zhang, S. Hirose, J.L. Curiel-Sosa, 3-D local mesh refinement XFEM with variable-node hexahedron elements for extraction of stress intensity factors of straight and curved planar cracks, COMPUT METHOD APPL M. 313 (2017) 375-405. [25] Z. Wang, T. Yu, Q.B. Tinh, A.T. Ngoc, T.H.L. Nguyen, D.D. Nguyen, D. Duc Hong, Numerical modeling of 3-D inclusions and voids by a novel adaptive XFEM, ADV ENG SOFTW. 102 (2016) 105-122. [26] G. Liu, D. Zhou, Y. Bao, J. Ma, Z. Han, Multiscale simulation of major crack/minor cracks interplay with the corrected XFEM, ARCH CIV MECH ENG. 17 (2017) 410-418. [27] G. Liu, D. Zhou, J. Ma, Z. Han, Numerical investigation of mixed-mode crack growth in ductile material using elastic–plastic XFEM, J BRAZ SOC MECH SCI. 38 (2016) 1689-1699. [28] T. Fries, A corrected XFEM approximation without problems in blending elements, INT J NUMER METH ENG. 75 (2008) 503-532. [29] N. Moës, M. Cloirec, P. Cartraud, J.F. Remacle, A computational approach to handle complex microstructure geometries, COMPUT METHOD APPL M. 192 (2003) 3163-3177. [30] M.J. Pais, F.A.C. Viana, N.H. Kim, Enabling high-order integration of fatigue crack growth with surrogate modeling, INT J FATIGUE. 43 (2012) 150-159. [31] G. Ventura, R. Gracie, T. Belytschko, Fast integration and weight function blending in the extended finite element method, INT J NUMER METH ENG. 77 (2009) 1-29. [32] S. Loehnert, T. Belytschko, Crack shielding and amplification due to multiple microcracks
interacting with a macrocrack, INT J FRACTURE. 145 (2007) 1-8. [33] S.X. Gong, Microcrack interaction with a finite main crack: an exact formulation, INT J FRACTURE. 66 (1994) R51-R56.
Highlights Multiscale numerical simulation of macro crack with micro defects problem; The simulation is conducted by using XFEM and multiscale projection method; The relation between the position of micro defects and their effects are revealed in detail.
S0167-8442(16)30382-2 http://dx.doi.org/10.1016/j.tafmec.2017.03.002 TAFMEC 1818
To appear in:
Theoretical and Applied Fracture Mechanics
Received Date: Revised Date: Accepted Date:
18 November 2016 22 February 2017 2 March 2017
Please cite this article as: G. Liu, D. Zhou, Y. Bao, J. Ma, Zhaolong. Han, Multiscale analysis of interaction between macro crack and microdefects by using multiscale projection method, Theoretical and Applied Fracture Mechanics (2017), doi: http://dx.doi.org/10.1016/j.tafmec.2017.03.002
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Multiscale analysis of interaction between macro crack and microdefects by using multiscale projection method a, b, c 1 Guangzhong Liu , Dai Zhou * , Yan Baoa, Jin Maa, Zhaolong. Hana a
a. Department of Civil Engineering, School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, No. 800, Dongchuan Road, Shanghai, 200240, China b. State Key Laboratory of Ocean Engineering, Shanghai Jiao Tong University, No. 800, Dongchuan Road, Shanghai, 200240, China c. Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration, No. 800, Dongchuan Road, Shanghai, 200240, China
Abstract: The presence of microdefects near a macro crack tip can induce either amplification or shielding of macro crack propagation. In the present work, microdefects (micro cracks, inclusions, voids) in vicinity of a solitary macro crack tip are numerical simulated by using multiscale projection method. In this method, macro crack and microdefects are efficiently simulated by two scale decomposition and parallelization. The discontinuities are taken into account in the framework of XFEM, and the influence of the microdefects is shown by the stress intensity factor (SIF) calculated at the macro crack tip. The domain contains a macro edge crack along with microdefects is modeled under mode I condition as well as mode II condition. The location of microdefects is varied to study the relation between the location and their effects on the SIF of macro crack. The influence of micro voids on the main crack tip is investigated thoroughly for the first time. The results obtained from these simulations are useful for accurately evaluating the residual strength, and deliberately suppress the main crack propagation. Keyword: multiscale projection method; XFEM; micro cracks; inclusions; voids 1. Introduction The evolution of localized phenomena such as cracks, voids, inclusions, shear bands have been recognized as main causes of material failure. Even in ductile materials, accumulation and development of such microdefects in the vicinity of a major crack tip can be observed [1]. The presence of microdefects near a solitary macro crack tip can greatly influence the propagation of the main crack. Each microdefect has either amplification or shielding effect on the macro crack, depending on the microdefect’s location, and geometry. Thus, the mechanism of main crack interacting with microdefects is an important problem in material strength analysis and fracture mechanics. Numerous analytical and numerical studies have been conducted to investigate the effect the micro cracks have on the macro crack. For analytical solutions, several studies using integral equations have been reported concerning interaction between a single macro crack and arrays of micro cracks.
*Correspondence to: Department of Civil Engineering, School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, No. 800, Dongchuan Road, Shanghai, 200240, China. Tel.:02134206195 E-mail address: [email protected] (Dai Zhou)
Due to the complexity of the integral equations, only approximate solutions of integral equations can be obtained for the special few cases under restrictive assumptions. Till today, analytical investigations on the influence of other kinds of micro defects (such as inclusions and voids) are quite limited and still under further investigation. A review of studies devoted to micro/macrocrack interaction problems is given in Tamuzs and Petrova [2]. For numerical analysis, two main different kinds of approaches can be distinguished. The first one is continuum approach, where the region containing microdefects were replaced by a homogeneous material with reduced stiffness. This approach is computational cost saving, but may lead to accuracy losing, as discussed in detail by Kachanov et al. [3], in that the presence of microdefects effect the macro crack not only by material softening, but also by local stress concentration. The second approach is discrete approach, in which the geometries of the microdefects are explicitly modelled, although the computational cost could be severely high due to the refined mesh. For example, the well-known FEM simulations of a semi-infinite crack interacting with microdefects obtained by Gong [4] and Meguid [5] were in good agreement with results calculated form the asymptotic formulas, in terms of the main crack SIF. Brencich and Carpinteri [6] numerically simulated the macro crack interacting with arrays of micro cracks by using displacement discontinuity boundary elements. It was found that, those micro cracks with tip inside the process zone ahead of the macro crack tip would amplify the SIF value of the macro crack. In recent years, various multiscale computational methods have been developed which offer great promise in modelling microdefects problem, in the sense that, considerably computational costs can be saved. In general, the multiscale method enables macro-micro scale decomposition and coupling, where the localized phenomena and the whole bulk are modelled by global and local meshes. When it comes to the interaction between different scales, the multiscale methods based on homogenization schemes [7-9] are only useful for the entire bulk of structure, but usually fail to predict the accurate stress field in the vicinity of microdefects. To overcome this difficulty, several concurrent multiscale methods based on decomposition technique were developed, such as Dirichlet projection method [10], the variational multiscale method (VMM) [11], domain decomposition methods [12], multiscale aggregating method [13], multiscale projection method [14]. Among them, the multiscale projection method proposed by S. Loehnert and T. Belytschko enables the parallel computations of the microscale and macroscale problems. In this method, the different scale decomposition is achieved by transition of field variables between different scales, where the effect of microfield in the macroscale formulation are emphasized. Therefore, this multiscale method is efficient and ideal for investigating interaction between macro crack and microdefects. The simulation of fracture problem with the finite element method (FEM) requires that, the mesh is conformal to each discontinuity and extremely refined, and special singular elements are created to capture the stress singularity. Recent years, a novel numerical method has been developed, known as extended finite element method (XFEM) [15], in which the discontinuities are modelled by adding enrichment functions to standard FEM approximate formulation, so that the mesh can be completely independent of discontinuity geometry. So far, the XFEM has been applied to a variety of materials and
discontinuities such as cracks [16], voids [17], inclusions [18], in both static and dynamic problems[19-21]. Alternative method is novel extended 4-node quadrilateral finite element method (XCQ4) proposed by Bui [22,23], which made full use of both consecutive-interpolation procedure and XFEM. The method has been successfully applied in dynamic and 3D problems [24,25]. Also, extended isogeometric analysis (XIGA) [19,20] is developed, which taking the advantages of high order NURBS basis functions and local enrichment methods. Till today, most of the research works mainly have been focused on the influence of micro cracks on a main crack, whereas the interaction between main crack and microscale inclusions, voids have not been given much attention. In our previous work [26], we modified the multiscale XFEM by employing corrected XFEM on the macroscale level to further save computational cost, the improvement of the method was demonstrated by microcracking examples. In present work, the interaction between microdefects and a solitary main crack was numerically investigated by using multiscale projection technique in conjugation with XFEM. The shielding or amplification effect of microdefects is evaluated in terms of the SIF of main crack tip, which is calculated by interaction integral method. A series of microdefects problems with different configurations are numerically simulated, the locations of the microdefects are altered to investigate the relation between the location and their effects. 2. The multiscale strategy and procedure In this section, we briefly introduce the basic implementation process of the multiscale projection method. For more details, the readers are referred to S. Loehnert [14]. For simplicity, only two scales problems are considered in the present paper, although the method can be applied to multiple scales. Consider a 2D material structure which contains a major crack, with multiple microdefects near the major crack tip. It is assumed that, the sizes of microdefects are much smaller than the size of macro crack. 2.1. Multiscale test functions Let the whole computational domain under consideration be Ω0, with a boundary ∂Ω0. A subdomain containing the macro crack tip and multiple microdefects is taken out for microscale analysis. The subdomain is denoted as Ω1, with a boundary ∂Ω1. The scale of macroscale crack is denoted as l0, and the scale of the microscale defects is denoted as l1. We assume that the geometry and response of the microscale defects is significantly smaller than that of macroscale crack, in that l0>> l1. It can be postulated that, along ∂Ω1, the fluctuations of the variables contributing to the existence of microstructure is negligible. According to the study made by Loehnert, the radius of Ω1 should be at least third times as big as the microdefects to ensure high accuracy.
Fig. 1. Multiscale domain with a main crack in presence of multiple microdefects In the subdomain Ω1, the displacement field is superposed as u1 u0 u 1 , where, the displacement u 0 and u 1 represent macro and micro displacement in the subdomain. The micro displacement field u 1 is the contribution of microdefects on the domain. When the size of microdefects is significantly smaller compared with the macro crack, the influence of microdefects on the domain Ω0\Ω1 is negligible, therefore we specify that u 1 0 in the domain Ω0\Ω1. For the whole domain Ω0, the weak form of governing equation can be obtained by multiplying equilibrium equation with macroscale test function η0,
div (u 0
0
u 1 ) 0 d
0
t 0d0 0 b 0d+ 0 t c D0 d
D
(1)
where, 0D denotes the macroscale discontinuity, tc is cohesive tension, the macroscale test function η0 consists of continuous part and discontinuous part:
0 ( x) C0 ( x) D0 ( x)
(2)
The macroscale test function η0 is continuous across all those micro defects within the subdomain. Thus, by the use of macroscale test function, the role of microscale discontinuities is eliminated from the coarse scale equilibrium equation. However, their effect is still included, because u 1 depends on the presence of microdefects. In the subdomain Ω1, the following equation can be obtained by multiplying equilibrium equation of subdomain with microscale test function, 1 1 1 div (u ) d b d 1
1
(3)
where, η1 is the microscale test function, which also consists of continuous part and discontinuous part. In subdomain, η1 accounts for both micro defects and main crack.
1 ( x) C1 ( x) 1D ( x)
(4)
The test function of both scale should be decomposable into continuous part and discontinuous part, in order to eliminate the role of microdefects in the right side of macroscale equilibrium equation. The serve this purpose, we employ XFEM to solve both scale problems. Because, XFEM models
discontinuities by adding enrichment functions into continuous finite element approximation, the enrichment functions serve as the discontinuous part of the test functions ηD(x), which will be elaborated in Section 3. 2.2. Multiscale projection technique In the multiscale projection method, the macroscale and microscale problem is solved with coarse scale mesh and fine scale mesh separately, as depicted in Fig. 2. The domain Ω0 is subdivided into macroscale elements, of which, those macroscale elements that constitutes the subdomain are called multiscale elements. Each multiscale element is subdivided into a mesh of microscale elements congruent to the coarse scale mesh.
Fig. 2. sketch of coarse scale mesh and fine scale mesh For the fine scale domain, the boundary conditions are u1 u0 on ∂Ω1
(5)
The main task of the multiscale projection is to enable transition of field variables between coarse and fine meshes, and meet the requirement of compatibility condition. In many multiscale methods in conjugation with FEM, the boundary condition is enforced by a curvilinear integral along the boundary ∂Ω1. However, in XFEM, since the enriched scheme at the coarse and fine scale meshes differ, the boundary condition above is enforced by a least-square projection on a specified domain 1 as:
( Nˆ uˆ
1 1
1
Nˆ 0uˆ 0 ) ( Nˆ 1 uˆ1 )dΩ 0
(6)
where, δu is the variational displacement. The domain 1 is chosen to be a strip of a width of one microscale element along the boundary ∂Ω1, as marked yellow in Fig. 2. N and u denote the enriched shape function and degree of freedom in the XFEM. The discretized form of Eq. (6) can be written as
P uˆ dˆ
(7)
Pij 1 Nˆ i ( Nˆ j )T d
(8)
dˆi 1 Nˆ i1 ( Nˆ 0uˆ 0 )d
(9)
The approximation of the coarse scale mesh may be relatively coarse, as it can not capture the multiple microscale discontinuities and a much coarser mesh. To avoid a stiff boundary condition on the microscale problem, the domain size of Ω1 should be adequately large, so that fluctuations of variable field due to the microstructure are relatively small on the boundary ∂Ω1.
2.3. Solution procedure The solution procedure of the multiscale method consists of a few sequential iterations, in which coarse scale and fine scale equation are solved by XFEM collaterally, as listed bellow. In the list, B denotes the enriched strain matrix in XFEM. The solution of coarse scale problem is corrected based on 0 the fine scale solution, by means of the internal force fint,k , in the following step 2. The present
multiscale method has excellent convergence behavior, such that merely two or three sequential iterative steps are needed. We restricted the present work to linear elastic material, even though the method can easily be extended to simulate plastic behavior by replacing the solution step 4 with a Newton iteration step. As for elastic-plastic XFEM, please refer to the literature [27]. 1.
Initialize variables
k 0 , uˆk0 0 , uˆ1k 0 and k0 0 , 1k 0 2.
Solve the macroscale problem 0 0 ˆk0 , with K 0 uˆk0 fext fint, k for u
Kij0 0 ( Bˆi0 )T D Bˆ 0j d , fi0,ext
0
0 0 T 0 1 Ni0td and fint, k 0 ( B ) (uk uk )d
Set uˆk01 uˆk0 uˆk0 3.
Based on uˆk01 , project the boundary conditions onto ∂Ω1 P uˆ1k 1 dˆk 1 for uˆ1k 1 on the boundary ∂Ω1 with Pij 1 Nˆ i1 ( Nˆ 1j )T d , dˆi , k 1 1 Nˆ i1 ( Nˆ 0 uˆ k01 )d .
4.
Solve the microscale problem 1 1 1 K 1 uˆ1k 1 fext, k 1 for uˆk 1 on entire microscale domain Ω with
1 1 1 ˆL,k 1 , and Kij1 0 ( Bˆi1 )T D Bˆ1j d , where L represents all the nodes on ∂Ω fext, k 1 K u
1
5.
If uˆ1k 1 uˆ1k tol1 or uˆk0 tol0 , then set k=k+1, and return to step 1, until the tolerances are reached.
3. Formulation of XFEM Based on partition unity conception, certain enriched test functions and enriched degree of freedom are added into classical finite element approximation, in order to model the discontinuities, such as cracks, inclusions, and voids. The generalized form of XFEM approximation for a 2D problem can be presented as:
uˆ( x ) N j ( x )u j N j ( x )g ( x )e j jA
(9)
jE
where, A is the set of all the nodes in the domain, Nj is the classical finite element shape function, and uj is the nodal displacement vector of all nodes. E is the set of enriched nodes which belonging to elements containing discontinuities or interfaces. g(x) is the enriched test functions introduced to capture the discontinuities, ej is the additional degree of freedom associated with E set of nodes. Different enriched test functions have been developed for cracks, inclusions, which will be detailed below. For simplicity of the notation, we define
N j N j ( x ) N ( x ) j g ( x ) u j u j
e j
T
(10)
T
(11)
3.1. XFEM approximation for cracks For modeling of 2D domain problem with cracks, the approximation takes the form as: uˆ( x ) N j ( x )u j jA
N j ( x)[ H ( x) H ( x j )]a j
jM
jI J
4
N j ( x ) R( x ) [ F ( x ) F ( x j )]bj
(12)
1
where, M is the set of nodes belonging to those elements which intersects with the crack; I is the set of nodes belonging to the elements containing the crack tip; J is the second layer nodes around the crack tip. All the nodal subsets is marked in Fig. 3. H(x) is the heaviside function used to model the discontinuity in displacement, which takes +1 on one side of the crack surface and -1 on the other side. aj are the nodal unknowns enriched on the M set of nodes. Fα(x) (α=1~4) are four crack tip enrichment test functions, used to incorporate the crack tip asymptotic displacement field. bj (α=1~4) are the nodal unknowns enriched on the I and J set of nodes.
Fig. 3. Enriched nodal subsets in XFEM for a crack Considering the local polar coordinates r and θ centered at the crack tip, these four branch functions can be written as: F ( x ), 1~4 r [cos
2
sin
2
cos sin 2
sin sin ] 2
(13)
R(x) is a ramp function introduced into the enrichment functions to cope with the blending element
problem, more details see the work of Fries [28]. R( x ) N j ( x )
(14)
jI
3.2. XFEM approximation for inclusions Consider a domain containing a circular inclusion, as marked with a red circle in Fig. 5. The approximation for 2D domain containing inclusions in XFEM can be written as uˆ( x ) N j ( x )u j N j ( x ) ( x )c j jA
(12)
jL
where, cj is the nodal unknowns enriched on the L set of nodes, which is marked in Fig. 4. γ(x) is the modified enrichment test function [29], defined as
( x) N j ( x) j ( x) jL
N j ( x) j ( x)
(13)
jL
where, φ(x) is a distance function defined by ( x ) x2 y 2 r , are the coordinates originated at the center of the circular inclusion.
Fig. 4. Enriched nodal subset in XFEM for void 3.3. XFEM approximation for voids Void enrichment function was first developed by Sukumar et al, by using a signed distance function φ(x). Unlike the use of enrichment functions presented before, the void enrichment function V(x) here does not require additional unknowns [30]. The approximation for voids can be written as uˆ( x ) V ( x ) N j ( x )u j
(14)
jA
The void enrichment function V(x) is a signed distance function, which takes +1 outside of the void, and takes 0 inside the void. In practice, integration is simply skipped where V(x)=0. Additionally, nodes whose supports are completely within the voids are considered as fixed DOF. 3.4. Integration schemes at discontinuities Due to the singularity of the enrichment functions, the ordinary gauss quadrature technique cannot
accurately produce integration result. Those elements containing interfaces are partitioned into several sub-quads or sub-triangles, as shown in Fig. 5. Within each sub-quads or sub-triangles, the integrands keep continuous, and standard gauss quadrature is used. Particularly, in the element containing crack tip, a singular mapping technique is used [31], where the gauss points of the each sub-triangles are mapped by those of a reference quadrangular cell, as the case of Fig. 5c.
b
a
c
Fig. 5. Division of enriched elements and gauss points distribution (a) Split element into quadrilaterals (b) Split element into triangles and (c) Tip element 3.5. Computation of stress Intensity Factors In this paper, an interaction integral approach is employed to compute the SIF at the tip of the main crack. In this method, J-integral is obtained by calculation of a domain based interaction integral. The interaction integral is given as q u (2) u (1) M (1,2) * ij(1) i ij(2) i W (1,2)1 j s dA, A x x 1 1 x j
(15)
where A* is the integral domain which contains the crack tip; qs is a smoothing function; W(1,2) is the interaction strain energy density. The superscripts 1 and 2 represent actual and auxiliary state respectively. Then, by intentionally selecting the auxiliary state as K(2) 1 , K
(2)
0 , and then K(2) 0 ,
K (2) 1 , the mixed SIFs can be obtained through the following equation: M (1,2)
2 ( K(1) K (2) K (1) K (2) ) E
(16)
4. Problem descriptions Loehnert and Belytschko[32] investigated the effect of large array of micro cracks by using single scale XFEM, and found out that, even large arrays of micro cracks have little additional effect beyond those which are in the immediate proximity of the crack. Therefore, the problems considered in the present paper are only closely located microdefects in pair, which make pure mode I or II problems. Four series of cases are considered, i.e. a solitary major edge crack along with symmetrically
located microscale cracks, circular inclusions under mode I condition, circular inclusions under mode Ⅱ condition and circular voids. Their locations are altered to explore the amplification or shielding effects on the SIF of major crack. A rectangular domain with the size 100 mm×200 mm of 2024-T4 aluminum alloy with E=73.1GPa and ν=0.33 is considered in the present study. The major crack located at the edge of the domain is with the length a=50mm. The material properties of the inclusions are taken as E=20GPa and ν=0.33. Uniform distributed traction is applied at the upper and bottom boundary of the domain along y-direction. Plane stress condition and linear elastic material are assumed. SIF is calculated at the tip of the major crack by interaction integral approach. 5. Results and discussions For demonstration of the present multiscale method, please refer to our previous work [26], in which, the present multiscale method was compared with one single XFEM, as well as analytical approach. In real metal material, during the fatigue life, multiple micro-cracks and micro-voids can be found randomly distributed around the crack tip. However, as superposition principle [2] stipulates, the SIF of the main crack can be evaluated as a sum of contributions of each micro-defects. (which is valid for distances between micro-defects are larger than half of their length). Therefore, in this section, only two symmetrically located microdefects are considered. The locations of the microdefects are altered to study the relation between the position of micro defects and their effects. 5.1. Case I: A major crack with symmetrically located micro cracks Two symmetrically located micro cracks are incorporated near a major edge crack tip, the dimension and boundary conditions are depicted at in Fig. 6. The geometric parameters that determine the effect of micro cracks, are micro crack orientation ψ, and the location of the micro cracks (r, θ). The distance from center of micro crack to macro crack tip r and the micro crack length 2c only effect the magnitude of the amplification or shielding. Therefore, we set r and 2c constantly at small values as r=4mm and 2c=2mm, and increase θ and ψ with increment of 15°.
Fig. 6. Geometric and boundary conditions for micro crack problem (unit: mm) As an example, the deformed geometry for the case of θ=60° and ψ=60° is shown in Fig. 7, by stepwise zooming. The Mises stress contour plot for the subdomain of case θ=30° and θ=120° is shown
in Fig. 8. By the comparison made in Fig. 8, it can be seen that, when the micro cracks is located in front of the major crack tip, the presence of micro cracks can increase the stress gradient, which is the case of θ=30°. On the contrary, the micro cracks behind the major crack tip can reduce the stress gradient.
a
b
c
Fig. 7. Deformed geometry of example case for micro crack problem
a
b
Fig. 8. Mises stress contour plot for case (a) θ=30° and (b) θ=120° (unit: MPa) The mode-I SIF of the major crack with micro flaws K MF was obtained for various cases of different θ and ψ. KI represents the SIF of the major crack tip in absence of micro crack. The influence factor KMF / K was plotted as contour θ and ψ shown in Fig. 9. In Fig. 9, the blue represent the
shielding zones where 1 , the red represent the amplification zones where 1 . The amplification-shielding zones are also influenced by the distance ρ=r/a. The contour result obtained by present numerical simulation reveals the same trend, as compared with semi-finite crack solution, which obtained by approximation solution of analytical formula. 150
1.1
120
θ
1.05 90 1 60 0.95 30 0
30
60
ψ
90
120
150
Fig. 9. Amplification (red) and shielding (blue) zones From Fig. 9, it can be observed that, the influence factor λ changes dramatically with variation of θ, when ψ=0°. When ψ approaches 90°, the presence of micro cracks has insignificant influence on the SIF of macro crack tip. The neutral location of the micro crack, which does not affect the SIF of macro crack is called neutral angle (θn) in the literature. Many researches predicted the neutral angle by approximation solution of analytical formula, and found out θn≈70°, for ψ=0°, and the value fluctuates depending on ρ=r/a. The result obtained by present study is θn≈76°, for ψ=0°, which is in general agreement with results by exact analytical methods [33]. 5.2. Case II: symmetrical micro inclusions under mode-I condition In this series of cases, two symmetrical micro inclusions near a major crack are considered, under mode-I condition. The boundary conditions and dimension sizes are same with those in case I, which are depicted in Fig. 10. The radii of two circular inclusions are kept constant as small as 1mm. The distance between major crack and micro inclusions is increased in step from 4mm to 6mm, and θ is increased from 0° to 150°, with increment of 30°.
Fig. 10. A major edge crack with micro inclusions under mode-I condition (unit: mm) The mode-I SIF is calculated at the tip of the major edge crack. In absence of micro flaws, the SIF is obtained as KI 338.2MPa mm . Fig. 11 shows the plots of the influence factor KMF / K for various θ and distances d. It can be observed from Fig. 11 that, the amplification effect reaches maximum at 30°, and the shielding effect reaches maximum at 120°. The neutral angle is estimated as θn≈80°, which is quite similar with the case of major crack along with two symmetrical micro cracks of ψ=0°. As expected, the distance only influences the magnitude of the amplification or shielding effect. d=4mm d=5mm d=6mm
1.08
λ
1.04 1.00 0.96 0
30
60
90
120
150
θ(°) Fig. 11. Mode-I SIF variation of a major edge crack with micro inclusions For examples, the Von Mises stress contour plots for the cases of θ=30° and θ=120° are presented in Fig. 12. Similar to the cases of micro cracks, when the micro inclusions are located in front of the major crack tip, the presence of micro inclusions can increase the stress gradient, and have amplification effect. It also can be observed that, the high stress is released within the micro inclusions, and stress concentration appears around the inclusions as expected.
b
a
Fig. 12. Von Mises stress plot for the subdomain of the case with (a) θ=30° and (b) θ=120° (unit: MPa) 5.3 Case III: symmetrical micro inclusions under mode-II condition In this case, two symmetrical micro inclusions are incorporated in a major edge crack mode-II problem. Fig. 13 depicts the geometric dimensions and boundary conditions for this series of problems. The material properties of the main plate and inclusions are the same as in case II. We considered the influence of distance, which are varied from 4mm to 6mm, and inclusion angle θ which increases from 0° to 150°.
Fig. 13. A major edge crack with micro inclusions under mode-II condition (unit: mm) The mode-II SIF is calculated at the tip of the major edge crack. In absence of the micro inclusions, the SIF of the major crack tip is obtained as KII 186.2MPa mm ． Fig. 14 shows the plots of the influence factor λ for various θ and distance d. When θ=60°, the influence factor λ reaches a polar point as, λ≈1. When θ ranges from 60° to 120°, the micro inclusions have amplification effect, which can be explained by the high stress zone shape under mode-II tearing condition. d=4mm d=5mm d=6mm
1.08
λ
1.04 1.00 0.96 0.92 0
30
60
90
θ(°)
120 150
Fig. 14. Mode-II SIF variation of a major edge crack with micro inclusions Fig. 15 provides the Mises stress contour plots for the cases of θ=30°, θ=60° and θ=150°. When θ=60°, the presence of the micro inclusions have minimal effect on the stress concentration of the major crack tip. When θ=150°, stress concentration appears around the circular inclusions, and the stress concentration of the major crack is relieved. When θ=30°, the stress concentration is enhanced.
a
b
c
Fig. 15. Von Mises stress plot for the subdomain of the case with (a) θ=30°, (b) θ=60° and (c) θ=150° (unit: MPa) 5.4 Case IV: edge crack in presence of symmetrical micro voids The effect of micro voids is quite like that of micro inclusions. The voids can be considered as inclusions with zero elastic modulus. The effects of the micro voids on the major crack tip have the same trend with those of micro inclusions, only with enhanced magnitude. Therefore, in this subsection is only provided the comparison of the influence factors. Both mode-I and mode-II problems are considered, the configuration is the same as that shown in Fig. 10 and Fig. 13. The deformed geometry is presented by stepwise zooming in Fig. 16 as an example.
b
a
c
Fig. 16. Deformed geometry of example case for micro voids problem In Fig. 17 is shown the comparison between the influence factors of micro inclusions and that of micro voids, with the same distance of d=4mm. As expected, the presences of micro voids have the
same kind of effect as that of micro inclusion, but with enhanced magnitude. The magnitude of the influence factor λ depends on elastic modulus of the inclusion and the distance d.
a inclusion void
1.1
λ
λ
1.1
b
1.0
1.0 0.9
inclusion void
0.9 0
30
60
90
θ(°)
120 150
0.8
0
30
60
90
θ(°)
120 150
Fig. 17. Comparison of influence factors under (a) mode-I configuration and (b) mode- II configuration In order to enhance the fatigue resistance of structures working with a main crack, symmetrical minor-voids of can be inserted in particular location to reduce the SIF of the main crack tip. Micro voids introduce smaller stress concentration in the vicinity of the defects than micro-cracks. According to the results present in this sub-section, symmetrical micro voids should be inserted at location as θ=60° for mode-I condition and θ=150° for mode- II condition. 6. Conclusions In the present work, the multiscale analysis of micro defects (minor cracks, inclusions, voids) in vicinity of a major crack tip is carried out by using multiscale XFEM method. The use of multiscale projection method makes it efficient to numerically simulate micro discontinuities with significantly small geometric dimension. The influences of the micro defects are evaluated by the SIF of major crack tip. Various configurations of micro defects are considered. On the basis of these simulations, the following conclusions can be drawn: (1) The present multiscale method has been found ideal for simulation of interaction between macrocrack and multiple micro defects in immediate proximity. (2) The SIF of the main crack tip is significantly affected due to the presence of micro defects, the magnitude of the effect the micro defects depend on their positions. (3) In general, those micro defects in front of the major crack tip have amplification effect, and behind the major crack tip have shielding effect. (4) For mode-I problem, the neutral angle that differentiate the amplification and shielding effect is found out to be θn≈70°~80°, which is in good agreement with that obtained by exact analytical method. (5) The relation between influence factor λ and the location of the micro defects are quite different for mode-Ⅰ problem and mode-Ⅱ problem. The failure of component is mainly decided on the SIF of main crack tip, which influenced by close micro defects. According to superposition principle, . The present study help to accurately predict the failure of crack propagation. Micro voids can be deliberately inserted in particular location
according to the results of present study, to reduce the SIF of the main crack tip, consequently enhance the fracture resistance of structures working with main crack. By the virtue of XFEM, the main crack propagation in presence of micro defects can also be easily simulated by the present method without any difficulty. The present method can be further extended to 3D cracking problem. Acknowledgements Support from National Natural Science Foundation of China (Nos. 51278297 and 51679139), the Major Program of the National Natural Science Foundation of China (No. 51490674), Research Program of Shanghai Leader Talent（No.20）and Doctoral Disciplinary Special Research Project of Chinese Ministry of Education (No.20130073110096) are acknowledged.
References: [1] A.P. Jivkov, Structure of micro-crack population and damage evolution in quasi-brittle media, THEOR APPL FRACT MEC. 70 (2014) 1-9. [2] V.P. Tamuzs, V.E. Petrova, On macrocrack-microdefect interaction, International Applied Mechanics. 38 (2002) 1157-1177. [3] M. Kachanov, J. Laures, Three-dimensional problems of strongly interacting arbitrarily located penny-shaped cracks, INT J FRACTURE. 41 (1989) 289-313. [4] S.X. Gong, On the formation of near-tip microcracking and associated toughening effects, ENG FRACT MECH. 50 (1995) 29-39. [5] S.A. Meguid, P.E. Gaultier, S.X. Gong, A comparison between analytical and finite element analysis of main crack-microcrack interaction, ENG FRACT MECH. 38 (1991) 451-465. [6] A. Brencich, A. Carpinteri, Interaction of a main crack with ordered distributions of microcracks: a numerical technique by displacement discontinuity boundary elements, INT J FRACTURE. 76 (1995) 373-389. [7] E. Svenning, M. Fagerstrom, F. Larsson, Computational homogenization of microfractured continua using weakly periodic boundary conditions, COMPUT METHOD APPL M. 299 (2016) 1-21. [8] G. Bhardwaj, I.V. Singh, Fatigue crack growth analysis of a homogeneous plate in the presence of multiple defects using extended isogeometric analysis, J BRAZ SOC MECH SCI. 37 (2015) 1065-1082. [9] S. Kumar, I.V. Singh, B.K. Mishra, A homogenized XFEM approach to simulate fatigue crack growth problems, COMPUT STRUCT. 150 (2015) 1-22. [10] J. Ortiz, A.P. Cisilino, J.L. Otegui, Effect of microcracking on the micromechanics of fatigue crack growth in austempered ductile iron, FATIGUE FRACT ENG M. 24 (2001) 591-605. [11] J. Mergheim, A variational multiscale method to model crack propagation at finite strains, INT J NUMER METH ENG. 80 (2009) 269-289. [12] X. Bian, Z. Li, G.E. Karniadakis, Multi-resolution flow simulations by smoothed particle hydrodynamics via domain decomposition, J COMPUT PHYS. 297 (2015) 132-155. [13] J. Song, Y. Yoon, Multiscale failure analysis with coarse-grained micro cracks and damage, THEOR APPL FRACT MEC. 72 (2014) 100-109.
[14] S. Loehnert, T. Belytschko, A multiscale projection method for macro/microcrack simulations, INT J NUMER METH ENG. 71 (2007) 1466-1482. [15] C. Ye, J. Shi, G.J. Cheng, An eXtended Finite Element Method (XFEM) study on the effect of reinforcing particles on the crack propagation behavior in a metal–matrix composite, INT J FATIGUE. 44 (2012) 151-156. [16] I.V. Singh, B.K. Mishra, S. Bhattacharya, R.U. Patil, The numerical simulation of fatigue crack growth using extended finite element method, INT J FATIGUE. 36 (2012) 109-119. [17] I.V. Singh, B.K. Mishra, S. Bhattacharya, XFEM simulation of cracks, holes and inclusions in functionally graded materials, INT J MECH MATER DES. 7 (2011) 199-218. [18] N. Sukumar, D.L. Chopp, N. Mo S, T. Belytschko, Modeling holes and inclusions by level sets in the extended finite-element method, COMPUT METHOD APPL M. 190 (2001) 6183-6200. [19] Q.B. Tinh, S. Hirose, C. Zhang, T. Rabczuk, C. Wu, T. Saitoh, J. Lei, Extended isogeometric analysis for dynamic fracture in multiphase piezoelectric/piezomagnetic composites, MECH MATER. 97 (2016) 135-163. [20] Q.B. Tinh, Extended isogeometric dynamic and static fracture analysis for cracks in piezoelectric materials using NURBS, COMPUT METHOD APPL M. 295 (2015) 470-509. [21] Q.B. Tinh, C. Zhang, Extended finite element simulation of stationary dynamic cracks in piezoelectric solids under impact loading, COMP MATER SCI. 62 (2012) 243-257. [22] Z. Kang, Q.B. Tinh, N. Du Dinh, T. Saitoh, S. Hirose, An extended consecutive-interpolation quadrilateral element (XCQ4) applied to linear elastic fracture mechanics, ACTA MECH. 226 (2015) 3991-4015. [23] Z. Kang, T.Q. Bui, T. Saitoh, S. Hirose, Quasi-static crack propagation simulation by an enhanced nodal gradient finite element with different enrichments, THEOR APPL FRACT MEC. (2017) 61-77. [24] Z. Wang, T. Yu, Q.B. Tinh, S. Tanaka, C. Zhang, S. Hirose, J.L. Curiel-Sosa, 3-D local mesh refinement XFEM with variable-node hexahedron elements for extraction of stress intensity factors of straight and curved planar cracks, COMPUT METHOD APPL M. 313 (2017) 375-405. [25] Z. Wang, T. Yu, Q.B. Tinh, A.T. Ngoc, T.H.L. Nguyen, D.D. Nguyen, D. Duc Hong, Numerical modeling of 3-D inclusions and voids by a novel adaptive XFEM, ADV ENG SOFTW. 102 (2016) 105-122. [26] G. Liu, D. Zhou, Y. Bao, J. Ma, Z. Han, Multiscale simulation of major crack/minor cracks interplay with the corrected XFEM, ARCH CIV MECH ENG. 17 (2017) 410-418. [27] G. Liu, D. Zhou, J. Ma, Z. Han, Numerical investigation of mixed-mode crack growth in ductile material using elastic–plastic XFEM, J BRAZ SOC MECH SCI. 38 (2016) 1689-1699. [28] T. Fries, A corrected XFEM approximation without problems in blending elements, INT J NUMER METH ENG. 75 (2008) 503-532. [29] N. Moës, M. Cloirec, P. Cartraud, J.F. Remacle, A computational approach to handle complex microstructure geometries, COMPUT METHOD APPL M. 192 (2003) 3163-3177. [30] M.J. Pais, F.A.C. Viana, N.H. Kim, Enabling high-order integration of fatigue crack growth with surrogate modeling, INT J FATIGUE. 43 (2012) 150-159. [31] G. Ventura, R. Gracie, T. Belytschko, Fast integration and weight function blending in the extended finite element method, INT J NUMER METH ENG. 77 (2009) 1-29. [32] S. Loehnert, T. Belytschko, Crack shielding and amplification due to multiple microcracks
interacting with a macrocrack, INT J FRACTURE. 145 (2007) 1-8. [33] S.X. Gong, Microcrack interaction with a finite main crack: an exact formulation, INT J FRACTURE. 66 (1994) R51-R56.
Highlights Multiscale numerical simulation of macro crack with micro defects problem; The simulation is conducted by using XFEM and multiscale projection method; The relation between the position of micro defects and their effects are revealed in detail.