- Email: [email protected]
Lattice spring model with angle spring and its application in fracture simulation of elastic brittle materials
Journal Pre-proofs Lattice spring model with angle spring and its application in fracture simulation of elastic brittle materials Tao Wang, Mao Zhou, Yongqiang Li, Yin Yu, Hongliang He PII: DOI: Reference:
S0167-8442(19)30479-3 https://doi.org/10.1016/j.tafmec.2019.102469 TAFMEC 102469
To appear in:
Theoretical and Applied Fracture Mechanics
Received Date: Revised Date: Accepted Date:
28 August 2019 28 December 2019 28 December 2019
Please cite this article as: T. Wang, M. Zhou, Y. Li, Y. Yu, H. He, Lattice spring model with angle spring and its application in fracture simulation of elastic brittle materials, Theoretical and Applied Fracture Mechanics (2020), doi: https://doi.org/10.1016/j.tafmec.2019.102469
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. 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.
© 2020 Elsevier Ltd. All rights reserved.
Lattice spring model with angle spring and its application in fracture simulation of elastic brittle materials Tao Wanga, Mao Zhoua, Yongqiang Lia,b, Yin Yuc, Hongliang Hec a. College of Science, Northeastern University, Shenyang 110819, China b. Key Laboratory of Ministry of Education on Safe Mining of Deep Metal Mines, Northeastern University, Shenyang, Liaoning 110819, China. c. National Key Laboratory of Shock Wave and Detonation Physics, Institute of Fluid Physics, CAEP, 621999 Mianyang, China
ABSTRACT Lattice Spring Model with Angle Spring (LSMA) is a novel type of lattice model. Compared with the general lattice spring model, LSMA adds the rotational and tangential degrees of freedom, whose spring parameters are obtained by finite element parameter mapping. The traditional lattice spring model only consider the normal interaction between two points (the Poisson’s ratio is fixed), which limits the application of the lattice spring model in a wider range of materials. Some scholars proposed to add shear spring into the traditional model to solve the problem of fixed Poisson’s ratio, but the rotation invariance could not be maintained. The main reason is that the shear spring can not distinguish the difference of tangential velocity (or displacement) between the two particles due to the common rotation or shear. Therefore, the rotation of the whole rigid body may cause incorrect generation of additional strain energy inside the shear spring. LSMA model is able to maintain rotation invariance and reproduce different Poisson’s ratio because of the spring parameters with rotational degrees of freedom and the rotation effect of lattice points. In this paper, the finite element stiffness matrix with rotational freedom is derived by using LSMA as the basic mechanical model, and the parameter mapping theory of spring stiffness coefficient is introduced. By means of numerical simulation, the deflection and rotation Angle of cantilever beam, shear wave and compressional wave velocity of stress wave and the simulation of compact tensile test are checked and calculated respectively, and the correctness of the model is verified. Keywords: Angle spring, Lattice spring model, Elastic brittle material, Poisson's ratio
1. Introduction Rock is a typical elastic brittle material, and the dynamic fracture is the core problem of earthquake, landslide, rock blasting and other mechanical behaviors. Therefore, it is of great theoretical and practical significance to simulate the dynamic fracture effectively. There are two reasons, first, on the level of mechanics, the dynamic fracture mechanism of rock is not very clear, there is not a perfect mechanical model can accurately describe the dynamic fracture mechanism of rock. Second, at the algorithm level, how to simulate dynamic crack generation, propagation and instability is still a very difficult problem. There are two main types of numerical simulation methods for rock dynamic fracture: continuous medium method and discrete medium method. Although the theory of dynamic fracture mechanics of continuous media has made good progress, the theoretical system still has many difficulties in simulating dynamic fracture of rock. Traditional continuum mechanics is a theoretical system based on the continuous assumption of field, regardless of the microstructure of the material. When simulating the fracture of material, the external fracture criterion should be used. The selection of dynamic fracture criteria is a very difficult problem. In addition, when there is crack propagation, The processing of finite element mesh is another difficult problem of the theoretical system. Compared with the continuous medium mechanics method, the discrete medium method has a unique advantage
Corresponding author at: College of Sciences, Northeastern University, Shenyang 110819, China. Tel.: +86- 24-
83672653. E-mail address: [email protected] (Yongqiang Li). 1 / 16
in simulating the dynamic fracture of rock. By using a discrete system instead of a continuous medium, the 3-D material failure problem is transformed into a one-dimensional problem. The problems of continuum fracture criterion and the transition from continuous medium to discontinuous medium are avoided, which greatly simplifies the fracture simulation. At present, the mainstream discrete medium models mainly include lattice model [1,2], spring model [3] and particle flow model [4]. Lattice model as a discrete medium method has a strong advantage in simulating dynamic fracture. Lattice model was first proposed by Hrennikoff [5] to solve classical elastic mechanics problems. In this method, the continuum material is equivalent to a discrete structural system, the generation and evolution mechanism of micro-damage in the material is regarded as the failure of a single bar element or spring element in the discrete structural system, the model elements adopt simple constitutive relation and strength criterion. In this way, the problem of determining the constitutive relation of materials and the selection of dynamic fracture criterion is avoided, which provides a simple and effective method for the simulation of dynamic fracture of materials. Early in the lattice model development [6, 7], connections in discrete structures are made of springs, due to the spring model only consider normal interactions between two points, so it limits the material Poisson's ratio (1/4 for plane strain, 1/3 for plane stress), which restricts the application of the model. In order to break through this limitation, some scholars have added additional shear springs into the traditional model [8-10]. The introduction of shear springs can simulate any Poisson’s ratio. However, introducing shear springs into the classical LSM may not maintain rotational invariance [11-13]. Jagota and Scherer[14] believe that the main reason is that the shear spring cannot distinguish the difference between the tangential velocity (or displacement) of the two particles due to their joint rotation or shear, as shown in Fig.1. As a result, global rigid body rotation may incorrectly cause additional strain energy to be generated within the shear spring.
A
ks
B
A
rotation
ks
B
shear
Fig.1 Schematic showing incapability of shear spring to distinguish differences in tangential displacements of two particles owing to rigid rotation or shear. The literature [16,17] used the beam element to replace the spring in the particle model, thus the particle turning effect is taken into account, so that the lattice model could simulate an arbitrary Poisson’s ratio. However, the calculation of rotation Angle is indirectly calculated by lattice point displacement, which requires a large amount of calculation. Zhao et al. [18] considered the normal deformation and shear deformation effects of the spring in the spring model, and solved the fixed Poisson's ratio problem of the lattice model. In this method, the shear deformation directly calculated by the particle displacement still included the rigid body rotation. The literature[18] used local strain method to eliminate the influence of rigid body rotation. In order to consider the influence of bond Angle, Wang [19,20] proposed the concept of angle spring, which can be used to reproduce variable Poisson’s ratio, so that the lattice model has a great advantage in simulating dynamic fracture. In the traditional grid method, the parameter calibration of the bond is related to the specific case. There is no uniform micro-macro parameter correspondence, and the lattices of different topologies have different micro-macro parameter relationships. These all make the parameter calibration problem of the three-dimensional lattice model very complicated. In addition, the traditional lattice model cannot reproduce the microscopic structure of rock. In recent years, the lattice model has been constantly improved in theory, and it has been widely used in the process of dynamic fracture simulation, such as [21-28]. 2 / 16
uA
vB
kn
A
B
uB
kr
ks
A
B
A
B
A
vA
normal
B
shear
angle
Fig.2 Lattice spring model with angular spring In this paper, the lattice spring model with Angle spring is established, as shown in Fig.2. In the figure, u and v are translational degrees of freedom along the normal and tangential directions, and φ is a rotational degree of freedom, with these three degrees of freedom coupling on the lattice. The quantitative mapping between the stiffness matrix of the finite element and the stiffness coefficient of the lattice spring model is realized by using the finite element of three joints and three degrees of freedom with rotational degrees of freedom as the medium and the model mesh shared by the lattice spring model and the finite element model. In this paper, the finite element stiffness matrix with rotational degrees of freedom is derived, the parameter mapping theory of spring stiffness coefficient is explained, and the fracture criterion of spring is introduced. To verify the correctness of the model, the deflection and rotation Angle of the cantilever beam, the shear and longitudinal wave velocities of the stress wave and the simulation of the compact tensile test are checked by means of numerical simulation.
2 3-node finite element stiffness matrix with rotational degrees of freedom In this section, the rotation angles of six-node triangular elements are introduced to generate interpolation functions of triangular elements with rotational degrees of freedom. The six-node linear strain triangular element is shown in Fig.3. Nodes 1, 2 and 3 are intermediate nodes of each side. The displacement of each node along the x and y axes is represented by u and v respectively. αjm, αmi and αij represent the angle between the outer normal of each edge and the x-axis. m um , vm n jm nmi
mi
m(U m ,n , U m , s )
jm
n jm
1 u1 , v1
jm
2 u2 , v2
1(U1,n ,U1,s ) j u j , v j
y
x
i ui , vi
3 u3 , v3
ij nij
Fig. 3 Six-node linear strain triangular element.
x1
y1
j(U j ,n ,U j ,s )
Fig. 4 The local coordinate system
For the six-node line strain triangle element shown in Fig.3, rotational degrees of freedom are introduced on this element and eliminate the middle node 1, 2 and 3 of the edge. According to the interpolation method adopted by Allman, the tangential displacement Us along each side of the triangle varies linearly, while the normal displacement Un along the side changes quadratically. Therefore, in order to establish the local coordinate system along tangential direction and normal direction, coordinate transformation is required before interpolation. The specific transformation method is as follows: the local coordinate yl axis direction obtained after transformation is consistent with the external normal direction of each side, and the local coordinate x1 axis direction is consistent with the tangent direction of each side. The coordinate transformation is shown in Fig. 3 and Fig. 4, the displacement relation of the normal displacement and tangential displacement transformation of each side in the triangular element under the rectangular coordinate system is 3 / 16
u cos v sin
sin U n , jm , ij , mi cos U s
(1)
After the coordinate transformation by Eq.(1), the corner nodes are used to represent the middle nodes of the edge and the middle nodes are eliminated, interpolating in local coordinates respectively. Since the tangential displacement varies linearly along the edge, the tangential displacement of the middle node of each edge is
U 1, s
U j ,s U m,s 2
, U 2, s
U i ,s U m,s 2
U i,s U j ,s
, U 3, s
(2)
2
The normal displacement is maintained as a quadratic interpolation along the edge, Let the interpolation function be a quadratic parabola.
U n = a1 x12 + a2 x1 + a3
(3)
The rotation (slope) on that side is φ, so
=
dU n = 2a1 x1 + a2 dx1
(4) um
m
m
vm
uj j
ui
i
i
j
vj
vi
Fig. 5. Three-node triangular element with nine degrees of freedom Let nodes i, j, m have rotation parameters φi, φj and φm, as shown in Fig.5. The second derivative of Un can also be represented by the rotation parameters of nodes, taking edge j-m in Fig.4 as an example
j - m d 2U n d = = 2 a = 1 dx1 l jm dx12
(5)
where ljm is the length of the side j-m, there are
a1 =
j - m
(6)
2l jm
Substituting the local coordinates (0,Uj,n) and (ljm,Umn) of node j and node m into Eq.(3) respectively, there are
U j,n = a3 , U m,n = a1l 2jm + a2l jm + a3
(7)
According to Eq.(6) and Eq.(7), the coefficients a1, a2 and a3 of Eq.(3) can be solved
a1 =
j - m 2l jm
, a2 =
U m,n - U j ,n l jm
-
j - m 2
, a3 = U j ,n 4 / 16
(8)
Substituting Eq. (8) into Eq. (3) and taking x1=ljm/2, we can obtain the normal displacement of the middle node 1 of edge j-m as U j ,n + U m,n
U1,n =
2
+
m - j 8
(9)
l jm
Similarly, the normal displacements of the intermediate nodes 2 and 3 of edges m-i and i-j can be obtained as follows
U 2, n
U m,n U i ,n 2
i m 6
lmi , U 3,n
U j ,n U i ,n 2
j i 6
lij
(10)
where lmi and lij are the side lengths of edges m-i and i-j respectively. The interpolation function of a six-node triangular element is 3 u N u N u N u N k uk i i j j m m k 1 3 v N v N v N v N v i i j j m m k k k 1
(11)
and
l jm sin jm x j xm ci ,i j m l jm cos jm ym y j bi
(12)
Based on Eq. (1), Eq. (2), Eq. (9) , Eq. (10) , Eq. (11) and Eq. (12), the interpolation function of triangular elements with rotational degrees of freedom can be expressed as
2 u Li ui L j u j Lmum 3 i j bm Li L j j m bi L j Lm m i b j Lm Li v L v L v L v 2 c L L c L L c L L i i j j m m i j m i j j m i j m m i j m i 3
(13)
where, Li, Lj and Lm are area coordinates. The three nodes of the trilateral midpoint are eliminated. The three-node triangular element with nine degrees of freedom, as shown in Fig.5, is obtained by introducing rotational degrees of freedom into the vertex of the triangular element. From Eq.(13) can be deduced
u
i
vi uj vj um vm u1 v1 u2 v2 u3 v3
T
(14)
Tui vi i uj vj j um vm m
T
where T is the conversion matrix of 12×9.
5 / 16
0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 T 12 0 0 0 bm 6 1 2 0 bm 6 0 0 1 2 am 6 0 0 0 0 1 2 am 6 0 0 0 12 0 bi 6 1 2 0 bi 6 0 0 0 0 1 2 ai 6 0 1 2 ai 6 0 0 12 0 bj 6 1 2 0 b j 6 0 0 1 2 a j 6 0 0 0 0 1 2 a j 6
(15)
e
Based on Eq.(14), relational expression of K 99 (stiffness matrix of triangular element with three nodes and e
nine degrees of freedom) and K 1212 (stiffness matrix of triangular element with three nodes and twelve degrees of freedom) can be expressed as e K 9e9 TT K12 12 T
where K 1212 t e
(16)
B DBdxdy , T
A
B Bi
Bj
Bm
B1 B 2
B3
3 Spring parameter mapping theory Whether the lattice point spring model can accurately evolve the response of material under external load depends on the selection of the mass of lattice point and the elastic coefficient of spring, so that the spring network is just consistent with the elastic property of the actual sample. According to the cell size and material density, the lattice mass can be determined by the method of equivalent mass. The core problem of lattice spring model modeling has been transformed into the mathematical mapping relationship between the elastic constitutive relation of actual material and the elastic coefficient of spring network. Gusev [29] reported a method of mapping finite element parameters to discrete element models. When material satisfies elastic deformation, the relationship between node force and node displacement in finite element models can be established by finite element theory.
F Kδ
is global stiffness matrix,
Where F Fi is node force vector, K k ij T
(17)
δ i is node displacement T
vector, i, j 1, 2,..., n , n is degree of freedom. According to Eq. (17), the force of a node i can be rewritten as n
Fi kij j ki11 ki 2 2 kii i kin n
(18)
j 1
The element kii on the principal diagonal of the global stiffness matrix is equal to the opposite number of all elements` sum except the element of principal diagonal in the same row[30], i.e., n
k ii k im
(19)
m 1 m i
6 / 16
Based on Eq. (19), Eq. (18) can be rewritten as n
Fi kim m i
(20)
m 1 mi
As can be seen from Eq.(20), The force of a node is equal to the sum of the product of the adjacent spring parameters and the displacement difference of the node. According to this principle, the stiffness coefficients in the finite element are mapped to the elastic coefficients of the springs in the lattice spring model, and the values of the spring network in the lattice spring model are determined. This parameter mapping method has a rigorous mathematical logic, to ensure that the model can quantitatively represent the mechanical properties of real materials in the stage of elastic deformation, which can solve the difficulty of lacking strict demonstration in the parameter selection of discrete element method model. At the same time, the model does not depend on the background grid when calculating, which avoids the difficulty of numerical calculation caused by the mesh distortion when solving the extremely large deformation and fracture problems by the finite element method. Lattice spring model can impose arbitrary boundary, which effectively solves the problem that mesh-less method can’t impose displacement boundary. The given size of the basic element model is small enough and the number of elements is large enough, the lattice spring model can characterize the dynamic response of the sample from both mesoscopic and macroscopic levels, providing the necessary physical knowledge for the damage evolution of materials.
4. Criterion for spring fracture of elastic brittle materials In order to use the LSM model to simulate the fracture process of the material, a failure criterion needs to be defined to identify the key elements to be removed from the model. Most failure criteria used so far in LSM have been classified into two categories: critical stress (force) criterion and critical strain (displacement) criterion. The critical strain criterion is adopted for calculation in this paper. Because rock and concrete are typical elastic-brittle materials, the constitutive relationship between force F and strain of spring in fracture simulation can be shown in Fig. 6.
F
o
t
Fig.6 Constitutive relationship between force F and strain of spring
t is the critical tensile strain of the spring. Before the spring breaks, the spring is linear elastic. When the principal strain of the spring exceeds t value, the spring breaks. Pure compression (such as hydrostatic compression of uniform media) does not cause material damage, so the spring compression failure is not considered in the calculation. When the fracture criterion is satisfied, the spring will be broken permanently, and the tension and shear resistance between the two lattices will be lost, which is equivalent to forming a small crack in the sample.
5. Kinematics integral In the process of model calculation, the translational and rotational accelerations of each lattice in the unit time step can be obtained, and then the displacement, velocity, angular velocity and angular displacement of each 7 / 16
particle in the next time step can be obtained by acceleration-time integral. The basic principle of the motion equation integration algorithm in lattice spring model and the molecular dynamics is to perform Taylor expansion on the motion equation, i.e.[31],
r (t t ) r (t ) t
r t 2 2 r t 3 3 r ... t 2! 2t 3! 3t
(21)
where r , t , t are the coordinates, time and time step respectively. Generally, when the order of the integration algorithm is expanded very high, the longer the time step of each iteration can be obtained, but at the same time, the more motion information need to be stored in integral operation. Verlet algorithm and its substitution are commonly used low-order algorithms. Although it is necessary to shorten the time step to ensure sufficient accuracy, the particle trajectories calculated by Verlet algorithm are usually time reversible, and the total energy of the system rarely drifts in long-term operation. This paper chooses an alternative method of Verlet algorithm-Beeman algorithm,
4 f (t ) f (t t ) 2 ( ) ( ) ( ) r t t r t v t t t 6m v(t t ) v(t ) 2 f (t t ) 5 f (t ) f (t t ) t 6m
(22)
Where f(t) is the resultant force of the particle at time t, m is the mass of the particle, v is the velocity of the particle. The particle trajectory calculated by Beeman algorithm is the same as Verlet algorithm, but it is more accurate in calculating particle velocity and the total energy conservation property of the model is better.
6. Verification 6.1 Deflection and angle deformation calculation of rectangular section cantilever beam In order to verify the correctness of the model, a rectangular cantilever beam structure is selected. The geometry and boundary conditions of the problem are described in Fig.7. The left end of the beam is fixed in all directions. A force of 10 N is applied at the right end of the beam. The elastic constants of the material are: Young’s modulus E=200GPa and Poisson’s ratio μ= 0.22. Based on the elasticity theory, the maximum deflection and rotation angles are obtained respectively，
u y max
Fl 3 , 3EI
max
Fl 2 2 EI
(23)
Where I bh3 /12 , the calculation results of the maximum deflection and the maximum rotation angle in the Y direction with different number of lattices are shown in Table 1. Suppose the number of x-axis lattice points is mx; y is my y
x
h b
l
F
Fig.7 cantilever beam structure diagram (l=1.0m,b=h=0.01m)
8 / 16
Table 1 Calculation results of maximum deflection and rotation angle of cantilever beam under different number of nuclear power plants max (rad)
uymax (m) Lattice (mx×my)
LSMA Calculate values
without angle spring Error(%)
Calculate values
Error(%)
LSMA Calculate values
without angle spring Error(%)
Calculate values
20×2
0.015
25.0
0.0018
91.0
0.025
16.7
-
100×2
0.017
15.0
0.0141
29.5
0.026
13.3
-
200×2
0.018
10.0
0.0179
10.5
0.028
6.7
-
200×4
0.019
5.0
0.0181
9.5
0.029
3.3
-
400×4
0.020
0.0
0.0194
3.0
0.030
0.0
-
Theoretical values
0.020
0.030
In order to test the accuracy of LSMA, the comparison results among the calculated values of LSMA, the calculated values of the model without angle spring and the theoretical calculated values are given in Table 1. The relative error in Table 1 refers to the error between the numerical simulation solution and the theoretical solution. As can be seen from the comparison results in Table 1, the calculation accuracy of the LSMA is much higher than the model without angle spring. Because the spring stiffness coefficient in the lattice spring method without angle spring is calculated by the 3-node 6-degree-of-freedom triangular finite element. The 3-node triangular finite element is a constant strain element, which has good reliability and low calculation accuracy as a coordination element. However, the spring stiffness coefficient in the LSMA is calculated by the 3-node, 9-degree-of-freedom uncoordinated triangular finite element, which is obtained from the original Allman element, keeping the accuracy of the original Allman element. It can also be seen from Table 1 that the calculation accuracy of lattice spring method increases greatly with the increase of lattice points. 6.2 Large deflection bending deformation calculation of slender cantilever beam The lattice points in the traditional lattice-spring model are only subject to normal forces (the force is only related to the relative distance between two lattice points), which is called the "Hookean model"[32]. In Hookean model, there is no resistance to shear deformation between two lattice points. As a result, Hookean model only has fixed Poisson's ratio, which cannot simulate various Poisson's ratio materials. To break this limitation, several improved models have been proposed, the most popular model of which is the "Born model"[32]. Born model adds shear forces between two lattice points, which can control the overall Poisson's ratio by adjusting the tangential spring stiffness coefficient. One disadvantage of the Born model is that it does not have the rotational conservation of strain energy[33]. This shortcoming is not prominent in the study of small deformation, but it will have an impact on large deformation phenomena. The LSMA proposed in this paper by introducing the angle spring, the Born model can obtain the rotational conservation of strain energy automatically[33]. To verify LSMA's ability to calculate large deformation, the slender cantilever beam with rectangular section is selected as shown in Fig.7. The cantilever beam, of length l=10m, is characterized by a rectangular cross-section, of thickness h=0.1m and of width b=1m. It is assumed to be constituted of a linear elastic material, exhibiting Young’s modulus E =1.2×106 and Poisson’s ratio μ=0. The loading history consists in the application of a uniform follower load q = q0EI/l3 , where I is the cross-section moment of inertia. q0 is the load coefficient, increasing from zero up to the value q0max = 40.
9 / 16
1.0
q0 15.0
0.5
q0 10.0
q0 25 .0
q0 30.0
q0 5.0
q0 35 .0
Tip deflections(/l)
q0 20.0
0.0
utip,x
-0.5
utip,y
-1.0
utip,y theoretical
utip,x theoretical
q0 40.0 y
-1.5
x
0
10
20
30
40
Applied load coefficient q0
(a)
(b)
Fig.8. Clamped cantilever subjected to uniformly distributed follower load: (a) reference and deformed configurations at eight uniformly spaced load levels and (b) load–displacement curve.
Fig.8(a) shows the large-displacement deformed configuration at eight uniformly spaced load levels, where
usum u x2 u y2 , ux and uy is the displacement in the x and y directions. Fig.8(b) shows the load-displacement curve corresponding to x-deflection utip,x and y-deflection utip,y are reported. The comparison of the numerical results for the considered quite slender structure with the reference solution reported in [34] proves the effectiveness of the present formulation. 6.3 Wave velocity test In solid materials, both shear and longitudinal waves can be effectively propagated. The propagation velocities e
of shear and longitudinal waves depend only on the properties of materials, which can be represented by CL and
CRe in the one-dimensional elastic strain wave, i.e.,
CLe
E 1 E e , CR 2 1 1 1 2 1
(24)
where, is the density of the material. The theoretical values of CL and CR can be calculated by Eq. (24). The calculated values of shear and longitudinal waves velocity calculated by the LSMA and without angle spring model in section 6.1 example are shown in Table 2. e
e
Table 2 Shear and longitudinal wave velocities of different materials Shear waves velocity (m/s)
Longitudinal waves velocity(m/s)
without angle without angle LSMA LSMA ρ spring spring (kg/m3) Theoretical Theoretical values values Calculated Error Calculated Error Calculated Error Calculated Error values (%) values (%) values values (%) (%)
E (GPa)
μ
70
0.33
2700
3121.95
3180.09
1.86
3201.25
2.54
6197.82
6167.34
0.49
6153.82
0.71
120
0.32
7800
2414.02
2390.63
0.97
2380.95
1.37
4692.02
4625.83
1.41
4601.46
1.93
150
0.22
8010
2770.35
2690.89
2.87
2655.10
4.16
4623.84
4626.3
0.05
4632.16
0.18
210
0.25
7800
3281.65
3207.89
2.25
3181.89
3.04
5683.99
5626.38
1.01
5572.58
1.96
10 / 16
The calculation results show that the maximum error between the LSMA model and the without angle spring model and the theoretical values is 2.87% and 4.16% respectively for the shear and longitudinal wave velocities of materials, indicating the accuracy of the LSMA model. 6.4 Compact Tensile Test Simulation Dynamic crack propagation is closely related to loading rate. In order to verify whether LSMA model can reproduce the effect of loading rate on dynamic crack growth, this model is used to simulate compact tensile test. Geometry of the compact tension specimen is length-width = 0.2m×0.2m. The length and width of the notches are 0.064 m and 0.018 m respectively. The left end of the notch is fixed by the steel sleeve, and the right end is loaded by the displacement control of the same steel sleeve. Similar experiments have been done with concrete in literature [35]. The results show that different failure modes appear with different loading rates. The cracking modes of specimens under different loading rates are simulated by this example. The whole specimen is simulated with 6468 springs. The geometric, LSMA model mesh and boundary conditions of the specimen are shown in Fig.9. In order to more accurately reproduce the contact between the specimen and the loader during the actual loading process, a steel plate with relatively high stiffness is installed in the loading area. The displacement control load is applied on the steel plate, and the contact between the steel plate and the specimen is controlled by the bonding interface unit.
v
0.2m
0.2m
Fig. 9 Geometry of the compact tension specimen, LSMA model mesh, loading (displacement control) and boundary conditions. According to literature[35], the numerical simulation is carried out for normal strength concrete with the following macroscopic properties: Young’s modulus E = 36 GPa, Poisson’s ratio μ = 0.18, the critical tensile strain t= 0.116×10-3. The specimen is analyzed assuming dynamic loading condition with the following displacement rates v=0.304m/s, 1.375m/s and 3.318m/s. The experimental results of the three loading rates are shown in Fig. 10(a), the simulation results are shown in Fig.10(b) by LSMA, and the simulation results are shown in Fig.10(c) without angle spring model. As can be seen from Fig.10, the simulation results by LSMA are basically consistent with the experimental results. When the loading rate is small (i.e. 0.304m/s), the crack propagates along the central line, and the failure mode is similar to that of quasi-static loading; when the loading rate increases to 1.375m/s, the crack propagates towards the loading end gradually; when the loading rate continues to increase to 3.318m/s, the crack propagates appeared bifurcation phenomenon. With the increase of loading rate, the fracture mode of specimens gradually transits from mode I to mixed-mode.
11 / 16
v=0.304 m/s
v=1.375 m/s
v=3.318 m/s
(a) experimental result[32]
(b) LSMA
(c) without angle spring
Fig.10 Comparison of experimental results and simulation results of crack path under different loading velocities
7 Conclusions In this paper, the LSMA model uses the spring parameters with rotational degrees of freedom and considers the angular effect of the lattice points. It can keep the rotation invariance and reproduce different Poisson's ratios. Based on the finite element model with three nodes and nine degrees of freedom including rotational degrees of freedom, by means of the mesh model shared by lattice spring model and finite element model, the stiffness matrix of finite element is mapped to the spring stiffness coefficient of lattice spring model, so that the selection of spring stiffness coefficient can be deduced mathematically. In order to verify the correctness of the model, the deflection and rotation of the cantilever beam, the shear and longitudinal wave velocities of stress wave and the compact tensile test simulation are checked by numerical simulation. The lattice spring model with angle spring is more accurate than the model without angle spring.
Acknowledgement This work was supported by the National Natural Science Foundation of China (No. 11772090).
References [1]P. Roncera, F. B. Le, Range of geometrical frustration in lattice spin models, Physical review E. 100(2019) 052150. [2]A. Mohammed, B. Naina, A machine learning approach for the identification of the lattice discrete particle model parameters, Engineering Fracture Mechanics. 197(2018) 160-175. 12 / 16
[3] K. Szajeka, W. Sumelka, Discrete mass-spring structure identification in nonlocal continuum space-fractional model, Eur. Phys. J. Plus,134 (2019)448. [4] Q.C. Sun, F. Jin, G. Q.Wang, et al., On granular elasticity. Scientific Reports, 5(2015)9652. [5]A. Hrennikoff, Solution of problems of elasticity by the framework method, J. Appl. Mech. (1941) A169-A175. [6] D. Greenspan, Supercomputer simulation of cracks and fractures by quasi-molecular dynamics, Journal of Physics & Chemistry of Solids. 50(12) (1989)1245-1249. [7] G. Wang, M. Ostoja-Starzewski, Particle modeling of dynamic fragmentation-I: theoretical considerations, Computational Materials Science. 33(4) (2005)429-442. [8] H. Chen, E. Lin, Y. Jiao, Y. Liu, A generalized 2D non-local lattice spring model for fracture simulation, Comput. Mech. 54(6) (2014)1541-58. [9] S. F. Zhao, G. F. Zhao, Implementation of a high order lattice spring model for elasticity, Int. J. Solids. Struct. 49(18) (2012)2568-81. [10] J. Ai, J. F. Chen, J. M. Rotter, Assessment of rolling resistance models in discrete element simulations, Powder Technology.206(3) (2011)269-282. [11] A. Jagota, S. J. Bennison, Spring-network and finite-element models for elasticity and fracture, Lecture Notes in Physics. 437(1994)186-201. [12] M. Lax, The relation between microscopic and macroscopic theories of elasticity. Solid State Commun. 1(6) (1963) 195-195. [13] P. N. Keating, Effect of invariance requirements on the elastic strain energy of crystals with application to the diamond structure, Phys Rev. 145(2) (1966)637-645. [14] A. Jagota, G. W. Scherer, Viscosities and sintering rates of a two-dimensional granular composite, J Am Ceram Soc. 76(12) (1993)3123-3135. [15] Z.C. Pan, R. J. Ma, D. L. Wang, A. R. Chen, A review of lattice type model in fracture mechanics: theory, applications, and perspectives, Engineering Fracture Mechanics. 190 (2018) 382-409. [16] S. Geer, J. R. Berger, W. J. Parnell, G. G. W. Mustoe, A comparison of discrete element and micromechanical methods for determining the effective elastic properties of geomaterials, Computers and Geotechnics. 87(2017)1-9. [17] M. Y. Zhang, Z. Y. Yang, Z. X. Lu, B. H. Liao, X. F. He, Effective elastic properties and initial yield surfaces of two 3D lattice structures, International Journal of Mechanical Sciences. 138(2018) 146-158. [18] G. F. Zhao, Z. Q. Deng, B. Zhang, Multi-body failure criterion for the four-dimensional lattice spring model, International Journal of Rock Mechanics and Mining Sciences. 123(2019)104-126. [19] G. Wang , A. Al-Ostaz, et al, Hybrid lattice particle modeling of wave propagation induced fracture of solids, Computer Methods in Applied Mechanics and Engineering. 199 (2009) 197-209. [20] G. Wang, A. Al-Ostaz, A. H. D. Cheng, P. R. Mantena, Hybrid lattice particle modeling: Theoretical considerations for a 2D elastic spring network for dynamic fracture simulations, Computational Materials Science. 44(2009)1126-1134. [21] R. K. Choubey, S. Kumar, M. C. Rao, Modeling of fracture parameters for crack propagation in recycled aggregate concrete, Construction and Building Materials. 106(2016) 168-178. [22] C. Jiang, G. F. Zhao, A coupling model of distinct lattice spring model and lattice boltzmann method for hydraulic fracturing, Rock Mechanics and Rock Engineering. 52(10) (2019)3675-3692. [23] M. Braun, M. P. Ariza, New lattice models for dynamic fracture problems of anisotropic materials, Composites Part B-Engineering. 172(2019)760-768. [24] M. Mungule, B. K. Raghuprasad, Meso-scale studies in fracture of concrete: A numerical simulation, Computers and Structures. 89 (2011)912-920. [25] P. R. Budarapu, B. Javvaji, J. Reinoso, A three dimensional adaptive multiscale method for crack growth in Silicon, Theoretical and Applied Fracture Mechanics. 96(2018)576-603. [26] G. Carta, I. S. Jones, et al, Crack propagation induced by thermal shocks in structured media, International Journal of Solids and Structures. (50) (2013) 2725-2736. 13 / 16
[27] A. Mohammadipour, K. Willam, Lattice simulations for evaluating interface fracture of masonry composites, Theoretical and Applied Fracture Mechanics. 82(2016)152-168. [28] H. Chen, E. Lin, Y. Liu, A novel volume-compensated particle method for 2D elasticity and plasticity analysis, International Journal of Solids and Structures. 51(2014)1819-1833. [29] A. A. Gusev, Finite element mapping for spring network representations of the mechanics of solids, Physical Review Letters. 93(3) (2004) 034302-10-15. [30] W. T. Ashurst, W. G. Hoover, Microscopic fracture studies in the two-dimensional triangular lattice, Physical Review B. 14(4) (1976) 1465-1473. [31] D. Frenkel, B. Smit, Understanding molecular simulation: from algorithms to applications, Academic Press, Inc. 1996. [32] G. A. Buxton, C. M. Care, D. J. Cleaver, A lattice spring model of heterogeneous materials with plasticity, Modelling and Simulation in Materials Science and Engineering. 9(6) (2001) 485-497. [33] Y. Yu, W. Wang , H. He, T. Lu, Modeling multiscale evolution of numerous voids in shocked brittle material, Physical Review E. 89(4) (2014) 043309. [34] N. A. Nodargi, F. Caselli, E. Artioli, et al, A mixed tetrahedral element with nodal rotations for large-displacement analysis of inelastic structures, Int. J. Numer. Meth. Engng. 108(2016)722-749. [35] J. Ožbolt, A. Sharma, H. W. Reinhardt, Dynamic fracture of concrete–compact tension specimen, International Journal of Solids and Structures. 48 (10) (2011) 1534-1543.
14 / 16
Conflict of interest statement We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work, there is no professional or other personal interest of any nature or kind in any product, service and/or company that could be construed as influencing the position presented in, or the review of, the manuscript entitled, “Lattice spring model with angle spring and its application in fracture simulation of elastic brittle materials”.
15 / 16
This paper presents a new lattice spring model with rotational degrees of freedom. In this paper, a stiffness matrix with rotational degrees of freedom is established. The model in this paper has high computational accuracy.
16 / 16
S0167-8442(19)30479-3 https://doi.org/10.1016/j.tafmec.2019.102469 TAFMEC 102469
To appear in:
Theoretical and Applied Fracture Mechanics
Received Date: Revised Date: Accepted Date:
28 August 2019 28 December 2019 28 December 2019
Please cite this article as: T. Wang, M. Zhou, Y. Li, Y. Yu, H. He, Lattice spring model with angle spring and its application in fracture simulation of elastic brittle materials, Theoretical and Applied Fracture Mechanics (2020), doi: https://doi.org/10.1016/j.tafmec.2019.102469
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. 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.
© 2020 Elsevier Ltd. All rights reserved.
Lattice spring model with angle spring and its application in fracture simulation of elastic brittle materials Tao Wanga, Mao Zhoua, Yongqiang Lia,b, Yin Yuc, Hongliang Hec a. College of Science, Northeastern University, Shenyang 110819, China b. Key Laboratory of Ministry of Education on Safe Mining of Deep Metal Mines, Northeastern University, Shenyang, Liaoning 110819, China. c. National Key Laboratory of Shock Wave and Detonation Physics, Institute of Fluid Physics, CAEP, 621999 Mianyang, China
ABSTRACT Lattice Spring Model with Angle Spring (LSMA) is a novel type of lattice model. Compared with the general lattice spring model, LSMA adds the rotational and tangential degrees of freedom, whose spring parameters are obtained by finite element parameter mapping. The traditional lattice spring model only consider the normal interaction between two points (the Poisson’s ratio is fixed), which limits the application of the lattice spring model in a wider range of materials. Some scholars proposed to add shear spring into the traditional model to solve the problem of fixed Poisson’s ratio, but the rotation invariance could not be maintained. The main reason is that the shear spring can not distinguish the difference of tangential velocity (or displacement) between the two particles due to the common rotation or shear. Therefore, the rotation of the whole rigid body may cause incorrect generation of additional strain energy inside the shear spring. LSMA model is able to maintain rotation invariance and reproduce different Poisson’s ratio because of the spring parameters with rotational degrees of freedom and the rotation effect of lattice points. In this paper, the finite element stiffness matrix with rotational freedom is derived by using LSMA as the basic mechanical model, and the parameter mapping theory of spring stiffness coefficient is introduced. By means of numerical simulation, the deflection and rotation Angle of cantilever beam, shear wave and compressional wave velocity of stress wave and the simulation of compact tensile test are checked and calculated respectively, and the correctness of the model is verified. Keywords: Angle spring, Lattice spring model, Elastic brittle material, Poisson's ratio
1. Introduction Rock is a typical elastic brittle material, and the dynamic fracture is the core problem of earthquake, landslide, rock blasting and other mechanical behaviors. Therefore, it is of great theoretical and practical significance to simulate the dynamic fracture effectively. There are two reasons, first, on the level of mechanics, the dynamic fracture mechanism of rock is not very clear, there is not a perfect mechanical model can accurately describe the dynamic fracture mechanism of rock. Second, at the algorithm level, how to simulate dynamic crack generation, propagation and instability is still a very difficult problem. There are two main types of numerical simulation methods for rock dynamic fracture: continuous medium method and discrete medium method. Although the theory of dynamic fracture mechanics of continuous media has made good progress, the theoretical system still has many difficulties in simulating dynamic fracture of rock. Traditional continuum mechanics is a theoretical system based on the continuous assumption of field, regardless of the microstructure of the material. When simulating the fracture of material, the external fracture criterion should be used. The selection of dynamic fracture criteria is a very difficult problem. In addition, when there is crack propagation, The processing of finite element mesh is another difficult problem of the theoretical system. Compared with the continuous medium mechanics method, the discrete medium method has a unique advantage
Corresponding author at: College of Sciences, Northeastern University, Shenyang 110819, China. Tel.: +86- 24-
83672653. E-mail address: [email protected] (Yongqiang Li). 1 / 16
in simulating the dynamic fracture of rock. By using a discrete system instead of a continuous medium, the 3-D material failure problem is transformed into a one-dimensional problem. The problems of continuum fracture criterion and the transition from continuous medium to discontinuous medium are avoided, which greatly simplifies the fracture simulation. At present, the mainstream discrete medium models mainly include lattice model [1,2], spring model [3] and particle flow model [4]. Lattice model as a discrete medium method has a strong advantage in simulating dynamic fracture. Lattice model was first proposed by Hrennikoff [5] to solve classical elastic mechanics problems. In this method, the continuum material is equivalent to a discrete structural system, the generation and evolution mechanism of micro-damage in the material is regarded as the failure of a single bar element or spring element in the discrete structural system, the model elements adopt simple constitutive relation and strength criterion. In this way, the problem of determining the constitutive relation of materials and the selection of dynamic fracture criterion is avoided, which provides a simple and effective method for the simulation of dynamic fracture of materials. Early in the lattice model development [6, 7], connections in discrete structures are made of springs, due to the spring model only consider normal interactions between two points, so it limits the material Poisson's ratio (1/4 for plane strain, 1/3 for plane stress), which restricts the application of the model. In order to break through this limitation, some scholars have added additional shear springs into the traditional model [8-10]. The introduction of shear springs can simulate any Poisson’s ratio. However, introducing shear springs into the classical LSM may not maintain rotational invariance [11-13]. Jagota and Scherer[14] believe that the main reason is that the shear spring cannot distinguish the difference between the tangential velocity (or displacement) of the two particles due to their joint rotation or shear, as shown in Fig.1. As a result, global rigid body rotation may incorrectly cause additional strain energy to be generated within the shear spring.
A
ks
B
A
rotation
ks
B
shear
Fig.1 Schematic showing incapability of shear spring to distinguish differences in tangential displacements of two particles owing to rigid rotation or shear. The literature [16,17] used the beam element to replace the spring in the particle model, thus the particle turning effect is taken into account, so that the lattice model could simulate an arbitrary Poisson’s ratio. However, the calculation of rotation Angle is indirectly calculated by lattice point displacement, which requires a large amount of calculation. Zhao et al. [18] considered the normal deformation and shear deformation effects of the spring in the spring model, and solved the fixed Poisson's ratio problem of the lattice model. In this method, the shear deformation directly calculated by the particle displacement still included the rigid body rotation. The literature[18] used local strain method to eliminate the influence of rigid body rotation. In order to consider the influence of bond Angle, Wang [19,20] proposed the concept of angle spring, which can be used to reproduce variable Poisson’s ratio, so that the lattice model has a great advantage in simulating dynamic fracture. In the traditional grid method, the parameter calibration of the bond is related to the specific case. There is no uniform micro-macro parameter correspondence, and the lattices of different topologies have different micro-macro parameter relationships. These all make the parameter calibration problem of the three-dimensional lattice model very complicated. In addition, the traditional lattice model cannot reproduce the microscopic structure of rock. In recent years, the lattice model has been constantly improved in theory, and it has been widely used in the process of dynamic fracture simulation, such as [21-28]. 2 / 16
uA
vB
kn
A
B
uB
kr
ks
A
B
A
B
A
vA
normal
B
shear
angle
Fig.2 Lattice spring model with angular spring In this paper, the lattice spring model with Angle spring is established, as shown in Fig.2. In the figure, u and v are translational degrees of freedom along the normal and tangential directions, and φ is a rotational degree of freedom, with these three degrees of freedom coupling on the lattice. The quantitative mapping between the stiffness matrix of the finite element and the stiffness coefficient of the lattice spring model is realized by using the finite element of three joints and three degrees of freedom with rotational degrees of freedom as the medium and the model mesh shared by the lattice spring model and the finite element model. In this paper, the finite element stiffness matrix with rotational degrees of freedom is derived, the parameter mapping theory of spring stiffness coefficient is explained, and the fracture criterion of spring is introduced. To verify the correctness of the model, the deflection and rotation Angle of the cantilever beam, the shear and longitudinal wave velocities of the stress wave and the simulation of the compact tensile test are checked by means of numerical simulation.
2 3-node finite element stiffness matrix with rotational degrees of freedom In this section, the rotation angles of six-node triangular elements are introduced to generate interpolation functions of triangular elements with rotational degrees of freedom. The six-node linear strain triangular element is shown in Fig.3. Nodes 1, 2 and 3 are intermediate nodes of each side. The displacement of each node along the x and y axes is represented by u and v respectively. αjm, αmi and αij represent the angle between the outer normal of each edge and the x-axis. m um , vm n jm nmi
mi
m(U m ,n , U m , s )
jm
n jm
1 u1 , v1
jm
2 u2 , v2
1(U1,n ,U1,s ) j u j , v j
y
x
i ui , vi
3 u3 , v3
ij nij
Fig. 3 Six-node linear strain triangular element.
x1
y1
j(U j ,n ,U j ,s )
Fig. 4 The local coordinate system
For the six-node line strain triangle element shown in Fig.3, rotational degrees of freedom are introduced on this element and eliminate the middle node 1, 2 and 3 of the edge. According to the interpolation method adopted by Allman, the tangential displacement Us along each side of the triangle varies linearly, while the normal displacement Un along the side changes quadratically. Therefore, in order to establish the local coordinate system along tangential direction and normal direction, coordinate transformation is required before interpolation. The specific transformation method is as follows: the local coordinate yl axis direction obtained after transformation is consistent with the external normal direction of each side, and the local coordinate x1 axis direction is consistent with the tangent direction of each side. The coordinate transformation is shown in Fig. 3 and Fig. 4, the displacement relation of the normal displacement and tangential displacement transformation of each side in the triangular element under the rectangular coordinate system is 3 / 16
u cos v sin
sin U n , jm , ij , mi cos U s
(1)
After the coordinate transformation by Eq.(1), the corner nodes are used to represent the middle nodes of the edge and the middle nodes are eliminated, interpolating in local coordinates respectively. Since the tangential displacement varies linearly along the edge, the tangential displacement of the middle node of each edge is
U 1, s
U j ,s U m,s 2
, U 2, s
U i ,s U m,s 2
U i,s U j ,s
, U 3, s
(2)
2
The normal displacement is maintained as a quadratic interpolation along the edge, Let the interpolation function be a quadratic parabola.
U n = a1 x12 + a2 x1 + a3
(3)
The rotation (slope) on that side is φ, so
=
dU n = 2a1 x1 + a2 dx1
(4) um
m
m
vm
uj j
ui
i
i
j
vj
vi
Fig. 5. Three-node triangular element with nine degrees of freedom Let nodes i, j, m have rotation parameters φi, φj and φm, as shown in Fig.5. The second derivative of Un can also be represented by the rotation parameters of nodes, taking edge j-m in Fig.4 as an example
j - m d 2U n d = = 2 a = 1 dx1 l jm dx12
(5)
where ljm is the length of the side j-m, there are
a1 =
j - m
(6)
2l jm
Substituting the local coordinates (0,Uj,n) and (ljm,Umn) of node j and node m into Eq.(3) respectively, there are
U j,n = a3 , U m,n = a1l 2jm + a2l jm + a3
(7)
According to Eq.(6) and Eq.(7), the coefficients a1, a2 and a3 of Eq.(3) can be solved
a1 =
j - m 2l jm
, a2 =
U m,n - U j ,n l jm
-
j - m 2
, a3 = U j ,n 4 / 16
(8)
Substituting Eq. (8) into Eq. (3) and taking x1=ljm/2, we can obtain the normal displacement of the middle node 1 of edge j-m as U j ,n + U m,n
U1,n =
2
+
m - j 8
(9)
l jm
Similarly, the normal displacements of the intermediate nodes 2 and 3 of edges m-i and i-j can be obtained as follows
U 2, n
U m,n U i ,n 2
i m 6
lmi , U 3,n
U j ,n U i ,n 2
j i 6
lij
(10)
where lmi and lij are the side lengths of edges m-i and i-j respectively. The interpolation function of a six-node triangular element is 3 u N u N u N u N k uk i i j j m m k 1 3 v N v N v N v N v i i j j m m k k k 1
(11)
and
l jm sin jm x j xm ci ,i j m l jm cos jm ym y j bi
(12)
Based on Eq. (1), Eq. (2), Eq. (9) , Eq. (10) , Eq. (11) and Eq. (12), the interpolation function of triangular elements with rotational degrees of freedom can be expressed as
2 u Li ui L j u j Lmum 3 i j bm Li L j j m bi L j Lm m i b j Lm Li v L v L v L v 2 c L L c L L c L L i i j j m m i j m i j j m i j m m i j m i 3
(13)
where, Li, Lj and Lm are area coordinates. The three nodes of the trilateral midpoint are eliminated. The three-node triangular element with nine degrees of freedom, as shown in Fig.5, is obtained by introducing rotational degrees of freedom into the vertex of the triangular element. From Eq.(13) can be deduced
u
i
vi uj vj um vm u1 v1 u2 v2 u3 v3
T
(14)
Tui vi i uj vj j um vm m
T
where T is the conversion matrix of 12×9.
5 / 16
0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 T 12 0 0 0 bm 6 1 2 0 bm 6 0 0 1 2 am 6 0 0 0 0 1 2 am 6 0 0 0 12 0 bi 6 1 2 0 bi 6 0 0 0 0 1 2 ai 6 0 1 2 ai 6 0 0 12 0 bj 6 1 2 0 b j 6 0 0 1 2 a j 6 0 0 0 0 1 2 a j 6
(15)
e
Based on Eq.(14), relational expression of K 99 (stiffness matrix of triangular element with three nodes and e
nine degrees of freedom) and K 1212 (stiffness matrix of triangular element with three nodes and twelve degrees of freedom) can be expressed as e K 9e9 TT K12 12 T
where K 1212 t e
(16)
B DBdxdy , T
A
B Bi
Bj
Bm
B1 B 2
B3
3 Spring parameter mapping theory Whether the lattice point spring model can accurately evolve the response of material under external load depends on the selection of the mass of lattice point and the elastic coefficient of spring, so that the spring network is just consistent with the elastic property of the actual sample. According to the cell size and material density, the lattice mass can be determined by the method of equivalent mass. The core problem of lattice spring model modeling has been transformed into the mathematical mapping relationship between the elastic constitutive relation of actual material and the elastic coefficient of spring network. Gusev [29] reported a method of mapping finite element parameters to discrete element models. When material satisfies elastic deformation, the relationship between node force and node displacement in finite element models can be established by finite element theory.
F Kδ
is global stiffness matrix,
Where F Fi is node force vector, K k ij T
(17)
δ i is node displacement T
vector, i, j 1, 2,..., n , n is degree of freedom. According to Eq. (17), the force of a node i can be rewritten as n
Fi kij j ki11 ki 2 2 kii i kin n
(18)
j 1
The element kii on the principal diagonal of the global stiffness matrix is equal to the opposite number of all elements` sum except the element of principal diagonal in the same row[30], i.e., n
k ii k im
(19)
m 1 m i
6 / 16
Based on Eq. (19), Eq. (18) can be rewritten as n
Fi kim m i
(20)
m 1 mi
As can be seen from Eq.(20), The force of a node is equal to the sum of the product of the adjacent spring parameters and the displacement difference of the node. According to this principle, the stiffness coefficients in the finite element are mapped to the elastic coefficients of the springs in the lattice spring model, and the values of the spring network in the lattice spring model are determined. This parameter mapping method has a rigorous mathematical logic, to ensure that the model can quantitatively represent the mechanical properties of real materials in the stage of elastic deformation, which can solve the difficulty of lacking strict demonstration in the parameter selection of discrete element method model. At the same time, the model does not depend on the background grid when calculating, which avoids the difficulty of numerical calculation caused by the mesh distortion when solving the extremely large deformation and fracture problems by the finite element method. Lattice spring model can impose arbitrary boundary, which effectively solves the problem that mesh-less method can’t impose displacement boundary. The given size of the basic element model is small enough and the number of elements is large enough, the lattice spring model can characterize the dynamic response of the sample from both mesoscopic and macroscopic levels, providing the necessary physical knowledge for the damage evolution of materials.
4. Criterion for spring fracture of elastic brittle materials In order to use the LSM model to simulate the fracture process of the material, a failure criterion needs to be defined to identify the key elements to be removed from the model. Most failure criteria used so far in LSM have been classified into two categories: critical stress (force) criterion and critical strain (displacement) criterion. The critical strain criterion is adopted for calculation in this paper. Because rock and concrete are typical elastic-brittle materials, the constitutive relationship between force F and strain of spring in fracture simulation can be shown in Fig. 6.
F
o
t
Fig.6 Constitutive relationship between force F and strain of spring
t is the critical tensile strain of the spring. Before the spring breaks, the spring is linear elastic. When the principal strain of the spring exceeds t value, the spring breaks. Pure compression (such as hydrostatic compression of uniform media) does not cause material damage, so the spring compression failure is not considered in the calculation. When the fracture criterion is satisfied, the spring will be broken permanently, and the tension and shear resistance between the two lattices will be lost, which is equivalent to forming a small crack in the sample.
5. Kinematics integral In the process of model calculation, the translational and rotational accelerations of each lattice in the unit time step can be obtained, and then the displacement, velocity, angular velocity and angular displacement of each 7 / 16
particle in the next time step can be obtained by acceleration-time integral. The basic principle of the motion equation integration algorithm in lattice spring model and the molecular dynamics is to perform Taylor expansion on the motion equation, i.e.[31],
r (t t ) r (t ) t
r t 2 2 r t 3 3 r ... t 2! 2t 3! 3t
(21)
where r , t , t are the coordinates, time and time step respectively. Generally, when the order of the integration algorithm is expanded very high, the longer the time step of each iteration can be obtained, but at the same time, the more motion information need to be stored in integral operation. Verlet algorithm and its substitution are commonly used low-order algorithms. Although it is necessary to shorten the time step to ensure sufficient accuracy, the particle trajectories calculated by Verlet algorithm are usually time reversible, and the total energy of the system rarely drifts in long-term operation. This paper chooses an alternative method of Verlet algorithm-Beeman algorithm,
4 f (t ) f (t t ) 2 ( ) ( ) ( ) r t t r t v t t t 6m v(t t ) v(t ) 2 f (t t ) 5 f (t ) f (t t ) t 6m
(22)
Where f(t) is the resultant force of the particle at time t, m is the mass of the particle, v is the velocity of the particle. The particle trajectory calculated by Beeman algorithm is the same as Verlet algorithm, but it is more accurate in calculating particle velocity and the total energy conservation property of the model is better.
6. Verification 6.1 Deflection and angle deformation calculation of rectangular section cantilever beam In order to verify the correctness of the model, a rectangular cantilever beam structure is selected. The geometry and boundary conditions of the problem are described in Fig.7. The left end of the beam is fixed in all directions. A force of 10 N is applied at the right end of the beam. The elastic constants of the material are: Young’s modulus E=200GPa and Poisson’s ratio μ= 0.22. Based on the elasticity theory, the maximum deflection and rotation angles are obtained respectively，
u y max
Fl 3 , 3EI
max
Fl 2 2 EI
(23)
Where I bh3 /12 , the calculation results of the maximum deflection and the maximum rotation angle in the Y direction with different number of lattices are shown in Table 1. Suppose the number of x-axis lattice points is mx; y is my y
x
h b
l
F
Fig.7 cantilever beam structure diagram (l=1.0m,b=h=0.01m)
8 / 16
Table 1 Calculation results of maximum deflection and rotation angle of cantilever beam under different number of nuclear power plants max (rad)
uymax (m) Lattice (mx×my)
LSMA Calculate values
without angle spring Error(%)
Calculate values
Error(%)
LSMA Calculate values
without angle spring Error(%)
Calculate values
20×2
0.015
25.0
0.0018
91.0
0.025
16.7
-
100×2
0.017
15.0
0.0141
29.5
0.026
13.3
-
200×2
0.018
10.0
0.0179
10.5
0.028
6.7
-
200×4
0.019
5.0
0.0181
9.5
0.029
3.3
-
400×4
0.020
0.0
0.0194
3.0
0.030
0.0
-
Theoretical values
0.020
0.030
In order to test the accuracy of LSMA, the comparison results among the calculated values of LSMA, the calculated values of the model without angle spring and the theoretical calculated values are given in Table 1. The relative error in Table 1 refers to the error between the numerical simulation solution and the theoretical solution. As can be seen from the comparison results in Table 1, the calculation accuracy of the LSMA is much higher than the model without angle spring. Because the spring stiffness coefficient in the lattice spring method without angle spring is calculated by the 3-node 6-degree-of-freedom triangular finite element. The 3-node triangular finite element is a constant strain element, which has good reliability and low calculation accuracy as a coordination element. However, the spring stiffness coefficient in the LSMA is calculated by the 3-node, 9-degree-of-freedom uncoordinated triangular finite element, which is obtained from the original Allman element, keeping the accuracy of the original Allman element. It can also be seen from Table 1 that the calculation accuracy of lattice spring method increases greatly with the increase of lattice points. 6.2 Large deflection bending deformation calculation of slender cantilever beam The lattice points in the traditional lattice-spring model are only subject to normal forces (the force is only related to the relative distance between two lattice points), which is called the "Hookean model"[32]. In Hookean model, there is no resistance to shear deformation between two lattice points. As a result, Hookean model only has fixed Poisson's ratio, which cannot simulate various Poisson's ratio materials. To break this limitation, several improved models have been proposed, the most popular model of which is the "Born model"[32]. Born model adds shear forces between two lattice points, which can control the overall Poisson's ratio by adjusting the tangential spring stiffness coefficient. One disadvantage of the Born model is that it does not have the rotational conservation of strain energy[33]. This shortcoming is not prominent in the study of small deformation, but it will have an impact on large deformation phenomena. The LSMA proposed in this paper by introducing the angle spring, the Born model can obtain the rotational conservation of strain energy automatically[33]. To verify LSMA's ability to calculate large deformation, the slender cantilever beam with rectangular section is selected as shown in Fig.7. The cantilever beam, of length l=10m, is characterized by a rectangular cross-section, of thickness h=0.1m and of width b=1m. It is assumed to be constituted of a linear elastic material, exhibiting Young’s modulus E =1.2×106 and Poisson’s ratio μ=0. The loading history consists in the application of a uniform follower load q = q0EI/l3 , where I is the cross-section moment of inertia. q0 is the load coefficient, increasing from zero up to the value q0max = 40.
9 / 16
1.0
q0 15.0
0.5
q0 10.0
q0 25 .0
q0 30.0
q0 5.0
q0 35 .0
Tip deflections(/l)
q0 20.0
0.0
utip,x
-0.5
utip,y
-1.0
utip,y theoretical
utip,x theoretical
q0 40.0 y
-1.5
x
0
10
20
30
40
Applied load coefficient q0
(a)
(b)
Fig.8. Clamped cantilever subjected to uniformly distributed follower load: (a) reference and deformed configurations at eight uniformly spaced load levels and (b) load–displacement curve.
Fig.8(a) shows the large-displacement deformed configuration at eight uniformly spaced load levels, where
usum u x2 u y2 , ux and uy is the displacement in the x and y directions. Fig.8(b) shows the load-displacement curve corresponding to x-deflection utip,x and y-deflection utip,y are reported. The comparison of the numerical results for the considered quite slender structure with the reference solution reported in [34] proves the effectiveness of the present formulation. 6.3 Wave velocity test In solid materials, both shear and longitudinal waves can be effectively propagated. The propagation velocities e
of shear and longitudinal waves depend only on the properties of materials, which can be represented by CL and
CRe in the one-dimensional elastic strain wave, i.e.,
CLe
E 1 E e , CR 2 1 1 1 2 1
(24)
where, is the density of the material. The theoretical values of CL and CR can be calculated by Eq. (24). The calculated values of shear and longitudinal waves velocity calculated by the LSMA and without angle spring model in section 6.1 example are shown in Table 2. e
e
Table 2 Shear and longitudinal wave velocities of different materials Shear waves velocity (m/s)
Longitudinal waves velocity(m/s)
without angle without angle LSMA LSMA ρ spring spring (kg/m3) Theoretical Theoretical values values Calculated Error Calculated Error Calculated Error Calculated Error values (%) values (%) values values (%) (%)
E (GPa)
μ
70
0.33
2700
3121.95
3180.09
1.86
3201.25
2.54
6197.82
6167.34
0.49
6153.82
0.71
120
0.32
7800
2414.02
2390.63
0.97
2380.95
1.37
4692.02
4625.83
1.41
4601.46
1.93
150
0.22
8010
2770.35
2690.89
2.87
2655.10
4.16
4623.84
4626.3
0.05
4632.16
0.18
210
0.25
7800
3281.65
3207.89
2.25
3181.89
3.04
5683.99
5626.38
1.01
5572.58
1.96
10 / 16
The calculation results show that the maximum error between the LSMA model and the without angle spring model and the theoretical values is 2.87% and 4.16% respectively for the shear and longitudinal wave velocities of materials, indicating the accuracy of the LSMA model. 6.4 Compact Tensile Test Simulation Dynamic crack propagation is closely related to loading rate. In order to verify whether LSMA model can reproduce the effect of loading rate on dynamic crack growth, this model is used to simulate compact tensile test. Geometry of the compact tension specimen is length-width = 0.2m×0.2m. The length and width of the notches are 0.064 m and 0.018 m respectively. The left end of the notch is fixed by the steel sleeve, and the right end is loaded by the displacement control of the same steel sleeve. Similar experiments have been done with concrete in literature [35]. The results show that different failure modes appear with different loading rates. The cracking modes of specimens under different loading rates are simulated by this example. The whole specimen is simulated with 6468 springs. The geometric, LSMA model mesh and boundary conditions of the specimen are shown in Fig.9. In order to more accurately reproduce the contact between the specimen and the loader during the actual loading process, a steel plate with relatively high stiffness is installed in the loading area. The displacement control load is applied on the steel plate, and the contact between the steel plate and the specimen is controlled by the bonding interface unit.
v
0.2m
0.2m
Fig. 9 Geometry of the compact tension specimen, LSMA model mesh, loading (displacement control) and boundary conditions. According to literature[35], the numerical simulation is carried out for normal strength concrete with the following macroscopic properties: Young’s modulus E = 36 GPa, Poisson’s ratio μ = 0.18, the critical tensile strain t= 0.116×10-3. The specimen is analyzed assuming dynamic loading condition with the following displacement rates v=0.304m/s, 1.375m/s and 3.318m/s. The experimental results of the three loading rates are shown in Fig. 10(a), the simulation results are shown in Fig.10(b) by LSMA, and the simulation results are shown in Fig.10(c) without angle spring model. As can be seen from Fig.10, the simulation results by LSMA are basically consistent with the experimental results. When the loading rate is small (i.e. 0.304m/s), the crack propagates along the central line, and the failure mode is similar to that of quasi-static loading; when the loading rate increases to 1.375m/s, the crack propagates towards the loading end gradually; when the loading rate continues to increase to 3.318m/s, the crack propagates appeared bifurcation phenomenon. With the increase of loading rate, the fracture mode of specimens gradually transits from mode I to mixed-mode.
11 / 16
v=0.304 m/s
v=1.375 m/s
v=3.318 m/s
(a) experimental result[32]
(b) LSMA
(c) without angle spring
Fig.10 Comparison of experimental results and simulation results of crack path under different loading velocities
7 Conclusions In this paper, the LSMA model uses the spring parameters with rotational degrees of freedom and considers the angular effect of the lattice points. It can keep the rotation invariance and reproduce different Poisson's ratios. Based on the finite element model with three nodes and nine degrees of freedom including rotational degrees of freedom, by means of the mesh model shared by lattice spring model and finite element model, the stiffness matrix of finite element is mapped to the spring stiffness coefficient of lattice spring model, so that the selection of spring stiffness coefficient can be deduced mathematically. In order to verify the correctness of the model, the deflection and rotation of the cantilever beam, the shear and longitudinal wave velocities of stress wave and the compact tensile test simulation are checked by numerical simulation. The lattice spring model with angle spring is more accurate than the model without angle spring.
Acknowledgement This work was supported by the National Natural Science Foundation of China (No. 11772090).
References [1]P. Roncera, F. B. Le, Range of geometrical frustration in lattice spin models, Physical review E. 100(2019) 052150. [2]A. Mohammed, B. Naina, A machine learning approach for the identification of the lattice discrete particle model parameters, Engineering Fracture Mechanics. 197(2018) 160-175. 12 / 16
[3] K. Szajeka, W. Sumelka, Discrete mass-spring structure identification in nonlocal continuum space-fractional model, Eur. Phys. J. Plus,134 (2019)448. [4] Q.C. Sun, F. Jin, G. Q.Wang, et al., On granular elasticity. Scientific Reports, 5(2015)9652. [5]A. Hrennikoff, Solution of problems of elasticity by the framework method, J. Appl. Mech. (1941) A169-A175. [6] D. Greenspan, Supercomputer simulation of cracks and fractures by quasi-molecular dynamics, Journal of Physics & Chemistry of Solids. 50(12) (1989)1245-1249. [7] G. Wang, M. Ostoja-Starzewski, Particle modeling of dynamic fragmentation-I: theoretical considerations, Computational Materials Science. 33(4) (2005)429-442. [8] H. Chen, E. Lin, Y. Jiao, Y. Liu, A generalized 2D non-local lattice spring model for fracture simulation, Comput. Mech. 54(6) (2014)1541-58. [9] S. F. Zhao, G. F. Zhao, Implementation of a high order lattice spring model for elasticity, Int. J. Solids. Struct. 49(18) (2012)2568-81. [10] J. Ai, J. F. Chen, J. M. Rotter, Assessment of rolling resistance models in discrete element simulations, Powder Technology.206(3) (2011)269-282. [11] A. Jagota, S. J. Bennison, Spring-network and finite-element models for elasticity and fracture, Lecture Notes in Physics. 437(1994)186-201. [12] M. Lax, The relation between microscopic and macroscopic theories of elasticity. Solid State Commun. 1(6) (1963) 195-195. [13] P. N. Keating, Effect of invariance requirements on the elastic strain energy of crystals with application to the diamond structure, Phys Rev. 145(2) (1966)637-645. [14] A. Jagota, G. W. Scherer, Viscosities and sintering rates of a two-dimensional granular composite, J Am Ceram Soc. 76(12) (1993)3123-3135. [15] Z.C. Pan, R. J. Ma, D. L. Wang, A. R. Chen, A review of lattice type model in fracture mechanics: theory, applications, and perspectives, Engineering Fracture Mechanics. 190 (2018) 382-409. [16] S. Geer, J. R. Berger, W. J. Parnell, G. G. W. Mustoe, A comparison of discrete element and micromechanical methods for determining the effective elastic properties of geomaterials, Computers and Geotechnics. 87(2017)1-9. [17] M. Y. Zhang, Z. Y. Yang, Z. X. Lu, B. H. Liao, X. F. He, Effective elastic properties and initial yield surfaces of two 3D lattice structures, International Journal of Mechanical Sciences. 138(2018) 146-158. [18] G. F. Zhao, Z. Q. Deng, B. Zhang, Multi-body failure criterion for the four-dimensional lattice spring model, International Journal of Rock Mechanics and Mining Sciences. 123(2019)104-126. [19] G. Wang , A. Al-Ostaz, et al, Hybrid lattice particle modeling of wave propagation induced fracture of solids, Computer Methods in Applied Mechanics and Engineering. 199 (2009) 197-209. [20] G. Wang, A. Al-Ostaz, A. H. D. Cheng, P. R. Mantena, Hybrid lattice particle modeling: Theoretical considerations for a 2D elastic spring network for dynamic fracture simulations, Computational Materials Science. 44(2009)1126-1134. [21] R. K. Choubey, S. Kumar, M. C. Rao, Modeling of fracture parameters for crack propagation in recycled aggregate concrete, Construction and Building Materials. 106(2016) 168-178. [22] C. Jiang, G. F. Zhao, A coupling model of distinct lattice spring model and lattice boltzmann method for hydraulic fracturing, Rock Mechanics and Rock Engineering. 52(10) (2019)3675-3692. [23] M. Braun, M. P. Ariza, New lattice models for dynamic fracture problems of anisotropic materials, Composites Part B-Engineering. 172(2019)760-768. [24] M. Mungule, B. K. Raghuprasad, Meso-scale studies in fracture of concrete: A numerical simulation, Computers and Structures. 89 (2011)912-920. [25] P. R. Budarapu, B. Javvaji, J. Reinoso, A three dimensional adaptive multiscale method for crack growth in Silicon, Theoretical and Applied Fracture Mechanics. 96(2018)576-603. [26] G. Carta, I. S. Jones, et al, Crack propagation induced by thermal shocks in structured media, International Journal of Solids and Structures. (50) (2013) 2725-2736. 13 / 16
[27] A. Mohammadipour, K. Willam, Lattice simulations for evaluating interface fracture of masonry composites, Theoretical and Applied Fracture Mechanics. 82(2016)152-168. [28] H. Chen, E. Lin, Y. Liu, A novel volume-compensated particle method for 2D elasticity and plasticity analysis, International Journal of Solids and Structures. 51(2014)1819-1833. [29] A. A. Gusev, Finite element mapping for spring network representations of the mechanics of solids, Physical Review Letters. 93(3) (2004) 034302-10-15. [30] W. T. Ashurst, W. G. Hoover, Microscopic fracture studies in the two-dimensional triangular lattice, Physical Review B. 14(4) (1976) 1465-1473. [31] D. Frenkel, B. Smit, Understanding molecular simulation: from algorithms to applications, Academic Press, Inc. 1996. [32] G. A. Buxton, C. M. Care, D. J. Cleaver, A lattice spring model of heterogeneous materials with plasticity, Modelling and Simulation in Materials Science and Engineering. 9(6) (2001) 485-497. [33] Y. Yu, W. Wang , H. He, T. Lu, Modeling multiscale evolution of numerous voids in shocked brittle material, Physical Review E. 89(4) (2014) 043309. [34] N. A. Nodargi, F. Caselli, E. Artioli, et al, A mixed tetrahedral element with nodal rotations for large-displacement analysis of inelastic structures, Int. J. Numer. Meth. Engng. 108(2016)722-749. [35] J. Ožbolt, A. Sharma, H. W. Reinhardt, Dynamic fracture of concrete–compact tension specimen, International Journal of Solids and Structures. 48 (10) (2011) 1534-1543.
14 / 16
Conflict of interest statement We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work, there is no professional or other personal interest of any nature or kind in any product, service and/or company that could be construed as influencing the position presented in, or the review of, the manuscript entitled, “Lattice spring model with angle spring and its application in fracture simulation of elastic brittle materials”.
15 / 16
This paper presents a new lattice spring model with rotational degrees of freedom. In this paper, a stiffness matrix with rotational degrees of freedom is established. The model in this paper has high computational accuracy.
16 / 16