- Email: [email protected]
The continuous-discontinuous cellular automaton method for elastodynamic crack problems
Accepted Manuscript The continuous-discontinuous cellular automaton method for elastodynamic problems Fei Yan, Peng-Zhi Pan, Xia-Ting Feng, Shao-Jun Li PII: DOI: Reference:
S0013-7944(18)30758-6 https://doi.org/10.1016/j.engfracmech.2018.10.025 EFM 6199
To appear in:
Engineering Fracture Mechanics
Received Date: Revised Date: Accepted Date:
28 July 2018 11 October 2018 24 October 2018
Please cite this article as: Yan, F., Pan, P-Z., Feng, X-T., Li, S-J., The continuous-discontinuous cellular automaton method for elastodynamic problems, Engineering Fracture Mechanics (2018), doi: https://doi.org/10.1016/ j.engfracmech.2018.10.025
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.
The continuous-discontinuous cellular automaton method for elastodynamic problems Fei Yan1, Peng-Zhi Pan1,2,3*, Xia-Ting Feng3, Shao-Jun Li1 1. State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China 2. University of Chinese Academy of Sciences, Beijing 100049, China 3. Northeastern University, Shenyang 110819, Liaoning, China
Abstract: In this paper, the recently proposed continuous-discontinuous cellular automaton method (CDCA) is further developed for elastodynamic problems in 2D elastic solids. A continuity-to-discontinuity enriched function technique is developed for elastodynamics and combined with a cellular automaton formulation; the cellular automaton theory for elastodynamics, in which fast adaptive updating rules for the cell and its neighbors are developed, and the time domain CDCA are proposed. In this method, the cracks are independent of the integral meshes, and no assembled global matrix is needed for the whole process. Time-dependent equations for the cells are discretized by the Newmark method, and rapidly solved by the cellular automaton method. To analyze the fracture behaviors of elastic solids with dynamic loads, mixed-mode dynamic stress intensity factors are obtained by the interaction integral method. Some numerical examples for elastodynamic problems are studied to validate the accuracy and efficiency of the present method, and solutions from analytical method and some other methods are employed to show the high accuracy of the CDCA.
Keywords: the time domain continuous-discontinuous cellular automaton; cellular automaton method; fast updating scheme; dynamic stress intensity factor; elastodynamic problem
*
Corresponding author. Email address: [email protected] (P.- Z. Pan) 1
1. Introduction The modeling of dynamic crack problems is of great importance in the simulation of structure failure. Because of the mathematical difficulties of time-dependent loading, analytical methods are often not feasible to solve general dynamic crack problems in practical engineering, and numerical methods are thus developed for solving such problems. However, accurate and highly efficient numerical simulations of dynamic fracture problems remain a challenging task for many practical applications in engineering, and many significant numerical methods have been developed for dynamic crack problems in the past few years. As the most widely used simulation method, the finite element method (FEM) is also widely used to solve dynamic fracture problems[1-2], in which a singular element is developed. Based on cell-based smoothing domains, a singular edge-based smoothed finite element method (sES-FEM) is proposed for stationary dynamic crack problems in 2D elastic solids[3]. FEM, with remeshing and a stable numerical scheme, has also been used for dynamic crack propagation[4], but there are some disadvantages of the FEM such as the need for discretizing the entire domain; this is especially disadvantageous when remeshing for crack propagation, and inaccuracies arise in the form of stress concentrations. As a boundary-type numerical method, the boundary element method (BEM) has been widely applied in dynamic crack problems, such as the dynamic analysis of crack by Chirino and Dominguez[5], dynamic analysis of nonlinear systems by Konroni and Beskos[6] and Burczynski and Adamczyk[7], inelastic dynamic analysis by Ahmad and Banerjee[8] and Kontoni and Bekos[9], and dynamic elastic-plastic structural problem by Panzeca et al. [10] and Carrer and Telles[11]. To overcome the required domain integrations, dual reciprocity method (DRM) is employed in BEM for solving the elasto-dynamics, such as transient elastodynamic and free vibration analysis by Agnantiaris et al. [12-13] and Samaan and Rashed[14] and dynamic stress intensity factor (DSIF) computations by Dominguez and Gallego[15]. Recently, some techniques have been developed to avoid remeshing and are 2
suitable for crack propagation, i.e., the extended finite element method (XFEM). Gao and Klein[16] first developed the idea of loss of hyperbolicity for analyzing dynamic crack propagations. Belytschko et al.[17] developed a new discontinuous enrichment for dynamic crack propagation based on the loss of hyperbolicity. Further, Chessa and Belytschko[18] presented a locally enriched space-time XFEM for solving hyperbolic problems. Menouillard et al.[19] introduced a lumped mass matrix for enriched element to allow the use of pure explicit formulation in XFEM. Elguedj et al.[20] proposed a mass lumping for arbitrary enrichment functions. Rethore et al.[21] employed an energy-conserving scheme for dynamic crack growth, and Wu et al.[22] developed an edge-based smoothed extended finite element method for dynamic fracture analysis. Duarte et al.[23] proposed a generalized finite element method (GFEM) for 3D dynamic crack propagation, and Li et al.[24] combined moving least squares and numerical manifold method (NMM) for impact analysis of cracked bodies. In recent years, meshless methods have also been developed for dynamic crack problems due to their non-meshing, such as element-free Galerkin method (EFG), dynamic fracture and wave propagation analysis by Belytschko and Tabbara[25] and further for 3D dynamic propagation by Krysl and Belytschko[26] and Lu et al. [27], mixed-mode stress intensity factors evaluation by Wen and Aliabadi[28], the cracking particles method (CPM) employed by Rabczuk et al. [29] to simulate dynamic fracture, the hybrid boundary node method (HBNM) and dual reciprocity hybrid boundary node method[30-33], free vibration by Yan et al. [34], and transient dynamic fracture problem by Miao et al. [35]. To overcome the defect of remeshing and to easily consider the plastic fracture propagation, based on the cellular automaton method[36] and discontinuous enrichment functions, Yan et al.[37] proposed a continuous-discontinuous cellular automaton method (CDCA), and employed this method to solve the static multiple cracks propagation[38], cohesive crack propagation[39], contact crack[40] and so on. The advantages of this method are: first, the cracks are independent of the calculation mesh, so no remeshing is needed with the cracks propagation; second, no assembled 3
global matrix is needed in the whole calculation, the computing memory can be greatly saved, and the calculation is located at the local cell, so the local properties can be easily considered; finally, with the fast updating method, the calculating efficiency can be greatly improved and can be easily applied for some large practical engineering. Based on the continuous-discontinuous cellular automaton method for static fracture problems, a time domain CDCA is proposed in this paper. First, the continuity to discontinuity enriched function technique is developed for dynamic fracture problems; second, a time domain cellular automaton theory for elastodynamics is proposed, in which fast adaptive updating rules based on time domain cellular automaton for the cell and its neighbors are developed; then, time-dependent equations for the cells are discretized individually by the Newmark method, and rapidly solved by the cellular automaton method; finally, to analyze the fracture behaviors of elastic solids with dynamic loads, mixed-mode dynamic stress intensity factors are obtained by the interaction integral method. Because of the usage of time domain cellular automaton theory and fast adaptive updating rules, the cracks are independent of the integral meshes, and no assembled global matrix is needed for the whole process. The computing memory can be greatly saved, and the calculation is located at the local cell, so the local properties can be easily considered. In addition, the calculating efficiency can be greatly improved. Finally, some numerical examples of elastodynamic problems are studied to validate the accuracy and efficiency of the present method, and solutions from analytical method and some other methods are employed to show the high accuracy of the CDCA.
2. The continuous-discontinuous cellular automaton 2.1 Governing equation of elastodynamics Consider a body
with an initial crack
in the state of dynamic equilibrium,
the governing elastodynamic equation can be given as σ f b u 4
(1)
in which σ is the stress tensor, f b is the body force, is the density of the
is the acceleration vector, which is on the following boundary material, and u conditions: u(x, t ) u (x, t )
on u
(2)
σ n f t
on t
(3)
σ n 0
on c
(4)
in which u , t and c are the displacement, traction and crack boundaries, respectively, and f t is the external traction vectors. The initial conditions can be written as
u(x, t 0) u (0)
(5)
u (x, t 0) u (0)
(6)
in which u (0) and u (0) are the initial value vector of displacement and velocity, respectively. According to variational theory, the variational formulation of this boundary value problem of Eq. (1) can be given as
ud σ εd f b ud f t ud u
(7)
2.2 Continuous-discontinuous theory for cellular automaton Based on static continuous-discontinuous cellular automaton theory, the Heaviside function and the exact near-tip asymptotic field functions are employed in this method to model the strong discontinuity and high gradient stress field around the crack tip. Then, the shape function at time t of the present method can be given as[37-39] n
m
u h (x, t ) N j (x)u j (t ) N k (x)H ( ) H ( k ) a k (t ) j 1 k 1 kP
N i (x) Fl (x) Fl (x i ) b li (t ) i 1 l 1 t
nf
iT
5
(8)
in which N j (x) is the traditional FEM shape function, u j (t ) is the nodal displacements at time t , a k (t ) is a vector of additional degree of nodal freedom for modeling strong discontinuity t , and b li (t ) is a vector of additional degrees of nodal freedom for modeling the crack tip stress field at time t . The signed function is chosen as the Heaviside enrichment function, which is given as
1 0 H ( ) sign ( ) 1 0
(9)
in which is the value of level set function, which can be seen in[37-39] and Fl (x) is the exact near-tip asymptotic field functions, which can be given as[37-39]
{Fl ( x), l 1 4} r sin( ), r cos( ), r sin( ) sin( ), r sin( ) cos( ) 2 2 2 2
(10)
in which r , are the polar coordinate of the crack tip polar coordinate system. To ensure the interpolation property is satisfied, the shifting amendments of H ( ) H ( k ) and Fl (x) Fl (xi ) are employed in Eq. (8). Following Eq. (8), the
acceleration of nodes can be given as n
m
h (x, t ) N j (x)u j (t ) N k (x)H ( ) H ( k ) ak (t ) u j 1 k 1 kP
l (t ) N i (x) Fl (x) Fl (x i ) b i i 1 l 1 t
nf
(11)
iT
Then, the discretized form of Eq. (7) using Eqs. (8) and (11) can be rewritten as
(t ) Ku(t ) F(t ) Mu
(12)
(t ) denote the vector of nodal parameters and its second time in which u(t ) and u derivative at time t , respectively, which include displacement u j (t ) , enrichment
j (t ) and second time derivative of a k (t ) . degrees of freedom a k (t ) , acceleration u The stiffness matrix K , mass matrix M and the external loading vector F 6
can be defined as K uu K ij [ ijau K ij M uu M ij [ ijau M ij
K ua ij ] K ijaa
(13)
M ua ij ] M ijaa
(14)
Fi (t ) Fiu (t ) Fia (t )
T
(15)
in which the stiffness matrix K and the external loading vector F can be seen in references[37-39], and the mass matrix can be given as Muu ij N i N j d
M ua ij N i N j [ H ( ) H ( j )]d
M ijaa Ni [ H ( ) H (i )]N j [ H ( ) H ( j )]d
(16) (17) (18)
in which N i is the traditional FEM shape function for node i . If the enrichment is the exact near-tip asymptotic field functions, Eq. (18) can be rewritten as M bb ij N i [ Fl (x) Fl (x i )]N j [ Fl (x) Fl (x j )]d
(19)
2.3 Cellular automaton theory for dynamic problem By using discontinuous cellular automaton, the equilibrium state of the object can be obtained through the self-organization phenomenon formed by the one-another transfer of information between cells. Its advantages are: first, the localization property of the calculating structure can be easily considered, and the behavior of the cell is thought to be essentially local; in other words, the state of one cell is only determined by the states of itself and its neighbors. Second, there is no need to assemble the overall matrix, and the different degrees may bring some difficulties for the assembling operation for the enriched cells. Finally, the large-scale simulation can be performed because of its easy implementation property of the parallel algorithm. 7
The CDCA model is composed of cell, cell space, cell state, crack, neighborhood, updating rules and so on, and the relation between those components can be seen in references[37-39]. Crack, cell, neighbors and their relations can be seen in Fig. 1.
strong discontinuous cell
neighborhood
crack tip cell
crack tip element through element
cell
crack continuous element
neighborhood
continuous cell
Fig. 1 Cell, neighbors, crack and their relations 2.3.1 Cell space and its states for elastodynamics Same as the static cellular automaton model, the 2D cell space can be rectangular, triangular, hexagon and so on. To describe the cell states, a series of physical and mechanical properties must be defined to determine its different states. A CDCA model
is
composed
of
the
degree
values
vector
of
nodal
freedom
u(t ) {u(t ), a(t ), b(t )} , in which u(t ) is the traditional degree of nodal freedom, a(t ) is the Heaviside enriched degree of nodal freedom and b(t ) is the crack tip field
function
enriched
degree
of
nodal
freedom,
time
derivative
u (t ) {u (t ), a (t ), b (t )} , which includes the velocity of the node and time derivative of
the enriched nodal freedom, a(t ) or b(t ) , and the second time derivative (t )} , which includes the acceleration of node and second time (t ) {u (t ), a(t ), b u
derivative of enriched nodal freedom a(t ) or b(t ) . 8
In addition, the cell state includes current time t t k , time step increment t , the material property of thickness h , Young’s modulus E , Poisson’s ratio and fracture toughness K IC ; the cell force, strain and stress state can be given in the cell nodal force vector f (t ) {fu (t ), fa (t ), fb (t )} , in which the subscripts u , a and b represent the traditional, Heaviside enriched and exact near-tip asymptotic field functions enriched by the degrees of nodal freedom, respectively; additionally, the cell state includes the cell strain ε , cell stress σ , and crack location information, i.e., the distance function , and so on. Based on these theories, the cell state can be given in Fig. 2.
displacement
u (t )
velocity u (t )
strain
(t )
stress (t )
time derivative of enriched freedom
a(t ) / b(t )
cell state
enriched nodal freedom
t tk
enrihed nodal force
a(t) / b(t)
Fa / b (t) traditional nodal force
acceleration
Fu (t)
u(t )
second time derivative of enriched freedom
distance function
a(t) / b(t)
(t )
Fig. 2 Cell state model for elastodynamics 2.3.2 Fast updating rules for elastodynamics The updating rules are the most important part of the CDCA model, because the updating rules will determine the displacement, velocity, acceleration, strain and stress state of a cell element. Consider a cellular node N i for the plane stress problem, assume that the 9
u(t ) {u(t ), a(t ), b(t )}
displacement
,
velocity
u (t ) {u (t ), a (t ), b (t )}
and
(t )} are known at time t t , then at time (t ) {u (t ), a(t ), b acceleration u k
t1 tk t , the displacement of this node can be obtained due to the effect of nodal
, fib (t1 )} according to Eq. (12) and all degrees of the force vector fi (t1 ) {fiu (t1 ), fia (t1 )， nodal freedom are restricted on its neighbor cell nodes N ik , which can be shown in Fig. 3. The relation between the incremental force and incremental deformation can be
reflected
into
two
steps.
First,
the
nodal
force
increment
fi (t1 ) {fiu (t1 ), fia (t1 ), fib (t1 )} will lead the cell node N i to produce the displacement increment ui (t1 ) {ui (t1 ), ai (t1 ), bi (t1 )} . In order to accelerate the updating, an accelerate factor is defined as , and ui (t1 ) ui (t ) ui (t1 ) . Via
ui (t1 ) and the accelerate factor, one can obtain the velocity increment u (t1 ) {u (t1 ), a (t1 ), b (t1 )}
and
the
acceleration
increment
(t )} from time t t to time t t t . Then, the (t1 ) {u (t1 ), a(t1 ), b u 1 k k 1
i (t1 ) and the accelerate factor on the displacement increments u i (t1 ) , u i (t1 ) , u
cell node N i will lead its neighboring cell nodes to produce the nodal force k
N increment fi i (t1 ) .
10
N i3 neighbor
k
fiNi (t1) {fu (t1), fa (t1), fb (t1)}
N
4 i
neighbor
crack
u i (t1 ) u i (t1 ), a i (t1 ), b i (t1 )
ui (t1) ui (t1), ai (t1), bi (t1)
N i2 neighbor
cell
Ni
Kih k i k a / b
(t ) i (t1 ) u i (t1 ), ai (t1 ), b u i 1
Mih mi ma / b
f i (t1 ) {f u (t1 ), f a (t1 ), f b (t1 )}
N i1 neighbor Fig. 3 Local updating model for a cell and it neighbors Therefore, the process of the CDCA updating rules is: at time t1 tk t , the increment of nodal force leads to the increment of nodal displacement u i (t1 ) , nodal i (t1 ) , and the increments of nodal velocity u i (t1 ) and nodal acceleration u i (t1 ) lead displacement u i (t1 ) , nodal velocity u i (t1 ) and nodal acceleration u k
N to the increment of nodal force of its neighboring nodes fi i (t1 ) , until the system
dynamic equilibrium is achieved. In other words, the self-organization phenomenon k
N of ui (t1 ) 0 , u i (t1 ) 0 , u i (t1 ) 0 and fi i 0 appears. Therefore, the
updating steps for elastodynamics can be given as: (1) Obtain the initial value of displacement and velocity at time t 0 according to Eqs. (5) and (6); obtain the cell mass matrix and stiffness matrix M i , K i for cell Ni .
(2) The equilibrium equation of the cell node N i can be described as i (t1 ) K i ui (t1 ) fi (t1 ) Mi u 11
(20)
i (t1 ) , u i (t1 ) , and fi (t1 ) are increments of second time derivative of in which u
degrees of nodal freedom, degrees of nodal freedom and nodal force, respectively, of cell node N i . (3) Calculate the increment of degrees of nodal freedom u ih (t1 ) via the increment of nodal force fi (t1 ) and Eq. (20). (4) Via the increment of degrees of nodal freedom u ih (t1 ) and the accelerate factor
, one can obtain the displacement vector u i (t1 ) , which is ui (t1 ) ui (t ) ui (t1 )
(21)
(5) According to u ih (t1 ) and the accelerate factor , the velocity and the acceleration of cell N i can be easily obtained and are
u i (t1 ) u i (t ) u i (t1 )
(22)
i (t1 ) u i (t ) u i (t1 ) u
(23)
(6) Restrict all degrees of nodal freedom at all neighboring cells N ik , which can be seen in Fig. 3. k
N (7) Obtain the nodal force increment fi i (t1 ) of the neighboring cell N ik via
ui (t1 ) and the accelerate factor according to the following equation
i (t1 ) fiNi (t1 ) K iNi ui (t1 ) MiNi u k
k
k
k
(24)
k
N N where K i i and M i i are the stiffness matrix and mass matrix of neighboring cell
N ik , respectively. (8) Finish the calculation of steps (2) – (7) at all cell nodes, until ui (t1 ) 0 , k
N u i (t1 ) 0 , u i (t1 ) 0 and fi i 0 appear.
12
3. Time integration and dynamic stress intensity factors In the present method, an implicit time integration scheme of Newmark method is employed for the time derivative, and the interaction integral theory is developed in this section for solving the dynamic stress intensity factor. 3.1 Implicit time integration scheme The Newmark time integration scheme is widely used in dynamic analysis, because it can be designed to remain unconditionally stable. Taking cell N i as the example, and the time step from time t n to time tn1 tn t , then the time integration scheme can be given as i (tn1 ) K i ui (tn1 ) fi (tn1 ) Mi u
(25)
1 i (t n )t 2 u i (t n1 )t 2 u i (t n1 ) u i (t n ) u i (t n )t ( )u 2
(26)
i (tn )t u i (tn1 )t u i (tn1 ) u i (tn ) (1 )u
(27)
in which and are time difference constants; when 0.5 2 , this method is unconditionally stable, and when 0.5 and 2 , the method is stable if t 1 /(max / 2 ) . According to Eqs. (25) - (27), one can obtain all constants
that the present method uses:
0.5 0.25(0.5 ) 2 2 a0 1 /( t ) a3 1 /( 2 ) 1 a2 1 /( t ) a6 t (1 ) a7 t
(28)
In this case, the equivalent stiffness matrix can be written as ~ K i K i a0M i
(29)
Same as the equivalent stiffness matrix, the equivalent nodal loading vector can be written as ~ i (tn )) fi (tn1 ) fi (tn1 ) M i (a0ui (tn ) a2u i (tn ) a3u
According to Eqs. (29) - (30), one can obtain the equivalent equation 13
(30)
~ ~ K i ui (tn1 ) fi (tn1 )
(31)
Solving Eq. (31), one can obtain the displacement on time t n 1 , and based on the displacement on time step t n 1 , the velocity and the acceleration can be given as i (tn1 ) a0 (ui (tn1 ) ui (tn )) a2u i (tn ) a3u i (tn ) u
(32)
i (tn ) a7u i (tn1 ) u i (tn1 ) u i (tn ) a6u
(33)
Combining dynamic cellular automaton in section 2.3 and Eqs. (29) - (33), one can easily obtain the displacement, velocity and acceleration of cell N i . The same operation can be taken on all other cells, which can be seen that no assembled global stiffness matrix and mass matrix are needed in the whole calculation. Then, the memory is greatly saved, and by employing the fast updating scheme, the calculating efficiency can also be greatly improved. 3.2 Interaction integral for dynamic stress intensity factors In this section, a combination of interaction integral method and equivalent domain integral is employed to evaluate the DSIF[1,3,41]. In this method, the assumption of traction-free boundary condition on the crack surface is considered, and the dynamic equilibrium equation and strain displacement relation are employed. The J-integral [3] can be given as
J [ ij A
u j x1
W1i ]
2u u u 2ui q dA [( 2 j f b ) j i ]qdA A xi t x1 t x1t
(34)
in which q is a weighting function, and it is unity on the internal integral arc 0 and it is equal to zero on the outer arc 1 , which can be seen in Fig. 4, and in the present method, on nodes which are located inside the circle q 1 , on nodes which are located outside the circle q 0 , and on the other location of the integral element 4
q(x) N i (x)qi ; W ij ij , and ij , ij , u j are stress, strain and displacement i 1
on time tn1 tn t . Consider two states for this cracked body, which are State 1: 14
ij(1) , ij(1) , u (j1) are corresponding to the actual states, and State 2: ij( 2 ) , ij( 2) , u (j2) are the auxiliary states derived from Westergaard and Williams[42] for modes I and II cracks. Employing the interaction integral method and equivalent domain integral, the interaction integral for States 1 and 2 can be given as
M
(1, 2 )
[ A
(1) ij
u (j2) x1
( 2) ij
u (j1)
2u (j1) u (j2) q ] dA qdA A x1 xi t 2 x1 (1) ( 2 ) ij ij 1i
(35)
The integral domain can be seen in Fig. 4, which is marked by the shaded element, and h is the average area of elements. is the parameter of the integral domain, and it is usually given as 1 5 . According to the interaction integral method, the left hand term of Eq. (35) can be rewritten as M (1, 2)
2( K I(1) K I( 2) K II(1) K II( 2) ) ~ E
(36)
~ ~ in which E E for the plane stress problem, and E E /(1 2 ) for the plane strain problem. As we know, when K I( 2) 1.0 and K I(I2) 0.0 for State 2, K I can be obtained, and when K I( 2) 0.0 and K II( 2) 1.0 for State 2, K II can be obtained. Eq. (36) can be rewritten as M
(1, 2 )
(1) u (y2) u (y1) u x( 2) q (1) ( 2 ) u x ( 2) [ xy x xy x(1) x( 2) xy(1) xy( 2) ] dA A x x x x x ( 2 ) ( 1 ) u u q u ( 2) u (1) (37) [ xy(1) x y(1) y xy( 2 ) x y( 2) y ] dA A x x x x x 2 (1) ( 2) 2u x(1) u x( 2) u y u y [ 2 ]qdA A t x t 2 x (1) x
15
A
h
crack
R h
crack tip element
0
1 Fig. 4 Integral domain
4. Numerical results and discussions In this section, some numerical examples for elastodynamics problems under dynamic loads with the assumption of plane strain condition of 2D elasticity are considered. Time differential parameters 0.25 and 0.5 are used in those examples. In addition, many results by some other methods are introduced in this paper for comparison. 4.1 Mode-I semi-infinite crack with Heaviside step loading An infinite plate with a semi-infinite mode-I crack is considered in this section, which is loaded by a tensile stress perpendicular to the crack surface. The analytical solution of DSIF for the present example by Freund[43] is employed to validate the accuracy and efficiency of the present method. The geometry of the model is as follows: plate length L 10m , crack length a 5m , and vertical half height of the plate H 2m . The geometry and boundary conditions can be seen in Fig. 5. The material properties are: Young's modulus E 210GPa , Poisson's ratio 0.3 , and density of material 8000kg / m3 . A tensile Heaviside step stress of 0 500MPa is applied on the top surface of the plate, which can be seen in Fig. 6. According to references [3, 20, 44], the total time for the simulation is limited to
t 3tc 1.009 103 s , which is the time for when the tensile stress wave is reflected 16
on the bottom side and reaches again to the crack tip, and tc H / cd , where cd is the dilatational wave speed. For convenience, normalized DSIF K I (t ) /( 0 H ) is used in this example, and the analytical equation can be given as
K I (t ) 2 0 1
0 cd (t t c )(1 2 )
t tc t tc
(38)
0
2H=4m a=5m
L=10m
Fig. 5 Geometry and boundary conditions of a semi-infinite crack In the present method, t 1.2 105 , and the number of time steps is 90. The present method results of mode-I DSIF compared with the results by sES-FEM and analytical solutions are plotted in Fig. 7, in which one can see that the results by the present method agree with the analytical solutions and results by sES-FEM[3], and high accuracy and efficiency can be achieved by the present method.
0
t Fig. 6 Heaviside step loading
17
1.5 CDCA Analytical sES-FEM[3]
1.2
KI(t) / σ0H1/2
0.9 0.6 0.3 0.0 -0.3 0.0
0.3
0.6
0.9
1.2
1.5
1.8
2.1
2.4
2.7
3.0
t / tc
Fig. 7 Normalized mode-I DSIF results by different methods 1.5 101 83 63 41 21 Analytical
KI(t) / σ0H1/2
1.2 0.9 0.6 0.3 0.0 -0.3 0.0
0.3
0.6
0.9
1.2
1.5
1.8
2.1
2.4
2.7
3.0
t / tc
Fig. 8 Normalized mode-I DSIF results by different meshing elements Results from different meshing densities with the present method are shown in Fig. 8, in which the results are much closer to the analytical solutions with the increase in the node numbers. Fig. 9 plots the convergence of the present method with different meshing densities, and it is shown that convergence can be achieved with the increase in the node numbers, and the relative errors become much smaller with much more node numbers.
18
0.5 10141 8333 6325 4117 219
||KI(t)-KIA(t)|| / σ0H1/2
0.4 0.3 0.2 0.1 0.0
-0.1 0.0
0.3
0.6
0.9
1.2
1.5
1.8
2.1
2.4
2.7
3.0
t / tc
Fig. 9 Convergence of DSIF by different meshing elements 400
CDCA XFEM
350
Memory (MB)
300 250 200 150 100 50 0 0
3000
6000
9000
12000
15000
Node numbers
Fig. 10 Memory expenses comparison between the present method and XFEM 1050
CDCA XFEM
900
CPU time (s)
750 600 450 300 150 0 0
3000
6000
9000
12000
15000
Node numbers
Fig. 11 CPU time comparison between the present method and XFEM
Fig. 10 shows the memory expenses by the present method with different node numbers, in which it can be seen that the calculating memories determined by the present method are much smaller than those determined by the XFEM because no assembled global stiffness and mass matrix are needed in the whole calculation, and 19
in XFEM, the recent half bandwidth technique is employed. Fig. 11 plots the CPU time expending by the present method with different node numbers, and one can see that the CPU time by the present method is close to XFEM with much smaller node numbers. However, the CPU time of the present method is much smaller than that of XFEM with much larger node numbers, and the CPU time increases much more gently with the increase in the node numbers. Thus, the advantage is very obvious in the present method for much larger node numbers.
5.2 A slanted edge crack plate with Heaviside step loading To illustrate the efficiency and reliability of the present time domain CDCA for simulating mixed mode dynamic crack problems, a slanted edge crack is considered in this section. A rectangular plate with an edge crack slanted at an angle 45 is employed, which is plotted in Fig. 12. The plate is supported with left, upper and bottom sides on normal directions, and the right side is subjected to a uniform tensile Heaviside step loading 0 , which can be seen in Fig. 6. The geometrical parameters of the plate are: b 44mm , h 32mm , c 6mm , and the crack length
a 22.63mm . The material properties can be given as: shear modulus G 29.4GPa , Poisson ratio 0.286 , and the density of the material 2.45 10 6 kg / mm3 [1,3]. For comparison with the reference [1,3], the total time of the simulation is t 25s . In addition, the normalized DSIF K I,II / 0 a is considered.
20
b
0
h a
c Fig. 12 Geometry and boundary conditions of a slanted edge crack plate 1.8 MFGM[28] CDCA DBEM[46] sES-FEM[3]
1.5
KI
KI,II(t) / σ0(πa)1/2
1.2 0.9 0.6
KII
0.3 0.0 0
5
10
15
20
25
Time t (μs)
Fig. 13 Normalized mixed modes DSIF results by different meshing elements 1.8
α=4.8 α=4.4 α=4.0 α=3.6 α=3.2 α=2.8 α=2.4 α=2.0 α=1.6
1.5
KI(t) / σ0(πa)1/2
1.2 0.9 0.6 0.3 0.0 0
5
10
15
20
25
Time t (μs)
Fig. 14 Normalized mode-I DSIF results by different integral parameters
21
1.2
α=4.8 α=4.4 α=4.0 α=3.6 α=3.2 α=2.8 α=2.4 α=2.0 α=1.6
KII(t) / σ0(πa)1/2
0.9
0.6
0.3
0.0 0
5
10
15
20
25
Time t (μs)
Fig. 15 Normalized mode-II DSIF results by different integral parameters In the present method, t 5.0 107 is used, and time step number is 50, and a total of 5840 nodes and 5785 elements are used in this example. A comparison of the normalized mode-I and mode-II DSIFs by different methods is plotted in Fig. 13, in which one can see that the present method results agree well with the mesh-free Galerkin method (MFGM) by Wen and Aliabadi[28], dual boundary element method (DBEM) by Fedelinski et al. [45] and sES-FEM by Liu et al. [3], and it is verified that a high accuracy of the present method can be achieved. In this example, the influence of integral domain parameter on the equivalent domain integral of DSIF is studied, which can be seen in Figs. 14 - 15, and it can be seen in those figures that a suitable value of the parameter cannot influence the accuracy of the DSIF of the present method, and for 1.6 4.8 , the results are almost the same with each other. Then, the conclusion is found that 1.6 4.8 should be used in order to preserve the accuracy of the present method.
5.3 An oriented central crack plate with Heaviside and blast loading In this example, a rectangular plate with an arbitrarily oriented central crack is considered, and the geometry of the model is shown in Fig. 16. The height of the model is 2H 40mm , the width of the model is 2W 20mm , and the crack length is
2a 4.8mm . Two types of loadings are considered in this section, which are Heaviside step loading and blast loading, and the latter can be seen in Fig. 17. Those loadings are uniform and subjected on the top and bottom sides of the model. The 22
material properties are: Young's modulus E 200GPa , Poisson's ratio 0.3 , and material density 5000kg / m3 . For comparison, the results from the sES-FEM[3] and solutions from the traditional XFEM are shown in the same figures as the present method results, and the total simulating time is t 20s . In addition, different angles of crack are considered, i.e., 15o ,30o ,45o ,60o ,75o . For convenience, the normalized DSIF K I, II / 0 a is considered. In the present method, t 3.0 107 is used, and time step number is 68. A total of 7616 nodes and 7437 elements are used in this example, and the integral parameter 3.2 for the equivalent domain integral is used in the present example.
0
y
2H
2a
x 2W
0 Fig. 16 Geometry and loading of a central crack plate
23
0
t t2
t1
Fig. 17 The blast loading 2.8 2.4
CDCA XFEM sES-FEM[3]
15
2.0
30
KI(t) / 0a)1/2
1.6 45
1.2 0.8 0.4 0.0
60
-0.4
75 0
4
8
12
16
20
Time t (s)
Fig. 18 Comparison of the normalized mode-I DSIF by XFEM, Ses-FEM and the present method for different inclination angles and Heaviside step loading 1.6
CDCA XFEM sES-FEM[3]
45
1.2
KII(t) / 0a1/2
30&60 0.8
0.4 15&75 0.0
-0.4 0
4
8
12
16
20
Time t (s)
Fig. 19 Comparison of the normalized mode-II DSIF by XFEM, sES-FEM and the present method for different inclination angles and Heaviside step loading
24
1.8
15
CDCA sES-FEM[3]
30 1.2 45
KII(t) / 0a)1/2
0.6 0.0
60 75
-0.6 -1.2 -1.8 0
4
8
12
16
20
Time t (s)
Fig. 20 Comparison of the normalized mode-I DSIF by XFEM, Ses-FEM and the present method for different inclination angles and blast loading 1.5
45
CDCA sES-FEM[3]
30
1.2 0.9
KII(t) / 0(a)1/2
0.6 0.3 15
0.0 -0.3 -0.6 -0.9 -1.2 0
4
8
12
16
20
Time t (s)
Fig. 21 Comparison of the normalized mode-II DSIF by XFEM, Ses-FEM and the present method for different inclination angles and blast loading In this example, two types of dynamic loadings are considered. First, the Heaviside step loading is employed, and the normalized mode-I and mode-II DSIFs for the five different inclination angles are presented in Fig. 18 and Fig. 19, respectively. For comparison, the results by sES-FEM[3] and XFEM are employed in those figures, and the results given in those figures generally show good agreement among the different methods utilized. It can be seen that the amplitude of normalized mode-I DSIF increases with the inclination angle, but the maximum amplitude of the normalized mode-II DSIF occurs at
. Fig. 19 reveals that the normalized
mode-II DSIF is almost the same for
and
and for
and
. In addition, a blast dynamic loading is also considered in Fig. 17, and the loading parameters are given as:
,
. The normalized mode-I and mode-II 25
DSIF for blast loading can be seen in Fig. 20 and Fig. 21, respectively, which show that the present method results are close to those by XFEM and sES-FEM[3], and the same phenomenon for different inclination angles can be obtained for the blast loading, in which the maximum amplitude of normalized mode-I DSIF is at and the minimum one is at mode-II DSIF is at
,
, and the maximum amplitude of normalized
.
5. Conclusions In this paper, a time domain continuous-discontinuous cellular automaton method has been proposed for elastodynamic problems and is based on the static continuous-discontinuous cellular automaton method. First, the continuity to discontinuity enriched function technique is developed for dynamic fracture problems, in order to ensure the interpolation property is satisfied. The shifting amendments are employed in the present method. Second, a time domain cellular automaton theory for elastodynamics has been proposed, in which new cell states for dynamic fracture problems, dynamic cell model and dynamic cellular updating scheme are developed. Besides, fast adaptive updating rules based on time domain cellular automaton for the cell and its neighbors are developed. Third, several time-dependent equations for cells are discretized individually by the Newmark method, and rapidly solved by the fast cellular automaton updating method. Furthermore, to analyze the fracture behaviors of elastic solids with dynamic loads, mixed-mode dynamic stress intensity factors are obtained by the interaction integral method and equivalent domain integral. Because of the usages of the discontinuous enriched function technique, time domain cellular automaton theory and fast adaptive updating rules, the cracks are independent of the integral and meshing calculations, and no remeshing is needed for crack dynamic propagation. In addition, no assembled global matrix is needed for the whole calculating process, so the required computing memory can be greatly reduced. Furthermore, the calculation is located at the cell locality, so the local properties can be easily considered. Additionally, the calculating efficiency can be greatly improved 26
via the fast adaptive updating scheme. Finally, some numerical examples for dynamic stress intensity factor problems are studied in this paper, which validate that both the accuracy and efficiency are greatly improved by using the present method. Solutions determined by the analytical method and some other methods are employed to show that a high accuracy can be achieved by the time domain CDCA, and the examples show that the suitable integral domain parameter is between 1.6 and 4.8.
Acknowledgments The work was financially supported by the State Key Research Development Program of China (No. 2017YFC0804203), the International Partnership Program of the Chinese Academy of Sciences (115242KYSB20160017)and the Youth Innovation Promotion Association of CAS(No. 2014304).
References [1] K. Kishimoto, S. Aoki, M. Sakata. Dynamic stress intensity factors using J-integral and finite element method. Engrg. Fract. Mech. 1980, 13: 387-394 [2] J.W. Kelley, C.T. Sun. A singular finite element for computing time dependent stress intensity factors. Engrg. Fract. Mech. 1979, 12: 13-22 [3] P. Liu, T.Q. Bui, C. Zhang, T.T. Yu, G.R. Liu, M.V. Golub. The singular edge-based smoothed finite element method for stationary dynamic crack problems in 2D elastic solids. Comput. Methods Appl. Mech. Engrg. 2012, 233-236: 68-80 [4] J. Rethore, A. Graviouil, A, Combescure. A stable numerical scheme for the finite element simulation of dynamic crack propagation with remeshing. Comput. Meth. Appl. Mech. Engrg. 2004, 193: 4493-4510 [5] F. Chirino, J. Dominguez J. Dynamic analysis of cracks using boundary element method. EngFract Mech. 1989, 34: 1051-1061 [6] D.P.N Kontoni, D.E. Beskos. Boundary element formulation for dynamic analysis of nonlinear systems. Engng Anal. 1988, 5: 114-125 [7] T. Burczynski, T. Adamczyk. Analysis of nonlinear systems in terms of the 27
boundary element method. Mech. Comput. 1988, 7: 149-164 [8]
S.
Ahmad,
P.K.
Banerjee.
Inelastic
transient
dynamic
analysis
of
three-dimensional problems by BEM. Int. J. Num. Meth. Engng. 1990, 29: 371-390 [9] D. P. N. Kontoni, D. E. Beskos. Inelastic dynamic analysis by the boundary element method. In Boundary Elements IX, Proceedings of the 9th International conference on Boundary Element Methods in Engineering(Vol. 2), ed. C.A. Brebbia, W.L. Wendland& G. Kuhn. Springer-Verlag, Berlin, Germany, 1987, 335-351 [10] T. Panzeca, C. Polizzotto, M. Zito. A BEM formulation of the dynamic elastic-plastic structural problem via variational principle, In Boundary Elements in Machanical and Electrical Engineering, Proceedings of the International Boundary Element Symposium, ed. C.A. Brebbia& A. Chaudouet-Miranda, Computational Mechanics Publications, Southanpton, UK and Springer-Verlag, Berlin, Germany, 1990, 193-204 [11] J.A.M Carrer, J.C.F. Telles. Transient dynamic elastoplastic analysis by the boundary element method. In Boundary Element Technology VI, Proceedings of the 6th Int. Conf. on Boundry Element Technology, ed. C.A. Brebbia. CMP, Southanpton, UK, and Elsevier, London, UK, 1991, 265-277 [12] J.P. Agnantiaris, D. Polyzos, D.E. Beskos. Some studies on dual reciprocity BEM for elastodynamic analysis. Comput. Mech. 1996, 17: 270-277 [13] J.P. Agnantiaris, D. Polyzos, D.E. Beskos. Free vibration analysis of non-axisymmetric and axisymmetric structures by the dual reciprocity BEM. Eng. Anal. Bound. Elem. 2001, 25: 713-723 [14] M. F. Samaan, Rashed Y. F. BEM for transient 2D elastodynamicsusing multiquadric functions. Int. J. Solids Struct. 2007, 44: 8517-8531 [15] J. Dominguez, R. Gallego. Time domain boundary element method for dynamic stress intensity factor computations. Int. J. Num. Meth. Eng. 1992, 33: 635-647 [16] H. Gao, P. Klein. Numerical simulation of crack growth in an isotropic solid with randomized internal cohesive bonds. J. Mech. Physics Solids. 1998, 42: 187-218 [17] T. Belytschko, H. Chen, J.X. Xu, G. Zi. Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment. Int. J. Numer. Meth. Eng. 28
2003, 58: 1873-1905 [18] J. Chessa, T. Belytschko. Arbitrary discontinuities in space-time finite element by level sets and X-FEM. Int. J. Numer. Meth. Engrg. 2004, 61: 2595-2614 [19] T. Menouillard, J. Rethore, A. Combescure, H. Bung. Efficient explicit time stepping for the extended finite element method(XFEM). Int. J. Numer. Meth. Engrg. 2006, 68: 911-939 [20] T. Elguedj, A. Gravouil, H. Maigre. An explicit dynamics extended finite element method. Part 1: Mass lumping for arbitrary enrichment functions. Comput. Meth. Appl. Mech. Eng. 2009, 198: 2297-2317 [21] J. Rethore, A. Gravouil, A. Combescure. An energy-conserving scheme for dynamic crack growth using the extended finite element method. Int. J. Numer. Meth. Engng. 2005, 63: 631-659 [22] L. Wu, P. Liu, C. Shi, Z. Zhang, T. Q. Bui, D. Jiao. Edge-based smoothed extended finite element method for dynamic fracture analysis. Appl. Math. Model. 2016, 40: 8564-8579 [23] C. A. Durate, O. N. Hamzeh, T. J. Liszka, W. W. Tworzydlo. A generalized finite element method for the simulation of three-dimensional dynamic crack propagation. Comput. Meth. Appl. Mech. Eng. 2001, 190: 2227-2262 [24] W. Li, H. Zheng, G. H. Sun. The moving least squares based numerical manifold method for vibration and impact analysis of cracked bodies. Eng. Fract. Mech. 2018, 190: 410-434 [25] T. Belytschko, M. Tabbara. Dynamic fracture using element-free Galerkin methods. Int. J. Numer. Meth. Eng. 1996, 39: 923-938 [26] P. Krysl, T. Belytschko. Then Element free Galerkin method for dynamic propagation of arbitrary 3-D cracks. Int. J. Numer. Meth. Engng. 1999, 44: 767-800 [27] Y. Lu, T. Belytschko, M. Tabbara. Element-free Galerkin method for wave propagation and dynamic fracture. Comput. Meth. Appl. Mech. Eng. 1995, 126: 131-153 [28] P. H. Wen, M. H. Aliabadi. Evaluation of mixed-mode stress intensity factors by mesh-free Galerkin method: static and dynamic. J. Strain Anal. Eng. 2009, 44: 29
273-286 [29] T. Rabczuk, J.H. Song, T. Belytschko. Simulations of instability in dynamic fracture by the cracking particles method. Eng. Fract. Mech. 2009, 76: 730-741 [30] F. Yan, Y.H. Wang, L.G. Tham, Y.K. Cheung. Dual reciprocity hybrid boundary node method for 2-D elasticity with body force. Eng. Anal. Bound. Elem. 2008, 32: 713-725 [31] F. Yan, Y. Miao, Q.N. Yang. Quasilinear hybrid boundary node method for solving nonlinear problems. CMES-Comput. Model. Eng. Sci. 2009, 46: 21-50 [32] F. Yan, X. T. Feng, H. Zhou. A dual reciprocity hybrid radial boundary node method based on radial point interpolation method. Comput. Mech. 2010, 45: 541-552 [33] F. Yan, X. T. Feng, J. H.Lv, P. Z. Pan. A new hybrid boundary node method based on Taylor expansion and Shepard interpolation method. Int. J.Numer. Meth. Eng. 2015, 102: 1488-1506 [34] F. Yan, Y.H. Wang, Y. Miao, L.G. Tham, Y.K. Cheung. Dual reciprocity hybrid boundary node method for free vibration analysis. J. Sound Vib. 2009, 321: 1036-1057 [35] Y. Miao, T. G. He, H Luo, H. P. Zhu. Dual hybrid boundary node method for solving transient dynamic fracture problems. CMES-Comput. Model. Eng. Sci. 2012, 85: 484-498 [36] X.T. Feng, P.Z. Pan, H. Zhou. Simulation of the crack microfracturing process under uniaxial compression using an elasto-plastic cellular automaton. Int. J. Rock Mech. Min. Sci. 2006, 43: 1091-1108. [37] F. Yan, X.T. Feng, P.Z. Pan, S.J. Li. Discontinuous cellular automaton method for crack growth analysis without remeshing. Appl. Math. Model. 2014, 38: 291-307 [38] F. Yan, X. T. Feng, P. Z. Pan, S. J. Li. A continuous-discontinuous cellular automaton method for cracks growth and coalescence in brittle material. Acta Mech.Sinica. 2014, 30: 73-83 [39] F. Yan, X. T. Feng, J. H.Lv, P. Z. Pan, S. J. Li. Continuous-discontinuous cellular automaton method for cohesive crack growth in rock. Eng.Fract. Mech. 2018, 188: 361-380 30
[40] F. Yan, X.T. Feng, P.Z. Pan, Y.P. Li. A continuous-discontinuous cellular automaton method for regular frictional contact problems. Arch. Appl. Mech. 2013, 83: 1239-1255. [41] S.H. Song, G. H. Paulino. Dynamic tress intensity factors for homogeneous and smoothly heterogeneous materials using the interaction integral method. Int. J. Solids Struct. 2006, 43: 4830-4866 [42] N. Moes, J. Dolbow, T. Belytschko. A finite element for crack growth without remeshing. Int. J. Numer. Mehods Engrg. 1999, 46: 131-150 [43] L.B. Freund. Dynamic fracture mechanics. Cambridge University Press, London, 1990 [44] T. Menoulillard, T. Belytschko. Dynamic fracture with meshfree enriched XFEM. Acta Mech. 2010, 213: 53-69 [45] P. Fedelinshi, M.H. Aliabadi, D.P. Rooke. The laplace transforms DBEM for mixed-mode dynamic crack analysis. Comput. struct. 1996, 59: 1021-1031
31
Highlights: A continuity-to-discontinuity enriched function technique is developed for elastodynamics. Fast adaptive cellular automaton updating rules for the cell and its neighbors are developed. A time domain cellular automaton theory for elastodynamics is proposed. The cracks are independent of the integral meshes and no global matrix is needed.
32
S0013-7944(18)30758-6 https://doi.org/10.1016/j.engfracmech.2018.10.025 EFM 6199
To appear in:
Engineering Fracture Mechanics
Received Date: Revised Date: Accepted Date:
28 July 2018 11 October 2018 24 October 2018
Please cite this article as: Yan, F., Pan, P-Z., Feng, X-T., Li, S-J., The continuous-discontinuous cellular automaton method for elastodynamic problems, Engineering Fracture Mechanics (2018), doi: https://doi.org/10.1016/ j.engfracmech.2018.10.025
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.
The continuous-discontinuous cellular automaton method for elastodynamic problems Fei Yan1, Peng-Zhi Pan1,2,3*, Xia-Ting Feng3, Shao-Jun Li1 1. State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China 2. University of Chinese Academy of Sciences, Beijing 100049, China 3. Northeastern University, Shenyang 110819, Liaoning, China
Abstract: In this paper, the recently proposed continuous-discontinuous cellular automaton method (CDCA) is further developed for elastodynamic problems in 2D elastic solids. A continuity-to-discontinuity enriched function technique is developed for elastodynamics and combined with a cellular automaton formulation; the cellular automaton theory for elastodynamics, in which fast adaptive updating rules for the cell and its neighbors are developed, and the time domain CDCA are proposed. In this method, the cracks are independent of the integral meshes, and no assembled global matrix is needed for the whole process. Time-dependent equations for the cells are discretized by the Newmark method, and rapidly solved by the cellular automaton method. To analyze the fracture behaviors of elastic solids with dynamic loads, mixed-mode dynamic stress intensity factors are obtained by the interaction integral method. Some numerical examples for elastodynamic problems are studied to validate the accuracy and efficiency of the present method, and solutions from analytical method and some other methods are employed to show the high accuracy of the CDCA.
Keywords: the time domain continuous-discontinuous cellular automaton; cellular automaton method; fast updating scheme; dynamic stress intensity factor; elastodynamic problem
*
Corresponding author. Email address: [email protected] (P.- Z. Pan) 1
1. Introduction The modeling of dynamic crack problems is of great importance in the simulation of structure failure. Because of the mathematical difficulties of time-dependent loading, analytical methods are often not feasible to solve general dynamic crack problems in practical engineering, and numerical methods are thus developed for solving such problems. However, accurate and highly efficient numerical simulations of dynamic fracture problems remain a challenging task for many practical applications in engineering, and many significant numerical methods have been developed for dynamic crack problems in the past few years. As the most widely used simulation method, the finite element method (FEM) is also widely used to solve dynamic fracture problems[1-2], in which a singular element is developed. Based on cell-based smoothing domains, a singular edge-based smoothed finite element method (sES-FEM) is proposed for stationary dynamic crack problems in 2D elastic solids[3]. FEM, with remeshing and a stable numerical scheme, has also been used for dynamic crack propagation[4], but there are some disadvantages of the FEM such as the need for discretizing the entire domain; this is especially disadvantageous when remeshing for crack propagation, and inaccuracies arise in the form of stress concentrations. As a boundary-type numerical method, the boundary element method (BEM) has been widely applied in dynamic crack problems, such as the dynamic analysis of crack by Chirino and Dominguez[5], dynamic analysis of nonlinear systems by Konroni and Beskos[6] and Burczynski and Adamczyk[7], inelastic dynamic analysis by Ahmad and Banerjee[8] and Kontoni and Bekos[9], and dynamic elastic-plastic structural problem by Panzeca et al. [10] and Carrer and Telles[11]. To overcome the required domain integrations, dual reciprocity method (DRM) is employed in BEM for solving the elasto-dynamics, such as transient elastodynamic and free vibration analysis by Agnantiaris et al. [12-13] and Samaan and Rashed[14] and dynamic stress intensity factor (DSIF) computations by Dominguez and Gallego[15]. Recently, some techniques have been developed to avoid remeshing and are 2
suitable for crack propagation, i.e., the extended finite element method (XFEM). Gao and Klein[16] first developed the idea of loss of hyperbolicity for analyzing dynamic crack propagations. Belytschko et al.[17] developed a new discontinuous enrichment for dynamic crack propagation based on the loss of hyperbolicity. Further, Chessa and Belytschko[18] presented a locally enriched space-time XFEM for solving hyperbolic problems. Menouillard et al.[19] introduced a lumped mass matrix for enriched element to allow the use of pure explicit formulation in XFEM. Elguedj et al.[20] proposed a mass lumping for arbitrary enrichment functions. Rethore et al.[21] employed an energy-conserving scheme for dynamic crack growth, and Wu et al.[22] developed an edge-based smoothed extended finite element method for dynamic fracture analysis. Duarte et al.[23] proposed a generalized finite element method (GFEM) for 3D dynamic crack propagation, and Li et al.[24] combined moving least squares and numerical manifold method (NMM) for impact analysis of cracked bodies. In recent years, meshless methods have also been developed for dynamic crack problems due to their non-meshing, such as element-free Galerkin method (EFG), dynamic fracture and wave propagation analysis by Belytschko and Tabbara[25] and further for 3D dynamic propagation by Krysl and Belytschko[26] and Lu et al. [27], mixed-mode stress intensity factors evaluation by Wen and Aliabadi[28], the cracking particles method (CPM) employed by Rabczuk et al. [29] to simulate dynamic fracture, the hybrid boundary node method (HBNM) and dual reciprocity hybrid boundary node method[30-33], free vibration by Yan et al. [34], and transient dynamic fracture problem by Miao et al. [35]. To overcome the defect of remeshing and to easily consider the plastic fracture propagation, based on the cellular automaton method[36] and discontinuous enrichment functions, Yan et al.[37] proposed a continuous-discontinuous cellular automaton method (CDCA), and employed this method to solve the static multiple cracks propagation[38], cohesive crack propagation[39], contact crack[40] and so on. The advantages of this method are: first, the cracks are independent of the calculation mesh, so no remeshing is needed with the cracks propagation; second, no assembled 3
global matrix is needed in the whole calculation, the computing memory can be greatly saved, and the calculation is located at the local cell, so the local properties can be easily considered; finally, with the fast updating method, the calculating efficiency can be greatly improved and can be easily applied for some large practical engineering. Based on the continuous-discontinuous cellular automaton method for static fracture problems, a time domain CDCA is proposed in this paper. First, the continuity to discontinuity enriched function technique is developed for dynamic fracture problems; second, a time domain cellular automaton theory for elastodynamics is proposed, in which fast adaptive updating rules based on time domain cellular automaton for the cell and its neighbors are developed; then, time-dependent equations for the cells are discretized individually by the Newmark method, and rapidly solved by the cellular automaton method; finally, to analyze the fracture behaviors of elastic solids with dynamic loads, mixed-mode dynamic stress intensity factors are obtained by the interaction integral method. Because of the usage of time domain cellular automaton theory and fast adaptive updating rules, the cracks are independent of the integral meshes, and no assembled global matrix is needed for the whole process. The computing memory can be greatly saved, and the calculation is located at the local cell, so the local properties can be easily considered. In addition, the calculating efficiency can be greatly improved. Finally, some numerical examples of elastodynamic problems are studied to validate the accuracy and efficiency of the present method, and solutions from analytical method and some other methods are employed to show the high accuracy of the CDCA.
2. The continuous-discontinuous cellular automaton 2.1 Governing equation of elastodynamics Consider a body
with an initial crack
in the state of dynamic equilibrium,
the governing elastodynamic equation can be given as σ f b u 4
(1)
in which σ is the stress tensor, f b is the body force, is the density of the
is the acceleration vector, which is on the following boundary material, and u conditions: u(x, t ) u (x, t )
on u
(2)
σ n f t
on t
(3)
σ n 0
on c
(4)
in which u , t and c are the displacement, traction and crack boundaries, respectively, and f t is the external traction vectors. The initial conditions can be written as
u(x, t 0) u (0)
(5)
u (x, t 0) u (0)
(6)
in which u (0) and u (0) are the initial value vector of displacement and velocity, respectively. According to variational theory, the variational formulation of this boundary value problem of Eq. (1) can be given as
ud σ εd f b ud f t ud u
(7)
2.2 Continuous-discontinuous theory for cellular automaton Based on static continuous-discontinuous cellular automaton theory, the Heaviside function and the exact near-tip asymptotic field functions are employed in this method to model the strong discontinuity and high gradient stress field around the crack tip. Then, the shape function at time t of the present method can be given as[37-39] n
m
u h (x, t ) N j (x)u j (t ) N k (x)H ( ) H ( k ) a k (t ) j 1 k 1 kP
N i (x) Fl (x) Fl (x i ) b li (t ) i 1 l 1 t
nf
iT
5
(8)
in which N j (x) is the traditional FEM shape function, u j (t ) is the nodal displacements at time t , a k (t ) is a vector of additional degree of nodal freedom for modeling strong discontinuity t , and b li (t ) is a vector of additional degrees of nodal freedom for modeling the crack tip stress field at time t . The signed function is chosen as the Heaviside enrichment function, which is given as
1 0 H ( ) sign ( ) 1 0
(9)
in which is the value of level set function, which can be seen in[37-39] and Fl (x) is the exact near-tip asymptotic field functions, which can be given as[37-39]
{Fl ( x), l 1 4} r sin( ), r cos( ), r sin( ) sin( ), r sin( ) cos( ) 2 2 2 2
(10)
in which r , are the polar coordinate of the crack tip polar coordinate system. To ensure the interpolation property is satisfied, the shifting amendments of H ( ) H ( k ) and Fl (x) Fl (xi ) are employed in Eq. (8). Following Eq. (8), the
acceleration of nodes can be given as n
m
h (x, t ) N j (x)u j (t ) N k (x)H ( ) H ( k ) ak (t ) u j 1 k 1 kP
l (t ) N i (x) Fl (x) Fl (x i ) b i i 1 l 1 t
nf
(11)
iT
Then, the discretized form of Eq. (7) using Eqs. (8) and (11) can be rewritten as
(t ) Ku(t ) F(t ) Mu
(12)
(t ) denote the vector of nodal parameters and its second time in which u(t ) and u derivative at time t , respectively, which include displacement u j (t ) , enrichment
j (t ) and second time derivative of a k (t ) . degrees of freedom a k (t ) , acceleration u The stiffness matrix K , mass matrix M and the external loading vector F 6
can be defined as K uu K ij [ ijau K ij M uu M ij [ ijau M ij
K ua ij ] K ijaa
(13)
M ua ij ] M ijaa
(14)
Fi (t ) Fiu (t ) Fia (t )
T
(15)
in which the stiffness matrix K and the external loading vector F can be seen in references[37-39], and the mass matrix can be given as Muu ij N i N j d
M ua ij N i N j [ H ( ) H ( j )]d
M ijaa Ni [ H ( ) H (i )]N j [ H ( ) H ( j )]d
(16) (17) (18)
in which N i is the traditional FEM shape function for node i . If the enrichment is the exact near-tip asymptotic field functions, Eq. (18) can be rewritten as M bb ij N i [ Fl (x) Fl (x i )]N j [ Fl (x) Fl (x j )]d
(19)
2.3 Cellular automaton theory for dynamic problem By using discontinuous cellular automaton, the equilibrium state of the object can be obtained through the self-organization phenomenon formed by the one-another transfer of information between cells. Its advantages are: first, the localization property of the calculating structure can be easily considered, and the behavior of the cell is thought to be essentially local; in other words, the state of one cell is only determined by the states of itself and its neighbors. Second, there is no need to assemble the overall matrix, and the different degrees may bring some difficulties for the assembling operation for the enriched cells. Finally, the large-scale simulation can be performed because of its easy implementation property of the parallel algorithm. 7
The CDCA model is composed of cell, cell space, cell state, crack, neighborhood, updating rules and so on, and the relation between those components can be seen in references[37-39]. Crack, cell, neighbors and their relations can be seen in Fig. 1.
strong discontinuous cell
neighborhood
crack tip cell
crack tip element through element
cell
crack continuous element
neighborhood
continuous cell
Fig. 1 Cell, neighbors, crack and their relations 2.3.1 Cell space and its states for elastodynamics Same as the static cellular automaton model, the 2D cell space can be rectangular, triangular, hexagon and so on. To describe the cell states, a series of physical and mechanical properties must be defined to determine its different states. A CDCA model
is
composed
of
the
degree
values
vector
of
nodal
freedom
u(t ) {u(t ), a(t ), b(t )} , in which u(t ) is the traditional degree of nodal freedom, a(t ) is the Heaviside enriched degree of nodal freedom and b(t ) is the crack tip field
function
enriched
degree
of
nodal
freedom,
time
derivative
u (t ) {u (t ), a (t ), b (t )} , which includes the velocity of the node and time derivative of
the enriched nodal freedom, a(t ) or b(t ) , and the second time derivative (t )} , which includes the acceleration of node and second time (t ) {u (t ), a(t ), b u
derivative of enriched nodal freedom a(t ) or b(t ) . 8
In addition, the cell state includes current time t t k , time step increment t , the material property of thickness h , Young’s modulus E , Poisson’s ratio and fracture toughness K IC ; the cell force, strain and stress state can be given in the cell nodal force vector f (t ) {fu (t ), fa (t ), fb (t )} , in which the subscripts u , a and b represent the traditional, Heaviside enriched and exact near-tip asymptotic field functions enriched by the degrees of nodal freedom, respectively; additionally, the cell state includes the cell strain ε , cell stress σ , and crack location information, i.e., the distance function , and so on. Based on these theories, the cell state can be given in Fig. 2.
displacement
u (t )
velocity u (t )
strain
(t )
stress (t )
time derivative of enriched freedom
a(t ) / b(t )
cell state
enriched nodal freedom
t tk
enrihed nodal force
a(t) / b(t)
Fa / b (t) traditional nodal force
acceleration
Fu (t)
u(t )
second time derivative of enriched freedom
distance function
a(t) / b(t)
(t )
Fig. 2 Cell state model for elastodynamics 2.3.2 Fast updating rules for elastodynamics The updating rules are the most important part of the CDCA model, because the updating rules will determine the displacement, velocity, acceleration, strain and stress state of a cell element. Consider a cellular node N i for the plane stress problem, assume that the 9
u(t ) {u(t ), a(t ), b(t )}
displacement
,
velocity
u (t ) {u (t ), a (t ), b (t )}
and
(t )} are known at time t t , then at time (t ) {u (t ), a(t ), b acceleration u k
t1 tk t , the displacement of this node can be obtained due to the effect of nodal
, fib (t1 )} according to Eq. (12) and all degrees of the force vector fi (t1 ) {fiu (t1 ), fia (t1 )， nodal freedom are restricted on its neighbor cell nodes N ik , which can be shown in Fig. 3. The relation between the incremental force and incremental deformation can be
reflected
into
two
steps.
First,
the
nodal
force
increment
fi (t1 ) {fiu (t1 ), fia (t1 ), fib (t1 )} will lead the cell node N i to produce the displacement increment ui (t1 ) {ui (t1 ), ai (t1 ), bi (t1 )} . In order to accelerate the updating, an accelerate factor is defined as , and ui (t1 ) ui (t ) ui (t1 ) . Via
ui (t1 ) and the accelerate factor, one can obtain the velocity increment u (t1 ) {u (t1 ), a (t1 ), b (t1 )}
and
the
acceleration
increment
(t )} from time t t to time t t t . Then, the (t1 ) {u (t1 ), a(t1 ), b u 1 k k 1
i (t1 ) and the accelerate factor on the displacement increments u i (t1 ) , u i (t1 ) , u
cell node N i will lead its neighboring cell nodes to produce the nodal force k
N increment fi i (t1 ) .
10
N i3 neighbor
k
fiNi (t1) {fu (t1), fa (t1), fb (t1)}
N
4 i
neighbor
crack
u i (t1 ) u i (t1 ), a i (t1 ), b i (t1 )
ui (t1) ui (t1), ai (t1), bi (t1)
N i2 neighbor
cell
Ni
Kih k i k a / b
(t ) i (t1 ) u i (t1 ), ai (t1 ), b u i 1
Mih mi ma / b
f i (t1 ) {f u (t1 ), f a (t1 ), f b (t1 )}
N i1 neighbor Fig. 3 Local updating model for a cell and it neighbors Therefore, the process of the CDCA updating rules is: at time t1 tk t , the increment of nodal force leads to the increment of nodal displacement u i (t1 ) , nodal i (t1 ) , and the increments of nodal velocity u i (t1 ) and nodal acceleration u i (t1 ) lead displacement u i (t1 ) , nodal velocity u i (t1 ) and nodal acceleration u k
N to the increment of nodal force of its neighboring nodes fi i (t1 ) , until the system
dynamic equilibrium is achieved. In other words, the self-organization phenomenon k
N of ui (t1 ) 0 , u i (t1 ) 0 , u i (t1 ) 0 and fi i 0 appears. Therefore, the
updating steps for elastodynamics can be given as: (1) Obtain the initial value of displacement and velocity at time t 0 according to Eqs. (5) and (6); obtain the cell mass matrix and stiffness matrix M i , K i for cell Ni .
(2) The equilibrium equation of the cell node N i can be described as i (t1 ) K i ui (t1 ) fi (t1 ) Mi u 11
(20)
i (t1 ) , u i (t1 ) , and fi (t1 ) are increments of second time derivative of in which u
degrees of nodal freedom, degrees of nodal freedom and nodal force, respectively, of cell node N i . (3) Calculate the increment of degrees of nodal freedom u ih (t1 ) via the increment of nodal force fi (t1 ) and Eq. (20). (4) Via the increment of degrees of nodal freedom u ih (t1 ) and the accelerate factor
, one can obtain the displacement vector u i (t1 ) , which is ui (t1 ) ui (t ) ui (t1 )
(21)
(5) According to u ih (t1 ) and the accelerate factor , the velocity and the acceleration of cell N i can be easily obtained and are
u i (t1 ) u i (t ) u i (t1 )
(22)
i (t1 ) u i (t ) u i (t1 ) u
(23)
(6) Restrict all degrees of nodal freedom at all neighboring cells N ik , which can be seen in Fig. 3. k
N (7) Obtain the nodal force increment fi i (t1 ) of the neighboring cell N ik via
ui (t1 ) and the accelerate factor according to the following equation
i (t1 ) fiNi (t1 ) K iNi ui (t1 ) MiNi u k
k
k
k
(24)
k
N N where K i i and M i i are the stiffness matrix and mass matrix of neighboring cell
N ik , respectively. (8) Finish the calculation of steps (2) – (7) at all cell nodes, until ui (t1 ) 0 , k
N u i (t1 ) 0 , u i (t1 ) 0 and fi i 0 appear.
12
3. Time integration and dynamic stress intensity factors In the present method, an implicit time integration scheme of Newmark method is employed for the time derivative, and the interaction integral theory is developed in this section for solving the dynamic stress intensity factor. 3.1 Implicit time integration scheme The Newmark time integration scheme is widely used in dynamic analysis, because it can be designed to remain unconditionally stable. Taking cell N i as the example, and the time step from time t n to time tn1 tn t , then the time integration scheme can be given as i (tn1 ) K i ui (tn1 ) fi (tn1 ) Mi u
(25)
1 i (t n )t 2 u i (t n1 )t 2 u i (t n1 ) u i (t n ) u i (t n )t ( )u 2
(26)
i (tn )t u i (tn1 )t u i (tn1 ) u i (tn ) (1 )u
(27)
in which and are time difference constants; when 0.5 2 , this method is unconditionally stable, and when 0.5 and 2 , the method is stable if t 1 /(max / 2 ) . According to Eqs. (25) - (27), one can obtain all constants
that the present method uses:
0.5 0.25(0.5 ) 2 2 a0 1 /( t ) a3 1 /( 2 ) 1 a2 1 /( t ) a6 t (1 ) a7 t
(28)
In this case, the equivalent stiffness matrix can be written as ~ K i K i a0M i
(29)
Same as the equivalent stiffness matrix, the equivalent nodal loading vector can be written as ~ i (tn )) fi (tn1 ) fi (tn1 ) M i (a0ui (tn ) a2u i (tn ) a3u
According to Eqs. (29) - (30), one can obtain the equivalent equation 13
(30)
~ ~ K i ui (tn1 ) fi (tn1 )
(31)
Solving Eq. (31), one can obtain the displacement on time t n 1 , and based on the displacement on time step t n 1 , the velocity and the acceleration can be given as i (tn1 ) a0 (ui (tn1 ) ui (tn )) a2u i (tn ) a3u i (tn ) u
(32)
i (tn ) a7u i (tn1 ) u i (tn1 ) u i (tn ) a6u
(33)
Combining dynamic cellular automaton in section 2.3 and Eqs. (29) - (33), one can easily obtain the displacement, velocity and acceleration of cell N i . The same operation can be taken on all other cells, which can be seen that no assembled global stiffness matrix and mass matrix are needed in the whole calculation. Then, the memory is greatly saved, and by employing the fast updating scheme, the calculating efficiency can also be greatly improved. 3.2 Interaction integral for dynamic stress intensity factors In this section, a combination of interaction integral method and equivalent domain integral is employed to evaluate the DSIF[1,3,41]. In this method, the assumption of traction-free boundary condition on the crack surface is considered, and the dynamic equilibrium equation and strain displacement relation are employed. The J-integral [3] can be given as
J [ ij A
u j x1
W1i ]
2u u u 2ui q dA [( 2 j f b ) j i ]qdA A xi t x1 t x1t
(34)
in which q is a weighting function, and it is unity on the internal integral arc 0 and it is equal to zero on the outer arc 1 , which can be seen in Fig. 4, and in the present method, on nodes which are located inside the circle q 1 , on nodes which are located outside the circle q 0 , and on the other location of the integral element 4
q(x) N i (x)qi ; W ij ij , and ij , ij , u j are stress, strain and displacement i 1
on time tn1 tn t . Consider two states for this cracked body, which are State 1: 14
ij(1) , ij(1) , u (j1) are corresponding to the actual states, and State 2: ij( 2 ) , ij( 2) , u (j2) are the auxiliary states derived from Westergaard and Williams[42] for modes I and II cracks. Employing the interaction integral method and equivalent domain integral, the interaction integral for States 1 and 2 can be given as
M
(1, 2 )
[ A
(1) ij
u (j2) x1
( 2) ij
u (j1)
2u (j1) u (j2) q ] dA qdA A x1 xi t 2 x1 (1) ( 2 ) ij ij 1i
(35)
The integral domain can be seen in Fig. 4, which is marked by the shaded element, and h is the average area of elements. is the parameter of the integral domain, and it is usually given as 1 5 . According to the interaction integral method, the left hand term of Eq. (35) can be rewritten as M (1, 2)
2( K I(1) K I( 2) K II(1) K II( 2) ) ~ E
(36)
~ ~ in which E E for the plane stress problem, and E E /(1 2 ) for the plane strain problem. As we know, when K I( 2) 1.0 and K I(I2) 0.0 for State 2, K I can be obtained, and when K I( 2) 0.0 and K II( 2) 1.0 for State 2, K II can be obtained. Eq. (36) can be rewritten as M
(1, 2 )
(1) u (y2) u (y1) u x( 2) q (1) ( 2 ) u x ( 2) [ xy x xy x(1) x( 2) xy(1) xy( 2) ] dA A x x x x x ( 2 ) ( 1 ) u u q u ( 2) u (1) (37) [ xy(1) x y(1) y xy( 2 ) x y( 2) y ] dA A x x x x x 2 (1) ( 2) 2u x(1) u x( 2) u y u y [ 2 ]qdA A t x t 2 x (1) x
15
A
h
crack
R h
crack tip element
0
1 Fig. 4 Integral domain
4. Numerical results and discussions In this section, some numerical examples for elastodynamics problems under dynamic loads with the assumption of plane strain condition of 2D elasticity are considered. Time differential parameters 0.25 and 0.5 are used in those examples. In addition, many results by some other methods are introduced in this paper for comparison. 4.1 Mode-I semi-infinite crack with Heaviside step loading An infinite plate with a semi-infinite mode-I crack is considered in this section, which is loaded by a tensile stress perpendicular to the crack surface. The analytical solution of DSIF for the present example by Freund[43] is employed to validate the accuracy and efficiency of the present method. The geometry of the model is as follows: plate length L 10m , crack length a 5m , and vertical half height of the plate H 2m . The geometry and boundary conditions can be seen in Fig. 5. The material properties are: Young's modulus E 210GPa , Poisson's ratio 0.3 , and density of material 8000kg / m3 . A tensile Heaviside step stress of 0 500MPa is applied on the top surface of the plate, which can be seen in Fig. 6. According to references [3, 20, 44], the total time for the simulation is limited to
t 3tc 1.009 103 s , which is the time for when the tensile stress wave is reflected 16
on the bottom side and reaches again to the crack tip, and tc H / cd , where cd is the dilatational wave speed. For convenience, normalized DSIF K I (t ) /( 0 H ) is used in this example, and the analytical equation can be given as
K I (t ) 2 0 1
0 cd (t t c )(1 2 )
t tc t tc
(38)
0
2H=4m a=5m
L=10m
Fig. 5 Geometry and boundary conditions of a semi-infinite crack In the present method, t 1.2 105 , and the number of time steps is 90. The present method results of mode-I DSIF compared with the results by sES-FEM and analytical solutions are plotted in Fig. 7, in which one can see that the results by the present method agree with the analytical solutions and results by sES-FEM[3], and high accuracy and efficiency can be achieved by the present method.
0
t Fig. 6 Heaviside step loading
17
1.5 CDCA Analytical sES-FEM[3]
1.2
KI(t) / σ0H1/2
0.9 0.6 0.3 0.0 -0.3 0.0
0.3
0.6
0.9
1.2
1.5
1.8
2.1
2.4
2.7
3.0
t / tc
Fig. 7 Normalized mode-I DSIF results by different methods 1.5 101 83 63 41 21 Analytical
KI(t) / σ0H1/2
1.2 0.9 0.6 0.3 0.0 -0.3 0.0
0.3
0.6
0.9
1.2
1.5
1.8
2.1
2.4
2.7
3.0
t / tc
Fig. 8 Normalized mode-I DSIF results by different meshing elements Results from different meshing densities with the present method are shown in Fig. 8, in which the results are much closer to the analytical solutions with the increase in the node numbers. Fig. 9 plots the convergence of the present method with different meshing densities, and it is shown that convergence can be achieved with the increase in the node numbers, and the relative errors become much smaller with much more node numbers.
18
0.5 10141 8333 6325 4117 219
||KI(t)-KIA(t)|| / σ0H1/2
0.4 0.3 0.2 0.1 0.0
-0.1 0.0
0.3
0.6
0.9
1.2
1.5
1.8
2.1
2.4
2.7
3.0
t / tc
Fig. 9 Convergence of DSIF by different meshing elements 400
CDCA XFEM
350
Memory (MB)
300 250 200 150 100 50 0 0
3000
6000
9000
12000
15000
Node numbers
Fig. 10 Memory expenses comparison between the present method and XFEM 1050
CDCA XFEM
900
CPU time (s)
750 600 450 300 150 0 0
3000
6000
9000
12000
15000
Node numbers
Fig. 11 CPU time comparison between the present method and XFEM
Fig. 10 shows the memory expenses by the present method with different node numbers, in which it can be seen that the calculating memories determined by the present method are much smaller than those determined by the XFEM because no assembled global stiffness and mass matrix are needed in the whole calculation, and 19
in XFEM, the recent half bandwidth technique is employed. Fig. 11 plots the CPU time expending by the present method with different node numbers, and one can see that the CPU time by the present method is close to XFEM with much smaller node numbers. However, the CPU time of the present method is much smaller than that of XFEM with much larger node numbers, and the CPU time increases much more gently with the increase in the node numbers. Thus, the advantage is very obvious in the present method for much larger node numbers.
5.2 A slanted edge crack plate with Heaviside step loading To illustrate the efficiency and reliability of the present time domain CDCA for simulating mixed mode dynamic crack problems, a slanted edge crack is considered in this section. A rectangular plate with an edge crack slanted at an angle 45 is employed, which is plotted in Fig. 12. The plate is supported with left, upper and bottom sides on normal directions, and the right side is subjected to a uniform tensile Heaviside step loading 0 , which can be seen in Fig. 6. The geometrical parameters of the plate are: b 44mm , h 32mm , c 6mm , and the crack length
a 22.63mm . The material properties can be given as: shear modulus G 29.4GPa , Poisson ratio 0.286 , and the density of the material 2.45 10 6 kg / mm3 [1,3]. For comparison with the reference [1,3], the total time of the simulation is t 25s . In addition, the normalized DSIF K I,II / 0 a is considered.
20
b
0
h a
c Fig. 12 Geometry and boundary conditions of a slanted edge crack plate 1.8 MFGM[28] CDCA DBEM[46] sES-FEM[3]
1.5
KI
KI,II(t) / σ0(πa)1/2
1.2 0.9 0.6
KII
0.3 0.0 0
5
10
15
20
25
Time t (μs)
Fig. 13 Normalized mixed modes DSIF results by different meshing elements 1.8
α=4.8 α=4.4 α=4.0 α=3.6 α=3.2 α=2.8 α=2.4 α=2.0 α=1.6
1.5
KI(t) / σ0(πa)1/2
1.2 0.9 0.6 0.3 0.0 0
5
10
15
20
25
Time t (μs)
Fig. 14 Normalized mode-I DSIF results by different integral parameters
21
1.2
α=4.8 α=4.4 α=4.0 α=3.6 α=3.2 α=2.8 α=2.4 α=2.0 α=1.6
KII(t) / σ0(πa)1/2
0.9
0.6
0.3
0.0 0
5
10
15
20
25
Time t (μs)
Fig. 15 Normalized mode-II DSIF results by different integral parameters In the present method, t 5.0 107 is used, and time step number is 50, and a total of 5840 nodes and 5785 elements are used in this example. A comparison of the normalized mode-I and mode-II DSIFs by different methods is plotted in Fig. 13, in which one can see that the present method results agree well with the mesh-free Galerkin method (MFGM) by Wen and Aliabadi[28], dual boundary element method (DBEM) by Fedelinski et al. [45] and sES-FEM by Liu et al. [3], and it is verified that a high accuracy of the present method can be achieved. In this example, the influence of integral domain parameter on the equivalent domain integral of DSIF is studied, which can be seen in Figs. 14 - 15, and it can be seen in those figures that a suitable value of the parameter cannot influence the accuracy of the DSIF of the present method, and for 1.6 4.8 , the results are almost the same with each other. Then, the conclusion is found that 1.6 4.8 should be used in order to preserve the accuracy of the present method.
5.3 An oriented central crack plate with Heaviside and blast loading In this example, a rectangular plate with an arbitrarily oriented central crack is considered, and the geometry of the model is shown in Fig. 16. The height of the model is 2H 40mm , the width of the model is 2W 20mm , and the crack length is
2a 4.8mm . Two types of loadings are considered in this section, which are Heaviside step loading and blast loading, and the latter can be seen in Fig. 17. Those loadings are uniform and subjected on the top and bottom sides of the model. The 22
material properties are: Young's modulus E 200GPa , Poisson's ratio 0.3 , and material density 5000kg / m3 . For comparison, the results from the sES-FEM[3] and solutions from the traditional XFEM are shown in the same figures as the present method results, and the total simulating time is t 20s . In addition, different angles of crack are considered, i.e., 15o ,30o ,45o ,60o ,75o . For convenience, the normalized DSIF K I, II / 0 a is considered. In the present method, t 3.0 107 is used, and time step number is 68. A total of 7616 nodes and 7437 elements are used in this example, and the integral parameter 3.2 for the equivalent domain integral is used in the present example.
0
y
2H
2a
x 2W
0 Fig. 16 Geometry and loading of a central crack plate
23
0
t t2
t1
Fig. 17 The blast loading 2.8 2.4
CDCA XFEM sES-FEM[3]
15
2.0
30
KI(t) / 0a)1/2
1.6 45
1.2 0.8 0.4 0.0
60
-0.4
75 0
4
8
12
16
20
Time t (s)
Fig. 18 Comparison of the normalized mode-I DSIF by XFEM, Ses-FEM and the present method for different inclination angles and Heaviside step loading 1.6
CDCA XFEM sES-FEM[3]
45
1.2
KII(t) / 0a1/2
30&60 0.8
0.4 15&75 0.0
-0.4 0
4
8
12
16
20
Time t (s)
Fig. 19 Comparison of the normalized mode-II DSIF by XFEM, sES-FEM and the present method for different inclination angles and Heaviside step loading
24
1.8
15
CDCA sES-FEM[3]
30 1.2 45
KII(t) / 0a)1/2
0.6 0.0
60 75
-0.6 -1.2 -1.8 0
4
8
12
16
20
Time t (s)
Fig. 20 Comparison of the normalized mode-I DSIF by XFEM, Ses-FEM and the present method for different inclination angles and blast loading 1.5
45
CDCA sES-FEM[3]
30
1.2 0.9
KII(t) / 0(a)1/2
0.6 0.3 15
0.0 -0.3 -0.6 -0.9 -1.2 0
4
8
12
16
20
Time t (s)
Fig. 21 Comparison of the normalized mode-II DSIF by XFEM, Ses-FEM and the present method for different inclination angles and blast loading In this example, two types of dynamic loadings are considered. First, the Heaviside step loading is employed, and the normalized mode-I and mode-II DSIFs for the five different inclination angles are presented in Fig. 18 and Fig. 19, respectively. For comparison, the results by sES-FEM[3] and XFEM are employed in those figures, and the results given in those figures generally show good agreement among the different methods utilized. It can be seen that the amplitude of normalized mode-I DSIF increases with the inclination angle, but the maximum amplitude of the normalized mode-II DSIF occurs at
. Fig. 19 reveals that the normalized
mode-II DSIF is almost the same for
and
and for
and
. In addition, a blast dynamic loading is also considered in Fig. 17, and the loading parameters are given as:
,
. The normalized mode-I and mode-II 25
DSIF for blast loading can be seen in Fig. 20 and Fig. 21, respectively, which show that the present method results are close to those by XFEM and sES-FEM[3], and the same phenomenon for different inclination angles can be obtained for the blast loading, in which the maximum amplitude of normalized mode-I DSIF is at and the minimum one is at mode-II DSIF is at
,
, and the maximum amplitude of normalized
.
5. Conclusions In this paper, a time domain continuous-discontinuous cellular automaton method has been proposed for elastodynamic problems and is based on the static continuous-discontinuous cellular automaton method. First, the continuity to discontinuity enriched function technique is developed for dynamic fracture problems, in order to ensure the interpolation property is satisfied. The shifting amendments are employed in the present method. Second, a time domain cellular automaton theory for elastodynamics has been proposed, in which new cell states for dynamic fracture problems, dynamic cell model and dynamic cellular updating scheme are developed. Besides, fast adaptive updating rules based on time domain cellular automaton for the cell and its neighbors are developed. Third, several time-dependent equations for cells are discretized individually by the Newmark method, and rapidly solved by the fast cellular automaton updating method. Furthermore, to analyze the fracture behaviors of elastic solids with dynamic loads, mixed-mode dynamic stress intensity factors are obtained by the interaction integral method and equivalent domain integral. Because of the usages of the discontinuous enriched function technique, time domain cellular automaton theory and fast adaptive updating rules, the cracks are independent of the integral and meshing calculations, and no remeshing is needed for crack dynamic propagation. In addition, no assembled global matrix is needed for the whole calculating process, so the required computing memory can be greatly reduced. Furthermore, the calculation is located at the cell locality, so the local properties can be easily considered. Additionally, the calculating efficiency can be greatly improved 26
via the fast adaptive updating scheme. Finally, some numerical examples for dynamic stress intensity factor problems are studied in this paper, which validate that both the accuracy and efficiency are greatly improved by using the present method. Solutions determined by the analytical method and some other methods are employed to show that a high accuracy can be achieved by the time domain CDCA, and the examples show that the suitable integral domain parameter is between 1.6 and 4.8.
Acknowledgments The work was financially supported by the State Key Research Development Program of China (No. 2017YFC0804203), the International Partnership Program of the Chinese Academy of Sciences (115242KYSB20160017)and the Youth Innovation Promotion Association of CAS(No. 2014304).
References [1] K. Kishimoto, S. Aoki, M. Sakata. Dynamic stress intensity factors using J-integral and finite element method. Engrg. Fract. Mech. 1980, 13: 387-394 [2] J.W. Kelley, C.T. Sun. A singular finite element for computing time dependent stress intensity factors. Engrg. Fract. Mech. 1979, 12: 13-22 [3] P. Liu, T.Q. Bui, C. Zhang, T.T. Yu, G.R. Liu, M.V. Golub. The singular edge-based smoothed finite element method for stationary dynamic crack problems in 2D elastic solids. Comput. Methods Appl. Mech. Engrg. 2012, 233-236: 68-80 [4] J. Rethore, A. Graviouil, A, Combescure. A stable numerical scheme for the finite element simulation of dynamic crack propagation with remeshing. Comput. Meth. Appl. Mech. Engrg. 2004, 193: 4493-4510 [5] F. Chirino, J. Dominguez J. Dynamic analysis of cracks using boundary element method. EngFract Mech. 1989, 34: 1051-1061 [6] D.P.N Kontoni, D.E. Beskos. Boundary element formulation for dynamic analysis of nonlinear systems. Engng Anal. 1988, 5: 114-125 [7] T. Burczynski, T. Adamczyk. Analysis of nonlinear systems in terms of the 27
boundary element method. Mech. Comput. 1988, 7: 149-164 [8]
S.
Ahmad,
P.K.
Banerjee.
Inelastic
transient
dynamic
analysis
of
three-dimensional problems by BEM. Int. J. Num. Meth. Engng. 1990, 29: 371-390 [9] D. P. N. Kontoni, D. E. Beskos. Inelastic dynamic analysis by the boundary element method. In Boundary Elements IX, Proceedings of the 9th International conference on Boundary Element Methods in Engineering(Vol. 2), ed. C.A. Brebbia, W.L. Wendland& G. Kuhn. Springer-Verlag, Berlin, Germany, 1987, 335-351 [10] T. Panzeca, C. Polizzotto, M. Zito. A BEM formulation of the dynamic elastic-plastic structural problem via variational principle, In Boundary Elements in Machanical and Electrical Engineering, Proceedings of the International Boundary Element Symposium, ed. C.A. Brebbia& A. Chaudouet-Miranda, Computational Mechanics Publications, Southanpton, UK and Springer-Verlag, Berlin, Germany, 1990, 193-204 [11] J.A.M Carrer, J.C.F. Telles. Transient dynamic elastoplastic analysis by the boundary element method. In Boundary Element Technology VI, Proceedings of the 6th Int. Conf. on Boundry Element Technology, ed. C.A. Brebbia. CMP, Southanpton, UK, and Elsevier, London, UK, 1991, 265-277 [12] J.P. Agnantiaris, D. Polyzos, D.E. Beskos. Some studies on dual reciprocity BEM for elastodynamic analysis. Comput. Mech. 1996, 17: 270-277 [13] J.P. Agnantiaris, D. Polyzos, D.E. Beskos. Free vibration analysis of non-axisymmetric and axisymmetric structures by the dual reciprocity BEM. Eng. Anal. Bound. Elem. 2001, 25: 713-723 [14] M. F. Samaan, Rashed Y. F. BEM for transient 2D elastodynamicsusing multiquadric functions. Int. J. Solids Struct. 2007, 44: 8517-8531 [15] J. Dominguez, R. Gallego. Time domain boundary element method for dynamic stress intensity factor computations. Int. J. Num. Meth. Eng. 1992, 33: 635-647 [16] H. Gao, P. Klein. Numerical simulation of crack growth in an isotropic solid with randomized internal cohesive bonds. J. Mech. Physics Solids. 1998, 42: 187-218 [17] T. Belytschko, H. Chen, J.X. Xu, G. Zi. Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment. Int. J. Numer. Meth. Eng. 28
2003, 58: 1873-1905 [18] J. Chessa, T. Belytschko. Arbitrary discontinuities in space-time finite element by level sets and X-FEM. Int. J. Numer. Meth. Engrg. 2004, 61: 2595-2614 [19] T. Menouillard, J. Rethore, A. Combescure, H. Bung. Efficient explicit time stepping for the extended finite element method(XFEM). Int. J. Numer. Meth. Engrg. 2006, 68: 911-939 [20] T. Elguedj, A. Gravouil, H. Maigre. An explicit dynamics extended finite element method. Part 1: Mass lumping for arbitrary enrichment functions. Comput. Meth. Appl. Mech. Eng. 2009, 198: 2297-2317 [21] J. Rethore, A. Gravouil, A. Combescure. An energy-conserving scheme for dynamic crack growth using the extended finite element method. Int. J. Numer. Meth. Engng. 2005, 63: 631-659 [22] L. Wu, P. Liu, C. Shi, Z. Zhang, T. Q. Bui, D. Jiao. Edge-based smoothed extended finite element method for dynamic fracture analysis. Appl. Math. Model. 2016, 40: 8564-8579 [23] C. A. Durate, O. N. Hamzeh, T. J. Liszka, W. W. Tworzydlo. A generalized finite element method for the simulation of three-dimensional dynamic crack propagation. Comput. Meth. Appl. Mech. Eng. 2001, 190: 2227-2262 [24] W. Li, H. Zheng, G. H. Sun. The moving least squares based numerical manifold method for vibration and impact analysis of cracked bodies. Eng. Fract. Mech. 2018, 190: 410-434 [25] T. Belytschko, M. Tabbara. Dynamic fracture using element-free Galerkin methods. Int. J. Numer. Meth. Eng. 1996, 39: 923-938 [26] P. Krysl, T. Belytschko. Then Element free Galerkin method for dynamic propagation of arbitrary 3-D cracks. Int. J. Numer. Meth. Engng. 1999, 44: 767-800 [27] Y. Lu, T. Belytschko, M. Tabbara. Element-free Galerkin method for wave propagation and dynamic fracture. Comput. Meth. Appl. Mech. Eng. 1995, 126: 131-153 [28] P. H. Wen, M. H. Aliabadi. Evaluation of mixed-mode stress intensity factors by mesh-free Galerkin method: static and dynamic. J. Strain Anal. Eng. 2009, 44: 29
273-286 [29] T. Rabczuk, J.H. Song, T. Belytschko. Simulations of instability in dynamic fracture by the cracking particles method. Eng. Fract. Mech. 2009, 76: 730-741 [30] F. Yan, Y.H. Wang, L.G. Tham, Y.K. Cheung. Dual reciprocity hybrid boundary node method for 2-D elasticity with body force. Eng. Anal. Bound. Elem. 2008, 32: 713-725 [31] F. Yan, Y. Miao, Q.N. Yang. Quasilinear hybrid boundary node method for solving nonlinear problems. CMES-Comput. Model. Eng. Sci. 2009, 46: 21-50 [32] F. Yan, X. T. Feng, H. Zhou. A dual reciprocity hybrid radial boundary node method based on radial point interpolation method. Comput. Mech. 2010, 45: 541-552 [33] F. Yan, X. T. Feng, J. H.Lv, P. Z. Pan. A new hybrid boundary node method based on Taylor expansion and Shepard interpolation method. Int. J.Numer. Meth. Eng. 2015, 102: 1488-1506 [34] F. Yan, Y.H. Wang, Y. Miao, L.G. Tham, Y.K. Cheung. Dual reciprocity hybrid boundary node method for free vibration analysis. J. Sound Vib. 2009, 321: 1036-1057 [35] Y. Miao, T. G. He, H Luo, H. P. Zhu. Dual hybrid boundary node method for solving transient dynamic fracture problems. CMES-Comput. Model. Eng. Sci. 2012, 85: 484-498 [36] X.T. Feng, P.Z. Pan, H. Zhou. Simulation of the crack microfracturing process under uniaxial compression using an elasto-plastic cellular automaton. Int. J. Rock Mech. Min. Sci. 2006, 43: 1091-1108. [37] F. Yan, X.T. Feng, P.Z. Pan, S.J. Li. Discontinuous cellular automaton method for crack growth analysis without remeshing. Appl. Math. Model. 2014, 38: 291-307 [38] F. Yan, X. T. Feng, P. Z. Pan, S. J. Li. A continuous-discontinuous cellular automaton method for cracks growth and coalescence in brittle material. Acta Mech.Sinica. 2014, 30: 73-83 [39] F. Yan, X. T. Feng, J. H.Lv, P. Z. Pan, S. J. Li. Continuous-discontinuous cellular automaton method for cohesive crack growth in rock. Eng.Fract. Mech. 2018, 188: 361-380 30
[40] F. Yan, X.T. Feng, P.Z. Pan, Y.P. Li. A continuous-discontinuous cellular automaton method for regular frictional contact problems. Arch. Appl. Mech. 2013, 83: 1239-1255. [41] S.H. Song, G. H. Paulino. Dynamic tress intensity factors for homogeneous and smoothly heterogeneous materials using the interaction integral method. Int. J. Solids Struct. 2006, 43: 4830-4866 [42] N. Moes, J. Dolbow, T. Belytschko. A finite element for crack growth without remeshing. Int. J. Numer. Mehods Engrg. 1999, 46: 131-150 [43] L.B. Freund. Dynamic fracture mechanics. Cambridge University Press, London, 1990 [44] T. Menoulillard, T. Belytschko. Dynamic fracture with meshfree enriched XFEM. Acta Mech. 2010, 213: 53-69 [45] P. Fedelinshi, M.H. Aliabadi, D.P. Rooke. The laplace transforms DBEM for mixed-mode dynamic crack analysis. Comput. struct. 1996, 59: 1021-1031
31
Highlights: A continuity-to-discontinuity enriched function technique is developed for elastodynamics. Fast adaptive cellular automaton updating rules for the cell and its neighbors are developed. A time domain cellular automaton theory for elastodynamics is proposed. The cracks are independent of the integral meshes and no global matrix is needed.
32