- Email: [email protected]
Edge-based smoothed extended finite element method for dynamic fracture analysis
Accepted Manuscript
Edge-based smoothed extended finite element method for dynamic fracture analysis L. Wu , P. Liu , C. Shi , Z. Zhang , T. Quoc Bui , D. Jiao PII: DOI: Reference:
S0307-904X(16)30277-3 10.1016/j.apm.2016.05.027 APM 11180
To appear in:
Applied Mathematical Modelling
Received date: Revised date: Accepted date:
12 November 2015 9 May 2016 19 May 2016
Please cite this article as: L. Wu , P. Liu , C. Shi , Z. Zhang , T. Quoc Bui , D. Jiao , Edge-based smoothed extended finite element method for dynamic fracture analysis, Applied Mathematical Modelling (2016), doi: 10.1016/j.apm.2016.05.027
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights DSIFs was computed with the domain-form of the interaction integral .
CR IP T
ES-XFEM is not sensitive to mesh distortion on solving dynamic fracture problems.
When smoothing subdomains nssd ≥ 3 , relative error of DSIFs tend
AC
CE
PT
ED
M
AN US
to be stable.
1
ACCEPTED MANUSCRIPT
Edge-based smoothed extended finite element method for dynamic fracture analysis L. Wu 1*, P. Liu 1, C. Shi 1*, Z. Zhang 2，T. Quoc Bui 3, D. Jiao 4
College of Civil Engineering, Hunan University, Changsha, Hunan 410082, China
2
Centre for Future Materials, The University of Southern Queensland, Toowoomba QLD 4350, Australia 3
CR IP T
1
Department of Mechanical and Environmental Informatics, Tokyo Institute of Technology, O-okayama Meguro, Tokyo 152-8552, Japan
State Key Laboratory of Green Building Materials, China Building Materials
AN US
4
Academy, Beijing 100024, China
M
Abstract
In this paper, the edge-based smoothed extended finite element method (ES-XFEM) is applied and
ED
extended to the dynamic fracture analysis of two-dimensional elastic solids. The dynamic stress
PT
intensity factor (DSIF) is deduced from J-integral with the train smoothing method. With the deduced DSIF, the implementation of numerical integration and sub-smoothing regions are
CE
verified in dynamic fracture behaviors of linear elastic solids, which can be concrete and rocks.
AC
The analysis shows that ES-XFEM is an effective and efficient numerical approach to simulate the dynamic fracture problems. Keywords: ES-XFEM; Dynamic fracture mechanics; Dynamic stress intensity factor; Strain smoothing
2
ACCEPTED MANUSCRIPT
Nomenclature
nj
outward normal vector;
ij
Cauchy stress tensor;
bi
physical;
ui , ti u , t ;
density;
ij , ij
stress and strain vectors, respectively;
Cijkl
elastic matrix; reference point；
x*
a point on crack which nearest to x point；
en
outward normal vector of crack element for x* ;
Ni , N j , N k
all nodes set in the grid；
K , K
ED
respectively means Heaviside Functions and crack-tip enriched node sets; nodal displacement vector for parts;
a j and b k
PT
ui
shape function for the node i, j, k；
M
I
AN US
x
CR IP T
displacement and force respectively on boundary ;
strengthen variables for the node; respectively represents the system kinetic energy, strain energy,
CE
T , U , Wd and Ws
AC
damping energy and potential energy force. strain vector; ε ζ stress vector; C is the elastic matrix.
Fb , Fs and F
is defined as enriched force, surface force and concentrated force
respectively; K stiffness matrix; M mass matrix; R load vectors. i
means area of smoothing boundary domain; 3
ACCEPTED MANUSCRIPT
NI
I -shape function node;
n(i ) ( x)
outward normal vector matrix of the boundary; total number of the boundary segment ji ；
M
period of midpoint coordinates (gauss point) in the border ji ；
n jhi 、 l j i
respectively means the outward normal vector and length of ji ;
N e( i )
smoothing area;
I
all nodes of the k set in smoothing zone；
K , K
CR IP T
x GP j
enriched node set respectively by Heaviside function and crack-tip function in a ;
Biu 、B j and Bbk a
AN US
smoothing domain
denote boundary smoothed strain matrix in conventional smoothing domain,
Heaviside function and enriched smoothed strain by crack-tip function;
Aks
k th smoothing domain;
h
separated into x, y ;
M
N seg
boundary number of the smoothing domain;
N gau
ED
Gaussian integral points on the boundary;
,x 、 , y
,x 、 ,x 1
0
a contour around the crack-tip;
AC
Ans
perform transformation between the whole and local coordinates;
CE
Aeffm
2
PT
need to calculate the crack-tip under local coordinate system;
area of the mth element in the equivalent integral domain N effd
;
area of the the mth element in the nth smoothing edge;
K dyn and Kexact
respectively normalized numerical solution and analytical solution of the
dynamic stress intensity factor
First author: Linmei Wu (1987- ), female, postgraduate. Correspondent author: Linmei Wu (1987- ), female, postgraduate. E-mail addresses: [email protected]; [email protected] 1. Introduction 4
ACCEPTED MANUSCRIPT
The modeling of fracture mechanics is important to describe the singular fields near the crack-tip of a solid. However, numerical simulation of dynamic fracture problems remains a challenge in many applications, and is more complex compared to the numerical simulations under static loading conditions. The studies of dynamic fracture mechanics usually give
CR IP T
consideration to the responses of time-dependent loads with inertial effects. However, due to the mathematical complexity, the deriving analytical models/approaches are hard to be used in solving the fracture problems in practice regarding general dynamic element analysis. Numerical methods
AN US
may solve the limitation of analytical model. A number of numerical methods, such as extended finite element method (XFEM) [1] and Meshfree method [2-4], have been applied to solve fracture problems in the last decades. The XFEM has become the most widely used approach to study the arbitrary crack growth. However, the computation accuracy of XFEM is dependent on
M
the mesh regulation as it is based on the level set method (LSM) [5] and partition of unity method
ED
(PUM) [6], which are similar to classical FEM.
PT
To improve the computation accuracy, , Liu et al. [8] proposed a cell-based finite element method (CS-FEM) in 2007, which is based on strain smoothing technique [7]. . Following this,
CE
smoothed finite element method (SFEM) was developed and up-grated to node-based smoothed
AC
finite element method (NS-FEM) [9], edge-based smoothed finite element method FEM (ES-FEM) [10], superconvergent alpha-based finite element method (SαFEM) [11] and face-based smoothed finite element method FEM (FS-FEM) [12] by the same team. Because the integral calculation of the stiffness matrix by SFEM is only at smoothed region boundary, therefore there is a significant reduction of the requirement of
the quality of mesh assemble but it still keeps high accuracy. In
brief, the advantages of SFEM include (1) singular integration will not occur in the calculation of 5
ACCEPTED MANUSCRIPT
crack-tip because no derivation of the shape function is used; (2) integration calculation is focused on smoothed region boundary, instead of a throughout mapping; and (3) insensitive to mesh distortion, which is suitable for the analysis of fracture problems under large deformation. Because of these advantages, smoothed XFEMs are extensively studied recently, particularly on
CR IP T
the fracture problems. Bordas et al. [13] combined CS-FEM with XFEM to contrive cell-based finite element method (CS-XFEM) and applied it to the analysis of the fracture problem of thin plate. The results by the authors showed that the stress intensity factor was affected by the mesh
AN US
distortion and the convergence was better than using conventional XFEM. Chen et al. [14] deduced the edge-based smoothed extended finite element method (ES-XFEM) and applied it to linear elastic crack propagation, and they suggested that the precision and convergence of
M
ES-XFEM were better than CS-XFEM. Jiang et al. [15] applied ES-XFEM to orthotropic material fracture mechanism, which proved again that
ES-XFEM was more accurate than XFEM.
ED
However, the effectiveness and efficiency of applying ES-XFEM to the analysis of dynamic
PT
fracture problems are rarely reported, remaining unknown. The aim of this study to examine the effectiveness and efficiency of ES-XFEM for the
CE
analysis of dynamic fracture problems in 2D elastic solids. To validate the accuracy of the method,
AC
numerical results from DSIFs are compared with the analytical results and other numerical solutions.
2. Solution of dynamic fracture 2.1 Problem description
6
ACCEPTED MANUSCRIPT
The diagram of a plane body with a crack is shown in Fig. 1. The boundary of is composed of the displacement boundary u , stress boundary t and the crack c . The whole domain is recorded as .
c
y u Γu
t
AN US
x
CR IP T
t
Fig.1. Scheme of crack problem
The motion and boundary conditions are given as follows:
ij , j bi ui
(1)
on t
(2)
ij n j 0
on c
(3)
ui ui
on u
(4)
PT
ED
ij n j ti
M
in
CE
where n j is outward normal vector; ij is the Cauchy stress tensor; bi is the body force; ui and ti
t , respectively; is the
AC
are the displacement and stress vectors on boundary u and density, is the domain of the body.
Considering small deformation, the geometric equation is expressed as follows:
ij
1 ui, j u j ,i 2
The constitutive equation is:
7
(5)
ACCEPTED MANUSCRIPT
ij Cijkl kl
(6)
where ij , ij are the stress and strain vectors, respectively; Cijkl are components of the elastic tensor. 2.2 Displacement mode
CR IP T
Compared with classical finite element method, the displacement mode of XFEM takes more functions to describe the models of the discontinuous region. Displacement discontinuity equation of crack elements is expressed demonstrated with enriched functions, such as the
AN US
Heaviside function H x [1]. Above the crack, H x is positively 1; while below the crack it is negative 1:
1 x x* en 0 H x * 1 x x en 0 *
represents a point on crack which is next to
M
In Equation (7), x means reference point; x
(7)
ED
x ; en is the outward normal vector of crack element of x* . In order to capture the singularity and discontinuity near crack-tip, the crack-tip enriched
PT
functions are used. For an isotropic elastic body, the crack-tip enriched functions usually use the
CE
branch crack-tip x [1], which is expressed as follows:
r, 4
r sin 2
r cos 2
r sin sin 2
r cos sin 2
AC
1
where r , is the crack-tip local coordinates in equation (8), and the first term
(8)
r sin 2
describes the discontinuity of crack surface, while the other three terms are employed to improve the singularity of crack-tip [16]. In this way, XFEM displacement mode can be expressed in the following form: 8
ACCEPTED MANUSCRIPT
uh x Ni x ui iI
4
N x H x a N x x b j
j
jK
k
(9)
k
1
k K
where Ni , N j , N k Ni, Nj and Nk are the shape functions for nodes i, j andk, respectively; I is the total node set in the mesh； K and K are Heaviside Functions and crack-tip enriched node sets, are two
CR IP T
respectively; ui are nodal displacement vectors of different parts; a j and b k strengthen variables of nodes.
According to the above XFEM displacement mode, it is not convenient to impose Dirichlet
enriched functions are proposed [17]:
uh x Ni x u i i I
AN US
boundary in the calculation of an unknown quantity . To solve this problem,
N x H x a j
j
K
j
k K
the shifted-basis
4
N kx x b 1
k
(10)
M
where, H x H x H x j ; x x xk . 2.3 The description of dynamic fracture
ED
According to Hamilton’s law of variational principle [17], the governing equation of the
PT
XFEM dynamics can be derived as follows Firstly, Lagrange functional system is expressed as:
L T U Wd Ws
(11)
CE
where T , U , Wd and Ws represents the system kinetic energy, strain energy, damping energy
AC
and potential energy of external force, respectively. The kinetic energy T is: T
1 δ δd 2
T
where δ and δ are displacement and velocity vector, respectively. Strain energy U can be calculated by: 9
(12)
ACCEPTED MANUSCRIPT
U
1 T 1 T ε ζd ε Cεd 2 2
(13)
where ε is strain vector, ζ is the stress vector and C is the elastic matrix. The damping potential energy Wd is:
1 c δ 2
Wd
T
δd
(14)
Potential energy of external force Ws is expressed by:
CR IP T
where c is the viscous damping coefficient.
Ws δ Fb d δ Fs d δ F T
T
T
(15)
AN US
where Fb , Fs and F are defined as enriched force, surface force and concentrated force vector, respectively.
According to the Hamilton's principle, the conservative system is expressed by: t2
Ldt 0
(16)
M
t1
Neglecting of damping, and combining the equations (11) to (15) with the XFEM
Mδ Kδ R
(17)
PT
where δ u a b
ED
displacement mode substitution equation (16), the equation of dynamic XFEM is expressed as:
T
;
K is stiffness matrix; M is the mass matrix and R is load vector.
CE
2.4 Solution of the dynamic equation dynamic equation can be classified as (1) implicit integration
AC
The methods to solve
(Newmark method, Wilson-method, etc.) and (2) explicit integration (center difference method, etc.)[18]. In this works Newmark method was adopted to resolve equation (17). At the n th time step, discrete simulation equations can be expressed as:
M t K δ 2
n
t 2 F K δn 1 tδn 1 1 2 δn 1 2 10
(18)
ACCEPTED MANUSCRIPT
δn δn1 1 δn1 δn t
(19)
1 δn δn 1 δn 1t δn 1 δn t 2 2
(20)
2 where t is the time increment; when 0.5 and 0.25( +0.5) , the Newmark method is
conditions, the effect of damping is therefore not considered
CR IP T
unconditionally stable. This study is to analyze the dynamic fracture of crack under the shock load
3. Edge-based smoothed extended finite element method (ES-XFEM) for crack problems
AN US
3.1 Smoothing domain of ES-FEM
Comparing with conventional finite element, the integral stiffness matrix of SFEM is based on the smoothing region, rather than anelement. Smoothing domain of ES-XFEM is based on edge
1
by
2
M
of an element (as shown in Fig.2). The entire smoothing domain is expressed Ns
, and i
j , i j , where N s represents the
number N s edge. The ABCD ( ) and EFG ( ) are the smoothing domains based on edge m
ED
k
BD ( ) and EF ( ), respectively. k
AC
CE
PT
k
Fig.2. Scheme of edge-based smoothing finite element method 11
ACCEPTED MANUSCRIPT
According to the smoothing operation [11], the boundary strain
i of the element i can be
obtained from the compatible strain ( x) in the region of the smoothing domain i :
i
( x) ( x)d u( x) ( x)d i
s
(21)
i
i
i
Smoothing function i ( x ) meets:
i ( x ) d 1
The smoothing function can be expressed as follows:
1 i Ai 0
x i
(23)
AN US
x
i
(22)
CR IP T
i
where i is area of smoothing boundary domain [2]. i
For smoothing area , the stiffness matrix is: i
i
M
k
B CB T I
J
d
(24)
i
ED
Under the bifurcation theory [2], the smoothed strain matrix B I of the node
I
in domain
CE
PT
i can be expressed as:
b Ix BI 0 b Iy
0 b Iy b Ix
(25)
AC
where
b Ih
1 A( i )
N I ( x )n i ( x )d
( h x, y )
(26)
i
(i ) and N I is I -shape function node; n ( x) is outward normal vector matrix of the boundary
i , which is represented as follows:
12
ACCEPTED MANUSCRIPT nx( i ) n(i ) ( x) 0 n(yi )
0 n(yi ) nx( i )
(27)
As linear displacement field is used at in , a Gaussian point is applied at . i
i
i
Equation (26) can be further expressed in the algebra form:
1
M
N A i
j 1
I
( x GP )n(jhi ) l ji j
( h x, y )
(28)
CR IP T
b Ih
i GP where M is the total number of the boundary segment j ; x j means the period of midpoint
i i i coordinates (gauss point) in the border j ; n jh and l j are the outward normal vector and
AN US
i length of j , respectively. Smoothing region of the stiffness matrix can be further expressed as
follows:
k
i
B I C B J Ai T
(29)
M
From equation (26), the stiffness matrix of the shape function does not require partial
ED
derivatives. The whole stiffness matrix K assembling each region stiffness matrix has beenset,
PT
which is:
N e( i )
K = k
i
(30)
i 1
CE
where N e( i ) is the smoothing area.
AC
3.2 Equations of ES-XFEM ES-XFEM is Deriving from equation (17), the assembled mass matrix and load vector for
ES-XFEM do not involve strain, therefore dynamic extension of finite element method (FEM) is needed. Similar to XFEM, the smoothed strain of the displacement mode of ES-XFEM in the smoothing domain can be written by: 13
ACCEPTED MANUSCRIPT
k Biu ui iI
jK
B aj a j
4
B b b
kK
1
k
k
(31)
in the smoothing zone； K and K mean enriched
where I means all nodes of the set k
node set by Heaviside function and crack-tip function in a smoothing domain, respectively; Biu ,
CR IP T
B aj and Bbk denote boundary smoothed strain matrix in conventional smoothing domain,
Heaviside function and enriched smoothed strain by crack-tip function, respectively [14]. To enrich the matrix, its specific expressions are as follows:
0 biyr bixr
where
ks
nh x Ni x d
1 Aks
ks
nh x Ni x H x d
ED
biha
1 Aks
1 Aks
ks
nh x Ni x x d
(33a)
(33b)
1 , 2 , 3 , 4
(33c)
PT
bihb
(32)
M
bihu
r u , a, b
AN US
bixr Bir 0 biyr
where Aks is the number k smoothing domain, h is separated into x, y .
AC
CE
Equation (33) can be presented in numerical forms:
bihu
biha
1 Aks
1 b s Ak b ih
1 Aks
N gau wm,n Ni xm n, nh xm n, m 1 n 1 N seg
(34a)
N gau wm,m Ni xm,n H xm,m nh xm ,n m 1 n 1
(34b)
N gau wm,m Ni xm,n xm,m nh xm,n m 1 n 1
(34c)
N seg
N seg
14
ACCEPTED MANUSCRIPT
where N seg is the boundary number of the smoothing domain; N gau represents Gaussian integral
M
AN US
CR IP T
point on the boundary.
Fig.3. Schematic diagram of edge-based smoothed extended finite element method
ED
The stiffness matrix kij is derived from the smoothed strain matrix:
CE
PT
B u T CB u d j ks i T kij s Bia CB uj d k B b T CB u d j ks i
B
B
CB bj d T T a a a u ks Bi CB j d ks Bi CB j d b T a u T u ks Bi CB j d ks Bi CB j d ks
u T i
CB aj d
ks
u T i
(35)
AC
Due to smoothed strain matrix in smoothing domain is constant, the equation (35) can be
written as
B u T CB u As j k i a T kij Bi CB uj Aks b T u s Bi CB j Ak
B B B
u T i
CB aj Aks
a T i
CB aj Aks
b T i
CB aj Aks
15
B B B
CB bj Aks T b s CB j Ak T CB bj Aks
u T i
a i b i
(36)
ACCEPTED MANUSCRIPT
e Using m to express the element mass matrix set:
mijuu mije mijau mijbu
mijub mijau mijbb
mijua mijaa mijba
(37)
where
mijuu e Ni N j d T
(38a)
mijua e Ni
T
mijub e Ni
T
HN d
mijaa e HNi
T
T
N j d,
(38c)
(38d)
i
AN US
mijbb e Ni
1,2,3,4
HN d
mijab e HNi N j d,
(38b)
j
N j d,
T
CR IP T
1,2,3,4
1,2,3,4
(38e) (38f)
Using ri e to express the element nodal load vector components:
rib1
rib 2
rib3
rib 4
(39)
ED
where
ri a
M
rie riu
(40)
CE
PT
riu e Ni Fb d Ni Fs d Ni F a ri e N i HFb d N i HFs d N i HF b ri e Ni Fb d Ni Fs d Ni F 1, 2,3, 4
AC
4. Integration implementation of ES-XFEM Because the stiffness matrix kij is derived from integration on the smoothing boundary, the
conventional XFEM integral solution is no longer applicable. A new integral formula has been proposed in [9], which is expressed as follows: (1) As shown in Fig.3, each boundary employs a Gaussian point in the conventional smoothing domain (such as ABCD). 16
ACCEPTED MANUSCRIPT
(2) Each boundary employs a Gaussian point in crack-free enhanced Heaviside function domain. (3) Each boundary adopts 3~5 Gaussian points in crack-free and crack-tip enriched smoothing domain.
CR IP T
(4) In the smoothing domain where the crack run through, it can be divided into two smoothing subdomain and one Gaussian point is used on the each smoothing boundary; As shown in Fig.4, irregular smoothing subdomain is usually divided into a triangle sub-domain in order to
PT
ED
M
AN US
facilitate programming.
CE
Fig.4. Partition of split smoothing domains in the ES-XFEM[14]
AC
(5) As shown in Fig.5, crack is extended to smoothing boundary of the region and then divide the crack-tip domain into multiple sub-domains, generally three to five ones. The boundary of the smoothing subdomain near the crack-tip need 5 Gaussian points, and other boundary needs 1 ~ 3 points.
17
ACCEPTED MANUSCRIPT
(6) As shown in Fig.6, in the integral of enriched function, the partition rule of smoothing sub-domain is same as above, but Gaussian points are not taken on the smoothing sub-domain boundary and while adopt the interior points of the smoothing sub-domain. This is because if the smoothing domain on the boundary points is taken as a Gaussian point, boundary AB, CD, EF,
CR IP T
GH are coincidence with crack surface and then cannot satisfy the accuracy of calculation. In general, the crack-tip items of enriched function utilize about 3~5 Gaussian points, otherwise enhanced Heaviside takes one Gaussian point. Sub-domain integration method is employed to
AN US
calculate both Mass matrix and stiffness matrix in conventional XFEM. It is needed to employ one Gaussian point in the internal unit of unenhanced nodes, 5 ~ 7 points in crack-tip sub-domains and
AC
CE
PT
ED
M
3 points in other domains.
18
ACCEPTED MANUSCRIPT
CR IP T
Fig.5. Partition of crack-tip smoothing domains in the ES-XFEM
AN US
Fig.6. Scheme of Gaussian point in split and crack-tip element
5. Dynamic stress intensity factor
M
Stress intensity factor can be calculated by displacement extrapolation method [19], virtual crack extension method[20], virtual crack closure legal [21], and interaction integral method[22].
ED
In this work, interaction integral method has been used. Considering the two dimensional crack,
PT
crack is described as a line, while a local coordinate is set up at crack-tip. As shown in Fig.7, 0 means a contour around the crack-tip.
CE
Contour integrals J near the crack-tip can be expressed as:
AC
J lim
0 0
0
W
1i
ij u j ,1 ni d lim
0 0
1 2
0
W
where W ij ij , which is the strain energy density [18].
19
1i
ij u j ,1 mi d
(43)
CR IP T
ACCEPTED MANUSCRIPT
AN US
Fig.7. Scheme of J-integral domain
Equation (43) cannot be directly applied to numerical calculation and then adopted by a 0 1 closed domain s 1 0 . The set weight factor q is 1 in , and is 0 in ,
M
Equation (43) can be rewritten as:
J ij u j ,1 W 1i qmi d s
2 j u j ,1d
ED
(44)
Using the divergence theorem, the integral form of the J integral is
PT
J ij u j , 1 W i 1 q i d,A ij u j i , 1W qd A , 1 A
A
(45)
CE
By incorporating the balance equations and geometric equation: J ij u j ,1 W 1i q,i dA uu j ,1 ui u,1 qdA A
A
(46)
Choose two independent equilibrium states: state 1 ( ij , ij , u i(1) ) as the real state, state 2
AC
(1)
(1)
( ij , ij , u i( 2 ) ) as the auxiliary condition. The third equilibrium is achieved by linear ( 2)
( 2)
1,2 superposition of the two equilibrium, with J integral can be rewriten by J :
J 1 , 2 J 1 J 2 I
1, 2
where I 1,2 means the interaction integral of state 1 and state 2: 20
(47)
ACCEPTED MANUSCRIPT
I 1,2 1u j 2,1 2u j1,1 W 1,21i q,i dA u 1u j2,1 qdA A
(48)
A
where W 1,2 ij 2 ij1 ij1 ij 2 is the interaction of strain energy density. According to [1], the two-dimensional equation of energy release rate can be expressed as [1]:
E with E E 1 v 2
K I2 K II2 E E
(49)
plane stress plane strain
.
CR IP T
J
For the third equilibrium, equation (49) can be reformed as follows:
J
1
J
2
2 K I1 K II 2 K I 2 K II1
AN US
J
1,2
E
(50)
Comparing with Equation (48) and Equation (50), it is expressed as:
2 K I1 K II 2 K I 2 K II1
E
(51)
M
I
1,2
ED
State 2 is chosen as asymptotic field of mode I and II, then:
E 1,mod eI I 2
K II1
E 1,mod eII I 2
(52)
PT
K I1
CE
As shown in Fig.7 (b), the interaction integrals are based on the smoothing boundary,no longer
AC
based on the element integrals for ES-XFEM, so equation (49) is rewritten as:
I
1,2
d Neff
m 1
m Aeff
ij
d Neff
u j ,1 ij u j ,1 ij ij 1i q,i dA
1 2
2 1
1 2
m 1
m Aeff
u 1u j2,1 qdA
(53)
m d where Aeff is area of the number m element in the equivalent integral domain N eff .
In equation (53), the first part of the right side is related to the displacement derivative which is based on integral of the smoothing domain; the second item of right side is independent of the
21
ACCEPTED MANUSCRIPT
displacement derivative term, which is based on the element integral.
Equation (53) is depicted
as:
I
1,2
d Neff
3
m 1 n 1
Ans
ij
d Neff
u j ,1 ij u j ,1 ij ij 1i q,i dA
1 2
2 1
1 2
m 1
m Aeff
u 1u j2,1 qdA
(54)
s
CR IP T
where An is area of the the m-th element in the n-th smoothing edge. 6. Numerical examples
6.1 Example 1: dynamic fracture of edge crack in the semi-infinite domain
Because of computational limits, a limited domain is modeled for analysis, as shown in Fig.8.
AN US
The size of the limited plate are L = 10 m, H = 2 m, crack length a = 5 m；The elastic constants of materials are: E = 210 GPa, v = 0.3, ρ =8000 kg/m3. σ0 = 500MPa. The top of the plate is subject to a distributed step load σ0 (t) = σ0 H(t). When t>0, H(t)=1. Computation time is 1×10-3 s, time step
M
is 2×10-5 s and the number of iterations is 50 times. The regular meshes use 65×25 triangular
ED
elements, as shown in Fig. 9 (a). In order to conform rough requirement of mesh in the smoothed expansion finite element method, different irregular meshes are adopted. The irregular mesh is
CE
PT
generated randomly by [8]:
x ' x x rd ir ' y y y rd ir
(55)
AC
Where x , y respectively represents node coordinates of initial rules mesh, as shown in Fig.9.
(a) ; Accordingly, x , y is the length of element respectively in the x, y direction； rd is random numbers in the domain[-1, 1]; ir is irregular factors. Regular mesh ( ir 0.00 )and irregular mesh( ir 0.12 ) are shown is shown in Fig. 9.
22
ACCEPTED MANUSCRIPT
CR IP T
Fig.8. Geometry of a semi-infinite crack
Since reflection occurs when the stress wave meets the crack-tip, calculation process does not consider the impact of reflection, and the dynamic stress intensity factor analytic formula is valid
H/cd , cd E(1 v ) / ( (1 v ) (1 2v ) ) ).
AN US
just as time t 3tc ( tc
Regardless of the crack propagation, dynamic stress intensity factor analytic is expressed as [23]:
2 0 cd t (1 2v) 1 v
M
K Idyn (t )
(56)
ED
The relative error of normalized dynamic stress intensity factor is defined as follows:
PT
err
K dyn Kexact K exact
100%
(57)
CE
where K dyn and K exact are normalized numerical solution and analytical solution of the dynamic stress intensity factor, respectively.
AC
In the following simulation figures, K Idyn / 0 H1/2 , t / tc are the dynamic stress intensity
factor of the normalized mode I and the normalized time which the dimension is 1, respectively.
23
AN US
Fig.9. Computational meshes
CR IP T
ACCEPTED MANUSCRIPT
6.1.1 Convergence analysis
It is researched on convergence of different computational meshes. As can be seen from
M
Fig.10., when t 2.0tc , relative error of dynamic stress intensity factor by 43×19 meshes
ED
decreased better than 34 × 14 meshes. 53 × 21 meshes and 65 × 25 meshes have similar
AC
CE
PT
performance with convergence. The following analysis is based on 65×25 meshes.
Fig.10. (a) normalized dynamic stress intensity factor; (b) relative error of dynamic stress intensity factor at various meshes 6.1.2 Selection of crack-tip smoothing subdomain 24
ACCEPTED MANUSCRIPT
As the crack-tip region is divided into different number of smoothing subdomain and the influence of subdomains number on the dynamic stress intensity factor is analyzed. As can be seen from Fig.11.(a), when number of smoothing subdomain n ssd 2 , the difference of dynamic stress intensity factor in normalization mode I is not obvious. However, as shown in Fig.11.(b),
CR IP T
when n ssd 3 , relative error of the normalized dynamic stress intensity factor in mode I tends to be stable. Therefore, in order to improve the computation efficiency and guarantee the calculation accuracy, in the crack-tip smoothing domain sequence number should be used as n ssd =3 and
AN US
employs five Gaussian points at each boundary. The smoothing domain with crack across do not
PT
ED
M
need partition with one Gaussian point is taken; other area also take one Gaussian point.
CE
Fig.11. (a) normalized dynamic stress intensity factor; (b) relative error of dynamic stress intensity factor with different smoothing subdomain
AC
6.1.3 Accuracy and mesh distortion As shown in Fig.12, it is compared the ES-XFEM results with the analytical solution, finite
element method of singular boundary smoothing [24] (sES-FEM) and extended finite element method (XFEM). It is demonstrated that when t 1.5tc , numerical error of three methods are all relative large; when t 1.5tc , the results of three kinds of numerical methods are very close to
25
ACCEPTED MANUSCRIPT
the analytical solution; when t tc , the error of the numerical solution is the largest and the maximum error is XFEM. With classical finite element method, mesh distortion may affect the computational accuracy. ES-XFEM is no longer restricted, while Fig.13. shows the calculation results of different mesh
CR IP T
irregular factors. It is illustrated that the ES-XFEM results affected by the quality of the mesh are very small. When the irregular factor takes 0.12, the mesh is extremely irregular. There are 74 angles greater than 120o , 372 angles less than 30o and 39 angles less than 20o (as shown in
AN US
Fig.9 (b)). In this case, the ES-XFEM can still get a stable result as the time 1.8tc t 3.0tc , the maximum relative error of the dynamic stress intensity factor in the normalized mode I is less than
PT
ED
M
3%.
intensity factor comparing with different method
AC
CE
Fig.12. (a) normalized dynamic stress intensity factor (b) relative error of dynamic stress
26
CR IP T
ACCEPTED MANUSCRIPT
Fig.13. (a) normalized dynamic stress intensity factor (b) relative error of dynamic stress intensity factor with different irregular factor
AN US
6.2 Example 2: dynamic fracture of the bevel edge crack plate
In order to prove the efficiency and stability of ES-XFEM for dynamic fracture problems, we further investigative the mixed-mode dynamic fracture. An inline edge crack of plate is shown in Fig.14, its geometry sizes are L = 44 mm, H = 32mm, D = 16 mm, crack length is a = 22.63 mm,
M
and the incline angle of crack α = 45o. Material parameters are set as: modulus of elasticity E =
ED
29.4 GPa, Poisson’s ratio v = 0.286, the density ρ=2450 kg/m3. The left, top and bottom of plate
PT
bear restriction of chain support in normal direction, the right end bears a step load distributionσ0 (t) = σ0 H(t). Computation time is 3×10-5 s, Time step is 6×10-7s and the number of iterations is
CE
50 times. Two kinds of computational meshes (66×48) with different quality (regular and irregular mesh) are revealed in Fig.15.
AC
Calculation results in the figure below shows that, K Idyn / 0 πa
1/2
, K / dyn II
πa
1/2
0
are the normalized dynamic stress intensity factor of mode I and II, respectively. In Fig.16, analyzes and compares the normalization of the dynamic stress intensity factor under different numerical methods with the regular mesh and the results of three methods are
27
ACCEPTED MANUSCRIPT
consistent and the peak value of dynamic stress intensity factor in normalized mode I is higher than that in normalized mode II. Fig.14 presents ES-XFEM results under different mesh distortion. When the irregular factor takes 0.2, maximum of the element angle is 154.3°, minimum is 10.6°. There are 172 elements
CR IP T
with maximum angles greater than 120°, 485 elements with the minimum angles less than 30°and 30 elements with the minimum angles less than 20°. As shown in Fig.17, mesh distortion effected on the dynamic stress intensity factor is small. It is demonstrated again that ES-XFEM is not
AN US
sensitive to the mesh quality, which will help to analyze large deformation of crack in propagation
AC
CE
PT
ED
M
process.
28
ACCEPTED MANUSCRIPT
PT
ED
M
AN US
Fig.15 Computational meshes
CR IP T
Fig.14 Geometry of a rectangular plate with an inclined crack
edge crack
AC
CE
Fig.16 Comparison of the normalized mixed-mode DSIFs for a rectangular plate with a slanted
29
CR IP T
ACCEPTED MANUSCRIPT
Fig.17 The normalized mixed-mode DSIFs for a rectangular plate with a slanted edge crack
AN US
under different mesh irregular factor
7 Conclusions
This paper applies the edge based smoothed extended finite element (ES-XFEM) to solve
M
stationary dynamic fracture problems subjected to impact loads. Dynamic stress intensity factor
ED
(DIF) is computed with the domain-form of the interaction integral in terms of the strain
PT
smoothing method. Two case studies including the mode I and mixed crack prove that the dynamic ES-XFEM is highly precise, and insensitive to mesh distortion. It is also proved that
3 , the relative error of dynamic stress intensity factor become
CE
when smoothing subdomain n ssd
AC
stable.
References [1] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. International journal for numerical methods in engineering. 1999;45:601-20. [2] Amiri F, Anitescu C, Arroyo M, Bordas S, Rabczuk T. XLME interpolants, a seamless bridge 30
ACCEPTED MANUSCRIPT
between XFEM and enriched meshless methods. Computational Mechanics. 2014;53:45-57. [3] Amiri F, Millán D, Shen Y, Rabczuk T, Arroyo M. Phase-field modeling of fracture in linear thin shells. Theoretical and Applied Fracture Mechanics. 2014;69:102-9. [4] Nguyen-Xuan H, Liu G, Bordas S, Natarajan S, Rabczuk T. An adaptive singular ES-FEM for
CR IP T
mechanics problems with singular field of arbitrary order. Computer Methods in Applied Mechanics and Engineering. 2013;253:252-73.
[5] Osher S, Sethian JA. Fronts propagating with curvature-dependent speed: algorithms based on
AN US
Hamilton-Jacobi formulations. Journal of computational physics. 1988;79:12-49.
[6] Duarte CA, Oden JT. An h-p adaptive method using clouds. Computer methods in applied mechanics and engineering. 1996;139:237-62.
M
[7] Chen J-S, Wu C-T, Yoon S, You Y. A stabilized conforming nodal integration for Galerkin mesh-free methods. International journal for numerical methods in engineering. 2001;50:435-66.
ED
[8] Liu G, Dai K, Nguyen T. A smoothed finite element method for mechanics problems.
PT
Computational Mechanics. 2007;39:859-77. [9] Liu G, Nguyen-Thoi T, Nguyen-Xuan H, Lam K. A node-based smoothed finite element
CE
method (NS-FEM) for upper bound solutions to solid mechanics problems. Computers &
AC
structures. 2009;87:14-26. [10] Liu G, Nguyen-Thoi T, Lam K. An edge-based smoothed finite element method (ES-FEM) for static, free and forced vibration analyses of solids. Journal of Sound and Vibration. 2009;320:1100-30. [11] Liu G, Nguyen-Xuan H, Nguyen-Thoi T, Xu X. A novel Galerkin-like weakform and a superconvergent alpha finite element method (SαFEM) for mechanics problems using triangular 31
ACCEPTED MANUSCRIPT
meshes. Journal of Computational Physics. 2009;228:4055-87. [12] Nguyen-Thoi T, Liu G, Lam K, Zhang G. A face-based smoothed finite element method (FS-FEM) for 3D linear and geometrically non-linear solid mechanics problems using 4-node tetrahedral elements. International Journal for Numerical Methods in Engineering. 2009;78:324.
in FEM and XFEM. Computers & structures. 2010;88:1419-43.
CR IP T
[13] Bordas SP, Rabczuk T, Hung N-X, Nguyen VP, Natarajan S, Bog T, et al. Strain smoothing
[14] Chen L, Rabczuk T, Bordas SPA, Liu G, Zeng K, Kerfriden P. Extended finite element
AN US
method with edge-based strain smoothing (ESm-XFEM) for linear elastic crack growth. Computer Methods in Applied Mechanics and Engineering. 2012;209:250-65.
[15] Jiang Y, Tay T, Chen L, Sun X. An edge-based smoothed XFEM for fracture in composite
M
materials. International Journal of Fracture. 2013;179:179-99.
[16] Ahmed A, Auricchio F. Extended finite element method (XFEM)-modeling arbitrary
ED
discontinuities and failure analysis. Research degree thesis. 2009.
PT
[17] Stąpór P. Application of xfem with shifted-basis approximation to computation of stress intensity factors. Archive of Mechanical Engineering. 2011;58:447-83.
CE
[18] Hughes TJ. The finite element method: linear static and dynamic finite element analysis:
AC
Courier Corporation; 2012. [19] Shih C, Lorenzi Hd, German M. Crack extension modeling with singular quadratic isoparametric elements. International Journal of Fracture. 1976;12:647-51. [20] Parks DM. A stiffness derivative finite element technique for determination of crack tip stress intensity factors. International Journal of fracture. 1974;10:487-502. [21] Rybicki EF, Kanninen M. A finite element calculation of stress intensity factors by a 32
ACCEPTED MANUSCRIPT
modified crack closure integral. Engineering Fracture Mechanics. 1977;9:931-8. [22] Yau J, Wang S, Corten H. A mixed-mode crack analysis of isotropic solids using conservation laws of elasticity. Journal of Applied Mechanics. 1980;47:335-41. [23] Menouillard T, Belytschko T. Dynamic fracture with meshfree enriched XFEM. Acta
CR IP T
mechanica. 2010;213:53-69. [24] Liu P, Bui T, Zhang C, Yu T, Liu G, Golub M. The singular edge-based smoothed finite element method for stationary dynamic crack problems in 2D elastic solids. Computer Methods in
AC
CE
PT
ED
M
AN US
Applied Mechanics and Engineering. 2012;233:68-80.
33
Edge-based smoothed extended finite element method for dynamic fracture analysis L. Wu , P. Liu , C. Shi , Z. Zhang , T. Quoc Bui , D. Jiao PII: DOI: Reference:
S0307-904X(16)30277-3 10.1016/j.apm.2016.05.027 APM 11180
To appear in:
Applied Mathematical Modelling
Received date: Revised date: Accepted date:
12 November 2015 9 May 2016 19 May 2016
Please cite this article as: L. Wu , P. Liu , C. Shi , Z. Zhang , T. Quoc Bui , D. Jiao , Edge-based smoothed extended finite element method for dynamic fracture analysis, Applied Mathematical Modelling (2016), doi: 10.1016/j.apm.2016.05.027
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Highlights DSIFs was computed with the domain-form of the interaction integral .
CR IP T
ES-XFEM is not sensitive to mesh distortion on solving dynamic fracture problems.
When smoothing subdomains nssd ≥ 3 , relative error of DSIFs tend
AC
CE
PT
ED
M
AN US
to be stable.
1
ACCEPTED MANUSCRIPT
Edge-based smoothed extended finite element method for dynamic fracture analysis L. Wu 1*, P. Liu 1, C. Shi 1*, Z. Zhang 2，T. Quoc Bui 3, D. Jiao 4
College of Civil Engineering, Hunan University, Changsha, Hunan 410082, China
2
Centre for Future Materials, The University of Southern Queensland, Toowoomba QLD 4350, Australia 3
CR IP T
1
Department of Mechanical and Environmental Informatics, Tokyo Institute of Technology, O-okayama Meguro, Tokyo 152-8552, Japan
State Key Laboratory of Green Building Materials, China Building Materials
AN US
4
Academy, Beijing 100024, China
M
Abstract
In this paper, the edge-based smoothed extended finite element method (ES-XFEM) is applied and
ED
extended to the dynamic fracture analysis of two-dimensional elastic solids. The dynamic stress
PT
intensity factor (DSIF) is deduced from J-integral with the train smoothing method. With the deduced DSIF, the implementation of numerical integration and sub-smoothing regions are
CE
verified in dynamic fracture behaviors of linear elastic solids, which can be concrete and rocks.
AC
The analysis shows that ES-XFEM is an effective and efficient numerical approach to simulate the dynamic fracture problems. Keywords: ES-XFEM; Dynamic fracture mechanics; Dynamic stress intensity factor; Strain smoothing
2
ACCEPTED MANUSCRIPT
Nomenclature
nj
outward normal vector;
ij
Cauchy stress tensor;
bi
physical;
ui , ti u , t ;
density;
ij , ij
stress and strain vectors, respectively;
Cijkl
elastic matrix; reference point；
x*
a point on crack which nearest to x point；
en
outward normal vector of crack element for x* ;
Ni , N j , N k
all nodes set in the grid；
K , K
ED
respectively means Heaviside Functions and crack-tip enriched node sets; nodal displacement vector for parts;
a j and b k
PT
ui
shape function for the node i, j, k；
M
I
AN US
x
CR IP T
displacement and force respectively on boundary ;
strengthen variables for the node; respectively represents the system kinetic energy, strain energy,
CE
T , U , Wd and Ws
AC
damping energy and potential energy force. strain vector; ε ζ stress vector; C is the elastic matrix.
Fb , Fs and F
is defined as enriched force, surface force and concentrated force
respectively; K stiffness matrix; M mass matrix; R load vectors. i
means area of smoothing boundary domain; 3
ACCEPTED MANUSCRIPT
NI
I -shape function node;
n(i ) ( x)
outward normal vector matrix of the boundary; total number of the boundary segment ji ；
M
period of midpoint coordinates (gauss point) in the border ji ；
n jhi 、 l j i
respectively means the outward normal vector and length of ji ;
N e( i )
smoothing area;
I
all nodes of the k set in smoothing zone；
K , K
CR IP T
x GP j
enriched node set respectively by Heaviside function and crack-tip function in a ;
Biu 、B j and Bbk a
AN US
smoothing domain
denote boundary smoothed strain matrix in conventional smoothing domain,
Heaviside function and enriched smoothed strain by crack-tip function;
Aks
k th smoothing domain;
h
separated into x, y ;
M
N seg
boundary number of the smoothing domain;
N gau
ED
Gaussian integral points on the boundary;
,x 、 , y
,x 、 ,x 1
0
a contour around the crack-tip;
AC
Ans
perform transformation between the whole and local coordinates;
CE
Aeffm
2
PT
need to calculate the crack-tip under local coordinate system;
area of the mth element in the equivalent integral domain N effd
;
area of the the mth element in the nth smoothing edge;
K dyn and Kexact
respectively normalized numerical solution and analytical solution of the
dynamic stress intensity factor
First author: Linmei Wu (1987- ), female, postgraduate. Correspondent author: Linmei Wu (1987- ), female, postgraduate. E-mail addresses: [email protected]; [email protected] 1. Introduction 4
ACCEPTED MANUSCRIPT
The modeling of fracture mechanics is important to describe the singular fields near the crack-tip of a solid. However, numerical simulation of dynamic fracture problems remains a challenge in many applications, and is more complex compared to the numerical simulations under static loading conditions. The studies of dynamic fracture mechanics usually give
CR IP T
consideration to the responses of time-dependent loads with inertial effects. However, due to the mathematical complexity, the deriving analytical models/approaches are hard to be used in solving the fracture problems in practice regarding general dynamic element analysis. Numerical methods
AN US
may solve the limitation of analytical model. A number of numerical methods, such as extended finite element method (XFEM) [1] and Meshfree method [2-4], have been applied to solve fracture problems in the last decades. The XFEM has become the most widely used approach to study the arbitrary crack growth. However, the computation accuracy of XFEM is dependent on
M
the mesh regulation as it is based on the level set method (LSM) [5] and partition of unity method
ED
(PUM) [6], which are similar to classical FEM.
PT
To improve the computation accuracy, , Liu et al. [8] proposed a cell-based finite element method (CS-FEM) in 2007, which is based on strain smoothing technique [7]. . Following this,
CE
smoothed finite element method (SFEM) was developed and up-grated to node-based smoothed
AC
finite element method (NS-FEM) [9], edge-based smoothed finite element method FEM (ES-FEM) [10], superconvergent alpha-based finite element method (SαFEM) [11] and face-based smoothed finite element method FEM (FS-FEM) [12] by the same team. Because the integral calculation of the stiffness matrix by SFEM is only at smoothed region boundary, therefore there is a significant reduction of the requirement of
the quality of mesh assemble but it still keeps high accuracy. In
brief, the advantages of SFEM include (1) singular integration will not occur in the calculation of 5
ACCEPTED MANUSCRIPT
crack-tip because no derivation of the shape function is used; (2) integration calculation is focused on smoothed region boundary, instead of a throughout mapping; and (3) insensitive to mesh distortion, which is suitable for the analysis of fracture problems under large deformation. Because of these advantages, smoothed XFEMs are extensively studied recently, particularly on
CR IP T
the fracture problems. Bordas et al. [13] combined CS-FEM with XFEM to contrive cell-based finite element method (CS-XFEM) and applied it to the analysis of the fracture problem of thin plate. The results by the authors showed that the stress intensity factor was affected by the mesh
AN US
distortion and the convergence was better than using conventional XFEM. Chen et al. [14] deduced the edge-based smoothed extended finite element method (ES-XFEM) and applied it to linear elastic crack propagation, and they suggested that the precision and convergence of
M
ES-XFEM were better than CS-XFEM. Jiang et al. [15] applied ES-XFEM to orthotropic material fracture mechanism, which proved again that
ES-XFEM was more accurate than XFEM.
ED
However, the effectiveness and efficiency of applying ES-XFEM to the analysis of dynamic
PT
fracture problems are rarely reported, remaining unknown. The aim of this study to examine the effectiveness and efficiency of ES-XFEM for the
CE
analysis of dynamic fracture problems in 2D elastic solids. To validate the accuracy of the method,
AC
numerical results from DSIFs are compared with the analytical results and other numerical solutions.
2. Solution of dynamic fracture 2.1 Problem description
6
ACCEPTED MANUSCRIPT
The diagram of a plane body with a crack is shown in Fig. 1. The boundary of is composed of the displacement boundary u , stress boundary t and the crack c . The whole domain is recorded as .
c
y u Γu
t
AN US
x
CR IP T
t
Fig.1. Scheme of crack problem
The motion and boundary conditions are given as follows:
ij , j bi ui
(1)
on t
(2)
ij n j 0
on c
(3)
ui ui
on u
(4)
PT
ED
ij n j ti
M
in
CE
where n j is outward normal vector; ij is the Cauchy stress tensor; bi is the body force; ui and ti
t , respectively; is the
AC
are the displacement and stress vectors on boundary u and density, is the domain of the body.
Considering small deformation, the geometric equation is expressed as follows:
ij
1 ui, j u j ,i 2
The constitutive equation is:
7
(5)
ACCEPTED MANUSCRIPT
ij Cijkl kl
(6)
where ij , ij are the stress and strain vectors, respectively; Cijkl are components of the elastic tensor. 2.2 Displacement mode
CR IP T
Compared with classical finite element method, the displacement mode of XFEM takes more functions to describe the models of the discontinuous region. Displacement discontinuity equation of crack elements is expressed demonstrated with enriched functions, such as the
AN US
Heaviside function H x [1]. Above the crack, H x is positively 1; while below the crack it is negative 1:
1 x x* en 0 H x * 1 x x en 0 *
represents a point on crack which is next to
M
In Equation (7), x means reference point; x
(7)
ED
x ; en is the outward normal vector of crack element of x* . In order to capture the singularity and discontinuity near crack-tip, the crack-tip enriched
PT
functions are used. For an isotropic elastic body, the crack-tip enriched functions usually use the
CE
branch crack-tip x [1], which is expressed as follows:
r, 4
r sin 2
r cos 2
r sin sin 2
r cos sin 2
AC
1
where r , is the crack-tip local coordinates in equation (8), and the first term
(8)
r sin 2
describes the discontinuity of crack surface, while the other three terms are employed to improve the singularity of crack-tip [16]. In this way, XFEM displacement mode can be expressed in the following form: 8
ACCEPTED MANUSCRIPT
uh x Ni x ui iI
4
N x H x a N x x b j
j
jK
k
(9)
k
1
k K
where Ni , N j , N k Ni, Nj and Nk are the shape functions for nodes i, j andk, respectively; I is the total node set in the mesh； K and K are Heaviside Functions and crack-tip enriched node sets, are two
CR IP T
respectively; ui are nodal displacement vectors of different parts; a j and b k strengthen variables of nodes.
According to the above XFEM displacement mode, it is not convenient to impose Dirichlet
enriched functions are proposed [17]:
uh x Ni x u i i I
AN US
boundary in the calculation of an unknown quantity . To solve this problem,
N x H x a j
j
K
j
k K
the shifted-basis
4
N kx x b 1
k
(10)
M
where, H x H x H x j ; x x xk . 2.3 The description of dynamic fracture
ED
According to Hamilton’s law of variational principle [17], the governing equation of the
PT
XFEM dynamics can be derived as follows Firstly, Lagrange functional system is expressed as:
L T U Wd Ws
(11)
CE
where T , U , Wd and Ws represents the system kinetic energy, strain energy, damping energy
AC
and potential energy of external force, respectively. The kinetic energy T is: T
1 δ δd 2
T
where δ and δ are displacement and velocity vector, respectively. Strain energy U can be calculated by: 9
(12)
ACCEPTED MANUSCRIPT
U
1 T 1 T ε ζd ε Cεd 2 2
(13)
where ε is strain vector, ζ is the stress vector and C is the elastic matrix. The damping potential energy Wd is:
1 c δ 2
Wd
T
δd
(14)
Potential energy of external force Ws is expressed by:
CR IP T
where c is the viscous damping coefficient.
Ws δ Fb d δ Fs d δ F T
T
T
(15)
AN US
where Fb , Fs and F are defined as enriched force, surface force and concentrated force vector, respectively.
According to the Hamilton's principle, the conservative system is expressed by: t2
Ldt 0
(16)
M
t1
Neglecting of damping, and combining the equations (11) to (15) with the XFEM
Mδ Kδ R
(17)
PT
where δ u a b
ED
displacement mode substitution equation (16), the equation of dynamic XFEM is expressed as:
T
;
K is stiffness matrix; M is the mass matrix and R is load vector.
CE
2.4 Solution of the dynamic equation dynamic equation can be classified as (1) implicit integration
AC
The methods to solve
(Newmark method, Wilson-method, etc.) and (2) explicit integration (center difference method, etc.)[18]. In this works Newmark method was adopted to resolve equation (17). At the n th time step, discrete simulation equations can be expressed as:
M t K δ 2
n
t 2 F K δn 1 tδn 1 1 2 δn 1 2 10
(18)
ACCEPTED MANUSCRIPT
δn δn1 1 δn1 δn t
(19)
1 δn δn 1 δn 1t δn 1 δn t 2 2
(20)
2 where t is the time increment; when 0.5 and 0.25( +0.5) , the Newmark method is
conditions, the effect of damping is therefore not considered
CR IP T
unconditionally stable. This study is to analyze the dynamic fracture of crack under the shock load
3. Edge-based smoothed extended finite element method (ES-XFEM) for crack problems
AN US
3.1 Smoothing domain of ES-FEM
Comparing with conventional finite element, the integral stiffness matrix of SFEM is based on the smoothing region, rather than anelement. Smoothing domain of ES-XFEM is based on edge
1
by
2
M
of an element (as shown in Fig.2). The entire smoothing domain is expressed Ns
, and i
j , i j , where N s represents the
number N s edge. The ABCD ( ) and EFG ( ) are the smoothing domains based on edge m
ED
k
BD ( ) and EF ( ), respectively. k
AC
CE
PT
k
Fig.2. Scheme of edge-based smoothing finite element method 11
ACCEPTED MANUSCRIPT
According to the smoothing operation [11], the boundary strain
i of the element i can be
obtained from the compatible strain ( x) in the region of the smoothing domain i :
i
( x) ( x)d u( x) ( x)d i
s
(21)
i
i
i
Smoothing function i ( x ) meets:
i ( x ) d 1
The smoothing function can be expressed as follows:
1 i Ai 0
x i
(23)
AN US
x
i
(22)
CR IP T
i
where i is area of smoothing boundary domain [2]. i
For smoothing area , the stiffness matrix is: i
i
M
k
B CB T I
J
d
(24)
i
ED
Under the bifurcation theory [2], the smoothed strain matrix B I of the node
I
in domain
CE
PT
i can be expressed as:
b Ix BI 0 b Iy
0 b Iy b Ix
(25)
AC
where
b Ih
1 A( i )
N I ( x )n i ( x )d
( h x, y )
(26)
i
(i ) and N I is I -shape function node; n ( x) is outward normal vector matrix of the boundary
i , which is represented as follows:
12
ACCEPTED MANUSCRIPT nx( i ) n(i ) ( x) 0 n(yi )
0 n(yi ) nx( i )
(27)
As linear displacement field is used at in , a Gaussian point is applied at . i
i
i
Equation (26) can be further expressed in the algebra form:
1
M
N A i
j 1
I
( x GP )n(jhi ) l ji j
( h x, y )
(28)
CR IP T
b Ih
i GP where M is the total number of the boundary segment j ; x j means the period of midpoint
i i i coordinates (gauss point) in the border j ; n jh and l j are the outward normal vector and
AN US
i length of j , respectively. Smoothing region of the stiffness matrix can be further expressed as
follows:
k
i
B I C B J Ai T
(29)
M
From equation (26), the stiffness matrix of the shape function does not require partial
ED
derivatives. The whole stiffness matrix K assembling each region stiffness matrix has beenset,
PT
which is:
N e( i )
K = k
i
(30)
i 1
CE
where N e( i ) is the smoothing area.
AC
3.2 Equations of ES-XFEM ES-XFEM is Deriving from equation (17), the assembled mass matrix and load vector for
ES-XFEM do not involve strain, therefore dynamic extension of finite element method (FEM) is needed. Similar to XFEM, the smoothed strain of the displacement mode of ES-XFEM in the smoothing domain can be written by: 13
ACCEPTED MANUSCRIPT
k Biu ui iI
jK
B aj a j
4
B b b
kK
1
k
k
(31)
in the smoothing zone； K and K mean enriched
where I means all nodes of the set k
node set by Heaviside function and crack-tip function in a smoothing domain, respectively; Biu ,
CR IP T
B aj and Bbk denote boundary smoothed strain matrix in conventional smoothing domain,
Heaviside function and enriched smoothed strain by crack-tip function, respectively [14]. To enrich the matrix, its specific expressions are as follows:
0 biyr bixr
where
ks
nh x Ni x d
1 Aks
ks
nh x Ni x H x d
ED
biha
1 Aks
1 Aks
ks
nh x Ni x x d
(33a)
(33b)
1 , 2 , 3 , 4
(33c)
PT
bihb
(32)
M
bihu
r u , a, b
AN US
bixr Bir 0 biyr
where Aks is the number k smoothing domain, h is separated into x, y .
AC
CE
Equation (33) can be presented in numerical forms:
bihu
biha
1 Aks
1 b s Ak b ih
1 Aks
N gau wm,n Ni xm n, nh xm n, m 1 n 1 N seg
(34a)
N gau wm,m Ni xm,n H xm,m nh xm ,n m 1 n 1
(34b)
N gau wm,m Ni xm,n xm,m nh xm,n m 1 n 1
(34c)
N seg
N seg
14
ACCEPTED MANUSCRIPT
where N seg is the boundary number of the smoothing domain; N gau represents Gaussian integral
M
AN US
CR IP T
point on the boundary.
Fig.3. Schematic diagram of edge-based smoothed extended finite element method
ED
The stiffness matrix kij is derived from the smoothed strain matrix:
CE
PT
B u T CB u d j ks i T kij s Bia CB uj d k B b T CB u d j ks i
B
B
CB bj d T T a a a u ks Bi CB j d ks Bi CB j d b T a u T u ks Bi CB j d ks Bi CB j d ks
u T i
CB aj d
ks
u T i
(35)
AC
Due to smoothed strain matrix in smoothing domain is constant, the equation (35) can be
written as
B u T CB u As j k i a T kij Bi CB uj Aks b T u s Bi CB j Ak
B B B
u T i
CB aj Aks
a T i
CB aj Aks
b T i
CB aj Aks
15
B B B
CB bj Aks T b s CB j Ak T CB bj Aks
u T i
a i b i
(36)
ACCEPTED MANUSCRIPT
e Using m to express the element mass matrix set:
mijuu mije mijau mijbu
mijub mijau mijbb
mijua mijaa mijba
(37)
where
mijuu e Ni N j d T
(38a)
mijua e Ni
T
mijub e Ni
T
HN d
mijaa e HNi
T
T
N j d,
(38c)
(38d)
i
AN US
mijbb e Ni
1,2,3,4
HN d
mijab e HNi N j d,
(38b)
j
N j d,
T
CR IP T
1,2,3,4
1,2,3,4
(38e) (38f)
Using ri e to express the element nodal load vector components:
rib1
rib 2
rib3
rib 4
(39)
ED
where
ri a
M
rie riu
(40)
CE
PT
riu e Ni Fb d Ni Fs d Ni F a ri e N i HFb d N i HFs d N i HF b ri e Ni Fb d Ni Fs d Ni F 1, 2,3, 4
AC
4. Integration implementation of ES-XFEM Because the stiffness matrix kij is derived from integration on the smoothing boundary, the
conventional XFEM integral solution is no longer applicable. A new integral formula has been proposed in [9], which is expressed as follows: (1) As shown in Fig.3, each boundary employs a Gaussian point in the conventional smoothing domain (such as ABCD). 16
ACCEPTED MANUSCRIPT
(2) Each boundary employs a Gaussian point in crack-free enhanced Heaviside function domain. (3) Each boundary adopts 3~5 Gaussian points in crack-free and crack-tip enriched smoothing domain.
CR IP T
(4) In the smoothing domain where the crack run through, it can be divided into two smoothing subdomain and one Gaussian point is used on the each smoothing boundary; As shown in Fig.4, irregular smoothing subdomain is usually divided into a triangle sub-domain in order to
PT
ED
M
AN US
facilitate programming.
CE
Fig.4. Partition of split smoothing domains in the ES-XFEM[14]
AC
(5) As shown in Fig.5, crack is extended to smoothing boundary of the region and then divide the crack-tip domain into multiple sub-domains, generally three to five ones. The boundary of the smoothing subdomain near the crack-tip need 5 Gaussian points, and other boundary needs 1 ~ 3 points.
17
ACCEPTED MANUSCRIPT
(6) As shown in Fig.6, in the integral of enriched function, the partition rule of smoothing sub-domain is same as above, but Gaussian points are not taken on the smoothing sub-domain boundary and while adopt the interior points of the smoothing sub-domain. This is because if the smoothing domain on the boundary points is taken as a Gaussian point, boundary AB, CD, EF,
CR IP T
GH are coincidence with crack surface and then cannot satisfy the accuracy of calculation. In general, the crack-tip items of enriched function utilize about 3~5 Gaussian points, otherwise enhanced Heaviside takes one Gaussian point. Sub-domain integration method is employed to
AN US
calculate both Mass matrix and stiffness matrix in conventional XFEM. It is needed to employ one Gaussian point in the internal unit of unenhanced nodes, 5 ~ 7 points in crack-tip sub-domains and
AC
CE
PT
ED
M
3 points in other domains.
18
ACCEPTED MANUSCRIPT
CR IP T
Fig.5. Partition of crack-tip smoothing domains in the ES-XFEM
AN US
Fig.6. Scheme of Gaussian point in split and crack-tip element
5. Dynamic stress intensity factor
M
Stress intensity factor can be calculated by displacement extrapolation method [19], virtual crack extension method[20], virtual crack closure legal [21], and interaction integral method[22].
ED
In this work, interaction integral method has been used. Considering the two dimensional crack,
PT
crack is described as a line, while a local coordinate is set up at crack-tip. As shown in Fig.7, 0 means a contour around the crack-tip.
CE
Contour integrals J near the crack-tip can be expressed as:
AC
J lim
0 0
0
W
1i
ij u j ,1 ni d lim
0 0
1 2
0
W
where W ij ij , which is the strain energy density [18].
19
1i
ij u j ,1 mi d
(43)
CR IP T
ACCEPTED MANUSCRIPT
AN US
Fig.7. Scheme of J-integral domain
Equation (43) cannot be directly applied to numerical calculation and then adopted by a 0 1 closed domain s 1 0 . The set weight factor q is 1 in , and is 0 in ,
M
Equation (43) can be rewritten as:
J ij u j ,1 W 1i qmi d s
2 j u j ,1d
ED
(44)
Using the divergence theorem, the integral form of the J integral is
PT
J ij u j , 1 W i 1 q i d,A ij u j i , 1W qd A , 1 A
A
(45)
CE
By incorporating the balance equations and geometric equation: J ij u j ,1 W 1i q,i dA uu j ,1 ui u,1 qdA A
A
(46)
Choose two independent equilibrium states: state 1 ( ij , ij , u i(1) ) as the real state, state 2
AC
(1)
(1)
( ij , ij , u i( 2 ) ) as the auxiliary condition. The third equilibrium is achieved by linear ( 2)
( 2)
1,2 superposition of the two equilibrium, with J integral can be rewriten by J :
J 1 , 2 J 1 J 2 I
1, 2
where I 1,2 means the interaction integral of state 1 and state 2: 20
(47)
ACCEPTED MANUSCRIPT
I 1,2 1u j 2,1 2u j1,1 W 1,21i q,i dA u 1u j2,1 qdA A
(48)
A
where W 1,2 ij 2 ij1 ij1 ij 2 is the interaction of strain energy density. According to [1], the two-dimensional equation of energy release rate can be expressed as [1]:
E with E E 1 v 2
K I2 K II2 E E
(49)
plane stress plane strain
.
CR IP T
J
For the third equilibrium, equation (49) can be reformed as follows:
J
1
J
2
2 K I1 K II 2 K I 2 K II1
AN US
J
1,2
E
(50)
Comparing with Equation (48) and Equation (50), it is expressed as:
2 K I1 K II 2 K I 2 K II1
E
(51)
M
I
1,2
ED
State 2 is chosen as asymptotic field of mode I and II, then:
E 1,mod eI I 2
K II1
E 1,mod eII I 2
(52)
PT
K I1
CE
As shown in Fig.7 (b), the interaction integrals are based on the smoothing boundary,no longer
AC
based on the element integrals for ES-XFEM, so equation (49) is rewritten as:
I
1,2
d Neff
m 1
m Aeff
ij
d Neff
u j ,1 ij u j ,1 ij ij 1i q,i dA
1 2
2 1
1 2
m 1
m Aeff
u 1u j2,1 qdA
(53)
m d where Aeff is area of the number m element in the equivalent integral domain N eff .
In equation (53), the first part of the right side is related to the displacement derivative which is based on integral of the smoothing domain; the second item of right side is independent of the
21
ACCEPTED MANUSCRIPT
displacement derivative term, which is based on the element integral.
Equation (53) is depicted
as:
I
1,2
d Neff
3
m 1 n 1
Ans
ij
d Neff
u j ,1 ij u j ,1 ij ij 1i q,i dA
1 2
2 1
1 2
m 1
m Aeff
u 1u j2,1 qdA
(54)
s
CR IP T
where An is area of the the m-th element in the n-th smoothing edge. 6. Numerical examples
6.1 Example 1: dynamic fracture of edge crack in the semi-infinite domain
Because of computational limits, a limited domain is modeled for analysis, as shown in Fig.8.
AN US
The size of the limited plate are L = 10 m, H = 2 m, crack length a = 5 m；The elastic constants of materials are: E = 210 GPa, v = 0.3, ρ =8000 kg/m3. σ0 = 500MPa. The top of the plate is subject to a distributed step load σ0 (t) = σ0 H(t). When t>0, H(t)=1. Computation time is 1×10-3 s, time step
M
is 2×10-5 s and the number of iterations is 50 times. The regular meshes use 65×25 triangular
ED
elements, as shown in Fig. 9 (a). In order to conform rough requirement of mesh in the smoothed expansion finite element method, different irregular meshes are adopted. The irregular mesh is
CE
PT
generated randomly by [8]:
x ' x x rd ir ' y y y rd ir
(55)
AC
Where x , y respectively represents node coordinates of initial rules mesh, as shown in Fig.9.
(a) ; Accordingly, x , y is the length of element respectively in the x, y direction； rd is random numbers in the domain[-1, 1]; ir is irregular factors. Regular mesh ( ir 0.00 )and irregular mesh( ir 0.12 ) are shown is shown in Fig. 9.
22
ACCEPTED MANUSCRIPT
CR IP T
Fig.8. Geometry of a semi-infinite crack
Since reflection occurs when the stress wave meets the crack-tip, calculation process does not consider the impact of reflection, and the dynamic stress intensity factor analytic formula is valid
H/cd , cd E(1 v ) / ( (1 v ) (1 2v ) ) ).
AN US
just as time t 3tc ( tc
Regardless of the crack propagation, dynamic stress intensity factor analytic is expressed as [23]:
2 0 cd t (1 2v) 1 v
M
K Idyn (t )
(56)
ED
The relative error of normalized dynamic stress intensity factor is defined as follows:
PT
err
K dyn Kexact K exact
100%
(57)
CE
where K dyn and K exact are normalized numerical solution and analytical solution of the dynamic stress intensity factor, respectively.
AC
In the following simulation figures, K Idyn / 0 H1/2 , t / tc are the dynamic stress intensity
factor of the normalized mode I and the normalized time which the dimension is 1, respectively.
23
AN US
Fig.9. Computational meshes
CR IP T
ACCEPTED MANUSCRIPT
6.1.1 Convergence analysis
It is researched on convergence of different computational meshes. As can be seen from
M
Fig.10., when t 2.0tc , relative error of dynamic stress intensity factor by 43×19 meshes
ED
decreased better than 34 × 14 meshes. 53 × 21 meshes and 65 × 25 meshes have similar
AC
CE
PT
performance with convergence. The following analysis is based on 65×25 meshes.
Fig.10. (a) normalized dynamic stress intensity factor; (b) relative error of dynamic stress intensity factor at various meshes 6.1.2 Selection of crack-tip smoothing subdomain 24
ACCEPTED MANUSCRIPT
As the crack-tip region is divided into different number of smoothing subdomain and the influence of subdomains number on the dynamic stress intensity factor is analyzed. As can be seen from Fig.11.(a), when number of smoothing subdomain n ssd 2 , the difference of dynamic stress intensity factor in normalization mode I is not obvious. However, as shown in Fig.11.(b),
CR IP T
when n ssd 3 , relative error of the normalized dynamic stress intensity factor in mode I tends to be stable. Therefore, in order to improve the computation efficiency and guarantee the calculation accuracy, in the crack-tip smoothing domain sequence number should be used as n ssd =3 and
AN US
employs five Gaussian points at each boundary. The smoothing domain with crack across do not
PT
ED
M
need partition with one Gaussian point is taken; other area also take one Gaussian point.
CE
Fig.11. (a) normalized dynamic stress intensity factor; (b) relative error of dynamic stress intensity factor with different smoothing subdomain
AC
6.1.3 Accuracy and mesh distortion As shown in Fig.12, it is compared the ES-XFEM results with the analytical solution, finite
element method of singular boundary smoothing [24] (sES-FEM) and extended finite element method (XFEM). It is demonstrated that when t 1.5tc , numerical error of three methods are all relative large; when t 1.5tc , the results of three kinds of numerical methods are very close to
25
ACCEPTED MANUSCRIPT
the analytical solution; when t tc , the error of the numerical solution is the largest and the maximum error is XFEM. With classical finite element method, mesh distortion may affect the computational accuracy. ES-XFEM is no longer restricted, while Fig.13. shows the calculation results of different mesh
CR IP T
irregular factors. It is illustrated that the ES-XFEM results affected by the quality of the mesh are very small. When the irregular factor takes 0.12, the mesh is extremely irregular. There are 74 angles greater than 120o , 372 angles less than 30o and 39 angles less than 20o (as shown in
AN US
Fig.9 (b)). In this case, the ES-XFEM can still get a stable result as the time 1.8tc t 3.0tc , the maximum relative error of the dynamic stress intensity factor in the normalized mode I is less than
PT
ED
M
3%.
intensity factor comparing with different method
AC
CE
Fig.12. (a) normalized dynamic stress intensity factor (b) relative error of dynamic stress
26
CR IP T
ACCEPTED MANUSCRIPT
Fig.13. (a) normalized dynamic stress intensity factor (b) relative error of dynamic stress intensity factor with different irregular factor
AN US
6.2 Example 2: dynamic fracture of the bevel edge crack plate
In order to prove the efficiency and stability of ES-XFEM for dynamic fracture problems, we further investigative the mixed-mode dynamic fracture. An inline edge crack of plate is shown in Fig.14, its geometry sizes are L = 44 mm, H = 32mm, D = 16 mm, crack length is a = 22.63 mm,
M
and the incline angle of crack α = 45o. Material parameters are set as: modulus of elasticity E =
ED
29.4 GPa, Poisson’s ratio v = 0.286, the density ρ=2450 kg/m3. The left, top and bottom of plate
PT
bear restriction of chain support in normal direction, the right end bears a step load distributionσ0 (t) = σ0 H(t). Computation time is 3×10-5 s, Time step is 6×10-7s and the number of iterations is
CE
50 times. Two kinds of computational meshes (66×48) with different quality (regular and irregular mesh) are revealed in Fig.15.
AC
Calculation results in the figure below shows that, K Idyn / 0 πa
1/2
, K / dyn II
πa
1/2
0
are the normalized dynamic stress intensity factor of mode I and II, respectively. In Fig.16, analyzes and compares the normalization of the dynamic stress intensity factor under different numerical methods with the regular mesh and the results of three methods are
27
ACCEPTED MANUSCRIPT
consistent and the peak value of dynamic stress intensity factor in normalized mode I is higher than that in normalized mode II. Fig.14 presents ES-XFEM results under different mesh distortion. When the irregular factor takes 0.2, maximum of the element angle is 154.3°, minimum is 10.6°. There are 172 elements
CR IP T
with maximum angles greater than 120°, 485 elements with the minimum angles less than 30°and 30 elements with the minimum angles less than 20°. As shown in Fig.17, mesh distortion effected on the dynamic stress intensity factor is small. It is demonstrated again that ES-XFEM is not
AN US
sensitive to the mesh quality, which will help to analyze large deformation of crack in propagation
AC
CE
PT
ED
M
process.
28
ACCEPTED MANUSCRIPT
PT
ED
M
AN US
Fig.15 Computational meshes
CR IP T
Fig.14 Geometry of a rectangular plate with an inclined crack
edge crack
AC
CE
Fig.16 Comparison of the normalized mixed-mode DSIFs for a rectangular plate with a slanted
29
CR IP T
ACCEPTED MANUSCRIPT
Fig.17 The normalized mixed-mode DSIFs for a rectangular plate with a slanted edge crack
AN US
under different mesh irregular factor
7 Conclusions
This paper applies the edge based smoothed extended finite element (ES-XFEM) to solve
M
stationary dynamic fracture problems subjected to impact loads. Dynamic stress intensity factor
ED
(DIF) is computed with the domain-form of the interaction integral in terms of the strain
PT
smoothing method. Two case studies including the mode I and mixed crack prove that the dynamic ES-XFEM is highly precise, and insensitive to mesh distortion. It is also proved that
3 , the relative error of dynamic stress intensity factor become
CE
when smoothing subdomain n ssd
AC
stable.
References [1] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. International journal for numerical methods in engineering. 1999;45:601-20. [2] Amiri F, Anitescu C, Arroyo M, Bordas S, Rabczuk T. XLME interpolants, a seamless bridge 30
ACCEPTED MANUSCRIPT
between XFEM and enriched meshless methods. Computational Mechanics. 2014;53:45-57. [3] Amiri F, Millán D, Shen Y, Rabczuk T, Arroyo M. Phase-field modeling of fracture in linear thin shells. Theoretical and Applied Fracture Mechanics. 2014;69:102-9. [4] Nguyen-Xuan H, Liu G, Bordas S, Natarajan S, Rabczuk T. An adaptive singular ES-FEM for
CR IP T
mechanics problems with singular field of arbitrary order. Computer Methods in Applied Mechanics and Engineering. 2013;253:252-73.
[5] Osher S, Sethian JA. Fronts propagating with curvature-dependent speed: algorithms based on
AN US
Hamilton-Jacobi formulations. Journal of computational physics. 1988;79:12-49.
[6] Duarte CA, Oden JT. An h-p adaptive method using clouds. Computer methods in applied mechanics and engineering. 1996;139:237-62.
M
[7] Chen J-S, Wu C-T, Yoon S, You Y. A stabilized conforming nodal integration for Galerkin mesh-free methods. International journal for numerical methods in engineering. 2001;50:435-66.
ED
[8] Liu G, Dai K, Nguyen T. A smoothed finite element method for mechanics problems.
PT
Computational Mechanics. 2007;39:859-77. [9] Liu G, Nguyen-Thoi T, Nguyen-Xuan H, Lam K. A node-based smoothed finite element
CE
method (NS-FEM) for upper bound solutions to solid mechanics problems. Computers &
AC
structures. 2009;87:14-26. [10] Liu G, Nguyen-Thoi T, Lam K. An edge-based smoothed finite element method (ES-FEM) for static, free and forced vibration analyses of solids. Journal of Sound and Vibration. 2009;320:1100-30. [11] Liu G, Nguyen-Xuan H, Nguyen-Thoi T, Xu X. A novel Galerkin-like weakform and a superconvergent alpha finite element method (SαFEM) for mechanics problems using triangular 31
ACCEPTED MANUSCRIPT
meshes. Journal of Computational Physics. 2009;228:4055-87. [12] Nguyen-Thoi T, Liu G, Lam K, Zhang G. A face-based smoothed finite element method (FS-FEM) for 3D linear and geometrically non-linear solid mechanics problems using 4-node tetrahedral elements. International Journal for Numerical Methods in Engineering. 2009;78:324.
in FEM and XFEM. Computers & structures. 2010;88:1419-43.
CR IP T
[13] Bordas SP, Rabczuk T, Hung N-X, Nguyen VP, Natarajan S, Bog T, et al. Strain smoothing
[14] Chen L, Rabczuk T, Bordas SPA, Liu G, Zeng K, Kerfriden P. Extended finite element
AN US
method with edge-based strain smoothing (ESm-XFEM) for linear elastic crack growth. Computer Methods in Applied Mechanics and Engineering. 2012;209:250-65.
[15] Jiang Y, Tay T, Chen L, Sun X. An edge-based smoothed XFEM for fracture in composite
M
materials. International Journal of Fracture. 2013;179:179-99.
[16] Ahmed A, Auricchio F. Extended finite element method (XFEM)-modeling arbitrary
ED
discontinuities and failure analysis. Research degree thesis. 2009.
PT
[17] Stąpór P. Application of xfem with shifted-basis approximation to computation of stress intensity factors. Archive of Mechanical Engineering. 2011;58:447-83.
CE
[18] Hughes TJ. The finite element method: linear static and dynamic finite element analysis:
AC
Courier Corporation; 2012. [19] Shih C, Lorenzi Hd, German M. Crack extension modeling with singular quadratic isoparametric elements. International Journal of Fracture. 1976;12:647-51. [20] Parks DM. A stiffness derivative finite element technique for determination of crack tip stress intensity factors. International Journal of fracture. 1974;10:487-502. [21] Rybicki EF, Kanninen M. A finite element calculation of stress intensity factors by a 32
ACCEPTED MANUSCRIPT
modified crack closure integral. Engineering Fracture Mechanics. 1977;9:931-8. [22] Yau J, Wang S, Corten H. A mixed-mode crack analysis of isotropic solids using conservation laws of elasticity. Journal of Applied Mechanics. 1980;47:335-41. [23] Menouillard T, Belytschko T. Dynamic fracture with meshfree enriched XFEM. Acta
CR IP T
mechanica. 2010;213:53-69. [24] Liu P, Bui T, Zhang C, Yu T, Liu G, Golub M. The singular edge-based smoothed finite element method for stationary dynamic crack problems in 2D elastic solids. Computer Methods in
AC
CE
PT
ED
M
AN US
Applied Mechanics and Engineering. 2012;233:68-80.
33