- Email: [email protected]
An ordinary state-based peridynamic modeling for dynamic fracture of laminated glass under low-velocity impact
Journal Pre-proofs An ordinary state-based peridynamic modeling for dynamic fracture of laminated glass under low-velocity impact Liwei Wu, Lei Wang, Dan Huang, Yepeng Xu PII: DOI: Reference:
S0263-8223(19)33108-3 https://doi.org/10.1016/j.compstruct.2019.111722 COST 111722
To appear in:
Composite Structures
Received Date: Revised Date: Accepted Date:
18 August 2019 2 October 2019 21 November 2019
Please cite this article as: Wu, L., Wang, L., Huang, D., Xu, Y., An ordinary state-based peridynamic modeling for dynamic fracture of laminated glass under low-velocity impact, Composite Structures (2019), doi: https://doi.org/ 10.1016/j.compstruct.2019.111722
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.
© 2019 Published by Elsevier Ltd.
An ordinary state-based peridynamic modeling for dynamic fracture of laminated glass under low-velocity impact Liwei Wu1, Lei Wang1, Dan Huang1 *, Yepeng Xu1 College of Mechanics and Materials, Hohai University, Nanjing, 211100, China.
Abstract: Laminated glass, as a typical kind of sandwich-structure composite, is widely used in various fields such as safety, vehicle and transportation engineering. Laminated glass mainly suffers from dynamic/impact loadings, and glass-ply cracking is the main failure pattern of laminated glass. This paper presents a non-local ordinary state-based peridynamic modeling and numerical approach for simulating the dynamic failure process of sandwiched laminated glass under low-velocity impact loading. The mechanical behavior of the PVB layer is simulated by reformulating classical visco-elastic model under the framework of ordinary state-based peridynamic theory, and the adhesion between the glass and PVB interlayer is described by using a penalty-based method. Fracture patterns of a laminated glass plate under drop-weight loading were investigated, and the numerical results compared well with experimental observations. Furthermore, a series of numerical simulations were conducted to in-depth analyze the effect of the thickness of the PVB interlayer and the glass layers on the fracture mode, initial locations and crack propagation speed of the laminated glass when subjected to impact. Keywords: Laminated glass; impact failure; peridynamics; non-local; visco-elastic
*Corresponding
Author: Dan Huang. Email: [email protected]
1. Introduction As a type of widely used sandwich composite structure, laminated glass not only owns the capacity of shielding for wind, rain, heat and sound, but also plays a significant role as a safety component in infrastructure as well as military fields. Typical laminated glass usually consists of two glass plies adhered by a soft interlayer, preferably polyvinyl butyral (PVB). Laminated glass suffers from low-velocity impact in various cases, such as impact of a small and low-velocity stone thrown by the wheel in automotive applications, and impact of wind-borne debris in architectural applications. Effectively declaring the failure mechanism of laminated glass when subjected to low-velocity impact is of important value for improving the resistance as well as safety of laminated glass. This motivates the present work to study the characteristics of dynamic fracture of the laminated glass under low-velocity impact. In the past decades, numerous experimental and numerical studies have been conducted on the mechanical behavior of laminated glass, particularly on the failure mechanism and fracture modeling. For example, Muralidhar et al. [1] and Xu et al. [2] experimentally investigated the debonding of glass-polymer interface under quasi-static and dynamic loading. Chen et al. [3, 4] experimentally studied the cracking behavior of the glass layer in a laminated glass panel. As to the modeling for failure of laminated glass, a series of numerical methods have been proposed in which the most commonly used is the element deletion or erosion contact technologies [5-7] based on finite element method, and incorporated continuum damage mechanics (CDM) method [8, 9]. Zang et al. [10] developed a three-dimensional discrete element
model to simulate the dynamic fracture of laminated glass, and Zhang and Lei [11] proposed a combined discrete element/finite element method (DEM/FEMs) based on a penalty function method to improve the calculation efficiency. Furthermore, Xu et al. [12] simulated performed a parametric study on the fracture of laminated glass plate subjected to low-velocity impact with the extended finite element method (XFEM), and Chen et al. [13] analyzed the crack propagation in the glass layer of a laminated glass plate under drop-weight impact by using the extrinsic cohesive zone model (CZM). In recent years, another alternative titled peridynamics (PD) [14] exhibits significant potential in handling fracture-concerned problems, in particular, simulating impact failure. The recently proposed non-local peridynamic theory [14] is a reformulation of the classical continuum mechanics particularly for handling problems with discontinuities and long-range forces. As opposed to conventional continuum mechanics theory, an integral-typed equation of motion without spatial derivatives is employed to describe the mechanical behavior of solids. The primary advantage of the peridynamic theory lies in its validity in the presence of discontinuities and its ability to model discontinuities naturally and spontaneously without any extra technical remedies employed in conventional continuous methods [15]. Owing to its unique advantages in naturally incorporating discontinuities within a single continuum framework, the peridynamic theory has been successfully applied to study various fracture problems, including quasi-static and dynamic fracture [16, 17], multi-crack propagation and crack branching [18, 19], and fatigue failure[20], etc. With respect to simulating
impact problems, Huang et al. [21] investigated the damage initiation, evolution and crack propagation in a concrete model subjected to out-of-plane impact loading. Xu et al. [22] employed the peridynamic model to simulate the delamination and matrix damage process in composite laminates under low-velocity impact. For hybrid composite, Askari et al. [23] investigated high- and low-energy hail impacts against a hybrid composite to explore the effect of introducing ply-level hybridization. Up to date, the peridyanmic theory and approach has not yet been employed to analyze the dynamic fracture of laminated glass composites. As an important part of laminated glass, the PVB interlayer is a type of strongly rate-, time- and temperature-dependent polymer material. In common studies, PVB is usually described by using elastic, hyper-elastic or visco-elastic models. With peridynamic modeling, Azizi’s work [24] shown that the bond-based peridynamic model was able to predict linear viscoelastic creep behavior including recovery by using the Burger’s viscoelastic equation. Under the framework of the ordinary state-based peridynamic theory, Mitchell [25] has successfully reformulated a viscoelastic model to investigate the relaxation response in solids. On the above basis, Delorme [26] proposed a generalization of the peridynamic viscoelastic model by taking into account the viscous effects both for spherical and deviatoric behaviors. Dorduncu [27] presented a peridynamic model with “truss” element for viscoelastic structures under mechanical and thermal loads, where the truss element enabled the peridynamic interactions among the finite element nodes. And in the study of Madenci [28], the peridynamic constitutive relations for viscoelastic deformation in
terms of Prony series under mechanical and thermal loads were presented. In addition, Weckner [29] applied the Fourier- and Laplace- transforms to local and nonlocal viscoelasticity, and derived integral-representation formulas using a Green’s function approach. Based on the previously available work in literature, the present work aims at developing a non-local model and numerical approach to describe the dynamic fracture of laminated glass under low-velocity impact such that inherits the advantages of peridynamics in modeling discontinuous problems. In particular, an ordinary state-based peridynamic visco-elastic model is developed to characterize the relaxation response of the PVB interlayer, while the ordinary state-based peridynamic elastic-brittle model is employed to describe the mechanical behavior of glass sheets. Furthermore, the adhesion between glass and PVB is of vital importance to the mechanical performance of laminated glass, and which is described by using a penalty-based method in the present work. The remainder of the paper is organized as follows: Section 2 briefly introduces the fundamentals of peridynamic theory, the viscoelastic constitutive modeling, failure and damage criterion, and numerical implementation. The modeling of numerical examples is presented in Section 3. After validating the proposed model and approach through examples, the effect of thickness of PVB interlayer and glass layers on the failure of laminated glass is discussed in Section 4. Finally, a few concluding remarks are summarized in Section 5.
2. Methodology 2.1. Fundamentals of Peridynamic theory In a peridynamic solid as shown in Fig. 1, the integral-typed equation of motion for any material point x in the reference configuration at time t t 0 is given as
x, t f ξ, dVx b x, t u
(1)
Hx
where and b denote the mass density and external applied body force density respectively, u x,t is the displacement field. f denotes the interacting force vector, and dVx denotes an infinitesimal volume of material point x . In the peridynamic theory, the motion of the body is analyzed by considering the interaction of a peridynamic material point x with any material point x , within a horizon H x defined as H x 0 x x . The connection between material points is referred to as the bond, and is denoted by ξ , which is defined as ξ x x . is the relative displacement of two interacting points.
u(x,t) x
y(x,t) ξ﹢η
ξ x'
u(x',t)
y(x',t)
Fig. 1. Kinematics of material points in the reference and current configuration.
In the state-based peridynamic theory, as a generalization of the bond-based peridynamic theory, the deformed bond is given by the deformation vector state Y , which is expressed as
Y x, t ξ y x ξ, t y x, t
(3)
The Equation of motion in the state-based peridynamic theory can be written as
x, t u
Hx
T x, t
x x T x, t x x dVx b x, t
(4)
where T x,t is the force vector state that represents the relationship between material points at time t . Generally, there exists three types of peridynamic models, namely bond-based, ordinary state-based and non-ordinary state-based peridynamic models, according to the description of interactive forces between material points, as shown in Fig. 2. f' x' f
ξ
x
Bond based peridynamics '
T' T' T
x'
T
x'
ξ
ξ x
x
Ordinary state-based peridynamics
Non-ordinary state-based peridynamics
Fig. 2. Schematics of bond based, ordinary and non-ordinary state-based peridyamics.
2.2 Reformulation of ordinary state-based viscoelastic model in pedidynamics In the state-based peridynamic theory, the non-local constitutive model provides the relationship between the force vector state and the deformation vector state of a material point. The general form of a constitutive model is written as
ˆ Y, T=T
(5)
where Tˆ is bounded and Riemann-integrable on H x , and where denotes all variables other than the current deformation vector state that T depend on for some particular material.
Let M denotes the unit state valued function, which is the deformed direction vector state. It can be defined as
M Y Dir Y Thus M Y x x
Y Y
(6)
is a unit vector that points from the deformed position of x
toward the deformed position of x . For the ordinary modeling, as shown in Fig. 2, there exists a scalar state t such that
T t M
(7)
where t is termed as the scalar force state. If material is elastic and there exists a differentiable scalar valued function such that
ˆ Y T=T
(8)
The extension scalar state e , the weighted volume m , and a scalar valued function termed as the dilatation ˆ are defined respectively as follows
e y x,
y Y,
x X
m x x
ˆ e
(9) (10)
3 x e m
(11)
where is influence function. Using the weighted volume m and a scalar valued function ˆ , the isotropic part and the deviatoric part of e are defined respectively by
ei
x 3
,
ed e ei ,
ˆ e
(12)
d To characterize the viscoelastic behavior of solids, the deviatoric extension state e
de db ( i ) d de db (i ) can be divided into elastic e and back extension e , e e e (shown as
Fig. 3). Using the above decomposition for the deviatoric extension state, the elastic energy is expressed as
W ,e ,e d
db i
k 2 d e e d i e d e dbi e d e dbi 2 2 2
(13)
The scalar volumetric and deviatoric force states t i and t d are derived by
t tˆ e ▽W e ,
ti
W W = , e i e i
td
W e d
(14)
And the scalar force state t can be expressed as 3 p (15) x i e d i e dbi m where p k is the peridynamic pressure and k is the bulk modulus. The elastic t ti td =
parameter and relaxation moduli i are defined as 15 (16) i m where is a positive constant, and is shear modulus of material. In order to
avoid unbounded creep of this model, it is assumed that 0 , i 0 . In addition, for simplification we define i = , where 0 1.
Fig. 3. Standard linear solid in viscoelasticity model.
The evolution equation for the scalar deviatoric back extension plays a significant
role in the viscoelasticity model, and which can be expressed based upon the following assumption t
where t
d i
i e
d i
db i
i e d e
db i
(17)
is the scalar deviatoric force state. The evolution equation for the
standard linear solid is given by db i e
b where i
1
b i
e
d
e
db i
(18)
i is a time-concerned constant. i
The deviatoric force state t d t d e d i e d i e
db i
can be expressed as
e d i e d e
db i
e d i e
db i
(19)
2.3 Damage and failure criterion For the failure criterion of bond in peridynamics, Silling and Askari [15] defined the stretch of bond as s
Y x x X x x e x x X x x X x x
(20)
When the stretch of bond beyond a predefined limit, the bond is broken. After bond failure, there is no force sustainable between two material points. Once the bond fails, there is no healing for the bond forever, making the model history-dependent. Defining x, ,t is a history-dependent scalar-valued function mapping the breakage of the bond and is expressed as
1 s t , s0 , 0 t t 0 else
x, , t
where s0 denotes the critical stretch or failure of the bond, s0
(21) 5G0 , where G0 is 9kδ
the energy release rate [15]. The bond is active when 1 , and it is broken when
0 . Fig. 4 shows the degradation of bond force. During the dynamic failure process of glass, due to the high pressure, high strain rate, and large deformation loading condition, compression damage takes an essential role. In this study, to take the different behavior and strength of glass when subjected to impact loading, we take tensile and compressive failure into account, as proposed in [30]. The critical bond stretch is divided into tensile and compressive critical stretch as s0fc c / E , ft s0 t / E ,
s0
(22)
s0
where c , t is the uniaxial compressive strength and tensile strength, respectively.
s 0 represents the compressive state and s 0 represents the tensile state.
f μ=0 s0fc
μ=1
μ=0 s0ft
s
Fig. 4. Relationship between bond stretch and pairwise force function.
Damage in the materials is expressed locally by the ratio of the number of broken bonds to the number of the total bonds in the field as
D x, t 1
Hx
x, , t dVx
Hx
dVx
(23)
where D x, t is a function of position and time. Note that 0 D 1 , with 0
representing original with no damage, and 1 representing complete disconnection of a material point from all other points connected in its horizon [15]. 2.4 Numerical implementation A standard, single-point-based spatial discretization of Eq. (3) is given by in T xin , t n x np xin T x np , t n xin x np V p b xin , t n u p
where n denotes the number of time steps, V p = x
3
(24)
is the involved volume of
in are the acceleration of point xi at time t n . material point x p . u
It is necessary to store the velocities and displacements for each material point at the current time-step and update them at the next time-step, employing explicit forward and backward difference techniques, which are given as
u n 1 un 1t u n
(25)
u n 1 u n 1t u n
(26)
where t is time step size. To obtain a stable numerical result, the time step t should satisfy [10], t
δ c
where δ denotes the size of horizon, c
(27)
2 /
represents the dilatational
wave speed, and are the Lame’s elastic constants of the material. The adhesion between glass plies and the PVB interlayer is of significant importance to the mechanical performance of laminated glass. If laminated glass exhibits a higher adhesion level, glass plies are adhered closer to the PVB film. On the contrary, the glass shards are more likely to splash from the PVB film during impact if the adhesion level is excessively low. In this study, the adhesion between the glass plies and the PVB film is described by taking use of a penalty-based method,
which resists the relative displacements between glass plies and the PVB interlayer by adding penalty springs (as shown in Fig. 5). The penalty force f g is calculated as follows [11] f g k u g u p
(28)
where k is a penalty factor, u g and u p are the displacements of the glass node and its counterpart node of the PVB respectively.
crack
penalty springs
Fig. 5. Schematic diagram of the penalty-based method for modeling of adhesion.
3. Case studies and results In this section, a numerical example was firstly conducted on a viscoelastic plate under bending loading to verify the deformation prediction of the ordinary state-based viscoelastic model. Then, two numerical cases were established to validate the capability of the proposed model and approach for simulating impact failure of laminated glass composites, and the influence of thickness of glass and PVB interlayer on the failure of the laminated glass was investigated further. For validation, a drop-weight experiment of a laminated glass plate was presented to investigate the dynamic fracture behavior, and the results were compared with experimental results in literature [13]. Furthermore, the impact failure behavior of a laminated glass beam was analyzed and the numerical results were compared with available experimental results and results from finite element analysis in literature [31].
3.1 Viscoelastic calculation: response of rectangular plate subjected to bending Considering a 3D rectangular plate with length L 0.2 m , width W 0.1 m , and thickness h 0.003 m (see Fig. 6). The viscoelastic constitutive model described in Section 2.2 was employed, with material property constants as described in [25]. The plate is fixed along the left end, and is subjected to vertically uniform stress along the right end (external loads acting on the rightest three layers of particles). The loading history in the creep-recovery test is divided into three stages, i.e. 0-t0, t0-t1, and after t1 ( t0 0.01 s , t1 0.02 s , t2 0.03 s ), as shown in Fig. 7. The load linearly increases to 0 =20 MPa , then keeps unchanged from t0 to t1, and decreases abruptly to zero at t1. In peridynamic simulation, the grid spacing size is 0.001 m and the horizon size is three times of the grid spacing size, and the time step size is 9.0×10-7 s after convergence and stability test.
σ
x 0.2
Fig. 6. The geometric model of the Plate.
0.1
y
σ σ0
t0
t1
t2
t
Fig. 7. Loading condition in the creep-recovery test.
Fig. 8 and Fig. 9 show the vertical displacement results from the ordinary state-based viscoelastic model and the numerical results by using finite element method. Firstly, the vertical displacement increases linearly with the increase of load. When external load remains constant, the vertical displacement has a slow and nonlinear increase, as referred to the creep process. Lastly, the displacement gradually recovers to zero after the unloading time. The peridynamic predictions successfully capture the elastic response as well as the creep and recovery responses, and quantitatively agree well with FEM results.
2.4x10-6
peridynamic results FEM results
-6
Vertical displacement (m)
2.1x10
1.8x10-6 1.5x10-6 1.2x10-6 9.0x10-7 6.0x10-7 3.0x10-7 0.0 0.000
0.005
0.010
0.020
0.015
0.025
0.030
Time (s) Fig. 8. Comparison of the vertically downward displacement of the central point at the right end of the plate versus time by using Peridynamic analysis and FEM.
(a) FEM results
(b) Peridynamic results
Fig. 9. Comparison of vertical displacement contours at t=0.002s by using Peridynamic analysis (right) and FEM (left).
3.2 Impact damage and crack propagation in a laminated glass plate The geometric model for the drop-weight experiment is shown in Fig. 10, which consists of an impactor, a laminated glass plate and supporter. The configuration of the specimen is 100 mm×40 mm, and the thicknesses of top-layer glass, PVB, and bottom-layer glass are 2 mm, 0.76 mm, and 3 mm, respectively. The initial velocity of
the impactor is 6.78 m/s, and the mass is 35.5 g. The impactor is simplified as a sphere in numerical simulation. For peridynamic modeling, the grid spacing size is 0.3 mm, and the horizon size is three times of the grid spacing size. This model includes 1,195,259 material points totally. For dynamic loading, the time step is set to be 8.0×10-8 s considering both calculation efficiency, convergence and stability. The critical bond stretch of the glass is 0.0001, and the PVB interlayer keeps visco-elastic deformation during the impact.
(2+0.76+3) mm
v=6.78 m/s m=35.5 g
40 mm
10 mm
200 mm
Fig. 10. The geometric model of the laminated glass plate.
The impactor in Fig. 10 is considered as a rigid body, and the glass plies are brittle linear elastic. The PVB interlayer is modeled by using the reformulated ordinary b state-based peridynamic viscoelastic constitutive model, where material constants i
and
are 0.2 ms and 0.5 respectively [25]. Additionally, the penalty factor
employed for adhesion modeling is k EPVB / (1 2 PVB ) [32] in this study, where
EPVB and PVB are Young’s modulus and Poisson’s ratio of the PVB. Other material parameters for this specimen are listed in Table 1.
Table 1 Material properties for laminated glass plate [31]. Material
Glass
PVB
Impactor
Support
Density (kg/m3)
2500
1000
—
1200
Young's modulus (GPa) 74.0
0.0158
210.0
0.042
Poisson’s ratio
0.49
0.27
0.49
0.20
As shown in Fig. 11, compressive waves are initially generated by the impactor in the contact area of the top-layer glass and subsequently spread toward all the directions of the top glass layer. The propagation of compressive waves is restricted due to the viscoelasticity property of the PVB interlayer. Hence, the compressive waves are reflected to be tensile waves when they shift vertically down and arrive at the boundary between the top-layer glass and PVB interlayer. This leads to the fracture of the glass layers. Cracks propagate gradually when the accumulated reflected tensile waves reach the fracture criterion and form radial cracks in the top-layer glass. Furthermore, the compressive waves arrive at the bottom-layer glass after passing through the PVB interlayer. Radial cracks will form and propagate in a manner similar to the cracks in the top-layer glass. In addition, the compressive waves are reflected to be tensile waves when they spread horizontally to the free boundary and eventually form the circular cracks in glass layers.
compressive waves
v
top-layer glass PVB
tensile waves
bottomlayer glass
Fig. 11. The schematic diagram of the impact waves spreading in the laminated glass plate.
The process of the dynamic fracture for the laminated glass plate in peridynamic simulation is shown in Fig. 12. As shown in Fig. 12 (a), when the compressive waves arrive at the bottom-layer glass and are reflect as tensile waves, the radial cracks of glass begin to propagate in different directions due to the tensile waves at t = 24 μs. Significant formation of circular cracks is observed when the compressive waves reach the free boundary (Fig. 12 (d)) at t =112 μs. After that, more and more radial cracks and circular cracks initiate and propagate.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Fig. 12. Dynamic fracture processes of a laminated glass plate subjected to impact loading (a): 24 μs (b): 56 μs (c): 80 μs (d): 112 μs (e): 144 μs (f): 168 μs (g): 200 μs (h): 232 μs.
Fig. 13 shows the corresponding experimental observations and finite element analysis results on fracture pattern of the laminated glass plate in literature [13]. Totally 22 radial cracks exist in Fig. 13 (b), and most of them propagate to the specimen boundary. From the peridynamic simulation results shown in Fig. 12 (h), we can see that all 22 radial cracks can be found, and most of them propagate close to the boundary. By comparing the numbers and final morphology of radial and circular cracks, it indicates that the fracture patterns of the laminated glass plate in peridynamic simulation keep in good agreement with the experimental results [13]. As shown in Fig. 13 (a), in finite element analysis, only some primary radial cracks and circular cracks of the laminated glass plate occur [13]. In addition, the impact force from the simulation results of peridynamic and finite element analysis are
compared with that measured in the experiment, as shown in Fig. 14. Comparing the numerical results by using two different methods, it indicates that the proposed peridynamic model and approach shows its capacity and advantages in handling impact failure of laminated glass both qualitatively and quantitatively.
(a) Finite element analysis results [13]
(b) Experimental results [13] Fig. 13. Fracture of the laminated glass plate in experiments and by using finite element analysis [13]. 2500
peridynamic simulation experiment finite element simulation
Force (N)
2000
1500
1000
500
0 0.0
0.1
0.2
0.3
0.4
0.5
Time (ms)
Fig. 14. The impact force-time curve of the impactor for different results.
3.3 Effects of the thicknesses of PVB interlayer/glass layer on impact failure of laminated glass To further demonstrate the capability of the proposed peridynamic model and approach in simulating impact failure of laminated glass, in the second example, we simulate the fracture process of a laminated glass beam under impact loading and compare the results with experimental observations and finite element analysis with cohesive zone model in literature [31]. And then, an in-depth investigation was conducted to analyze the effect of thickness of the PVB interlayer and glass layer on glass-ply cracking, since available researches [13, 31] have shown that the PVB thickness may take great effects on damage and crack propagation from top-layer to bottom-layer glass. The model is shown in Fig. 15, which includes an impactor, a laminated glass beam and supporter. In peridynamic simulation, the model is discretized to 23608 material points with a grid spacing of 1.5 mm. The time step is 4.0×10-8 s, and the horizon size is 4.5 mm. Velocity of the impactor is 3.13 m/s. Other material parameters keep the same as listed in Table 1.
v=3.13 m/s m=0.1 kg
(10+4+10) mm
10 mm
10 mm
10 mm
200 mm
Fig. 15. The geometric model of the laminated glass beam.
The dynamic fracture process of the laminated glass beam by using the proposed peridynamic approach and stress contours and crack propagation results from finite element analysis on the same model in reference [31] are shown in Fig. 16. The results show that due to the accumulation of the compressive waves in top-layer glass, the crack occurs first from the bottom of the top-layer glass, because the accumulated tensile stress there reaches the tensile strength of glass firstly. Then, the crack propagates upwardly until reaching the boundary of the top-layer glass. In peridynmaic simulation, the first crack occurs at 52 μs, and penetrates through the top-layer glass at 110 μs. While in finite element analysis, the crack initiates at about 59 μs propagates to the top surface at 133μs. In peridynamic simulation, same as experimental observations, the compressive waves will pass through the PVB film, and arrive at the bottom-layer glass. The crack in the bottom-layer glass occurs at about 175 μs, and eventually penetrates through the bottom-layer glass at 280 μs. In the present low-velocity impact loading cases, the principal crack pattern in the laminated glass beam is vertical crack in both top and
bottom glasses, and delamination between glass and PVB does not occur. Fig. 17 shows the corresponding experimental observations on crack initiation and propagation in the top-layer glass of the laminated glass beam with a time interval of 20 μs. It evidently declares that the crack initiates from the bottom side of the top-layer glass and propagates upwards to the top surface of the beam. By comparison, we can find again that the crack morphology, initial locations and propagation obtained by using the proposed peridynamic approach keep good agreement with the results of corresponding experiment and FEM results.
(a1)
(a2)
(b1)
(b2)
(c1)
(c2)
Fig. 16. Dynamic fracture process of a laminated glass beam subjected to impact loading in peridynmaic simulation (left), and finite element simulation (right). (a1). 52 μs, (b1). 60 μs, (c1). 205 μs; (a2). 55 μs, (b2). 65 μs, (c2). 234 μs.
Fig. 17. Experimental observations: crack initiation and propagation in the top-layer glass (with time interval of 20 μs) [31].
After qualitative and quantitative verification, five cases with different PVB interlayer thicknesses: 2 mm, 3 mm, 4mm, 5 mm, and 6 mm were performed to investigate the effect of the thickness of the PVB interlayer on impact failure of the laminated glass. The final fracture modes of the specimen are shown in Fig. 18. The final crack patterns of five cases are almost identical irrespective of the thickness of the PVB interlayer, thereby indicating that the variation of the thicknesses of the PVB interlayer takes little effect on the fracture mode of the glass under low-velocity impact.
(a)
(b)
(c)
(d)
(e) Fig. 18. Damage contours in the laminated glass beam with different PVB interlayer thickness (a): 2 mm, (b): 3 mm, (c): 4 mm, (d): 5 mm, (e): 6 mm.
To further investigate the effect of the thickness of the PVB interlayer, the plots of time variation versus crack initiation and growth in the bottom-layer glass with different PVB interlayer thicknesses are shown in Fig. 19. It indicates that the PVB thickness significantly influences the crack initiation time and has a slight influence on crack propagation velocity in the bottom-layer glass. The crack in the bottom-layer glass with a thinner PVB layer will initiate much earlier, demonstrating that the speed with which compressive waves pass through the PVB interlayer decreases when the PVB interlayer thickness increases, which results in the differences of crack initiation time in bottom-layer glass. However, in bottom-layer glass, the influence on the crack propagation velocity is much slighter. Fig. 19 shows that before the crack length reaches 8 mm, the crack propagation velocity with different PVB interlayer thicknesses is almost identical. After that, the crack propagates a little faster with
thinner PVB interlayer. It is worth to note that a special phenomenon can be found in Fig. 19. When the crack length reaches about 7 mm, the speed of crack propagation will decrease, no matter the thickness of the PVB interlayer changes.
360
2 mm 3 mm 4 mm 5 mm 6 mm
330 300
Time (s)
270 240 210 180 150 120 0
2
4
6
8
10
12
Crack length (mm)
Fig. 19. Plots of time versus crack initiation and growth with different PVB interlayer thickness.
Finally, five cases are conducted to analyze the effect of the thickness of glass layer. For the sake of simplicity, the thickness of the bottom-layer glass keeps 10 mm and the thickness of the top-layer glass changes from 6 mm, 8 mm, 10 mm, 12 mm to 14 mm. The thickness of the PVB interlayer keeps 4 mm, and the impactor velocity is 3.13 m/s. The final fracture of the laminated glass beam and the time variation versus crack length are plotted in Fig. 20 and Fig. 21, respectively.
(a)
(b)
(c)
(d)
(e) Fig. 20. Damage contours in the laminated glass beam with different top-layer glass thickness, (a): 6 mm, (b): 8 mm, (c): 10 mm, (d): 12 mm, (e): 14 mm.
360
6 mm 8 mm 10 mm 12 mm 14 mm
330 300
Time (s)
270 240 210 180 150 120 90 0
2
4
6
8
10
12
Crack length(mm)
Fig. 21. Plots of time variation versus crack length with different top-layer glass thickness.
As shown in Fig. 20, the final fracture modes for the former four cases are nearly identical to those previously. In the last case, the fracture of the bottom-layer glass is different from that of others. Because of the ultra-high value of the thickness of the top-layer glass (14mm), the crack in bottom-layer glass stops propagating at a length of approximately 3.3 mm. With respect to other four cases, the thickness of the top-layer glass affects the crack initiation time obviously in Fig. 21. The crack in the bottom-layer glass with a thinner top-layer glass initiates much earlier, which is similar to the influence of the thickness of the PVB interlayer. However, the crack propagation speed keeps almost unchanged when the thickness of the top-layer glass increases. Similarly, just like the observations mentioned above, the crack propagation speed will significantly decrease when the crack length reaches about 7 mm in call cases.
4. Summary
In this study, the dynamic fracture behavior of laminated glass composites under low-velocity impact is modeled under the framework of the state-based peridynamic theory for the first time. The relaxation response of the PVB interlayer is described by using a peridynamic reformulation of viscoelasticity. A penalty-based method is employed to characterize the adhesion between the glass plies and PVB interlayer, which is of key importance to the mechanical performance of laminated glass. The proposed model and numerical approach were validated through simulating the failure of a laminated glass plate under drop-weight impact loading, and comparing the
numerical results with available experimental observations. After model validation, the effects of the PVB interlayer thickness and glass layer thickness on glass-ply cracking were investigated through the simulation of a laminated glass beam. We compare the crack morphology, initial locations and propagation sequence with experiment results and finite element model results. In various cases, the glass layer keeps the same failure mode. Nevertheless, both the PVB thickness and glass thickness take effects on the crack initiation time in the bottom-layer glass. The crack in bottom-layer glass with thinner PVB interlayer or thinner top-layer glass will initiate significantly earlier. Furthermore, the speed of the crack propagation keeps nearly unchanged when the PVB thicknesses and glass thicknesses increase. In all simulation results for laminated glass with different PVB interlayer thicknesses or glass layer thicknesses, it shows that the crack propagation speed will significantly decrease when the crack length reaches some critical value (e.g. 7 mm under the working cases in this study). Further studies will focus on the mechanism of the crack growth.
Acknowledgements: The authors acknowledge the support from the Fundamental Research Funds for the Central Universities in China (No. 2017B13014, 2019B66014), the Research & Innovation Project of Jiangsu Province (No. SJKY19_0419), the National Natural Science Foundation of China (No. 51679077, 11932006), and the National Key R&D Program of China (No. 2018YFC0406703).
References [1]
Muralidhar S, Jagota A, Bennison SJ, Saigal S. Mechanical behaviour in tension of cracked glass bridged by an elastomeric ligament. Acta Mater 2000, 48(18-19):4577-88.
[2]
Xu J, Li YB, Liu BH, Zhu MY, Ge DY. Experimental study on mechanical behavior of PVB laminated glass under quasi-static and dynamic loadings. Compos B 2011;42(2):302-8.
[3]
Chen JJ, Xu J, Liu BH, Yao XF, Li YB. Quantity effect of radial cracks on the
cracking propagation behavior and the crack morphology. PLOS ONE
2014;9(7): e98196. [4]
[5]
Chen JJ, Xu J, Yao XF, Xu XQ, Liu BH, Li YB. Different driving mechanisms
of in-plane cracking on two brittle layers of laminated glass. Int
J Impact Eng
2014; 69:80-5.
Zhang XH, Hao H, Ma G. Laboratory test and numerical simulation of laminated glass window vulnerability to debris impact. Int J Impact Eng 2013; 55:49-62.
[6]
Pyttel T, Liebertz H, Cai J. Failure criterion for laminated glass under impact loading and its application in fnite element simulation, Int J Impact Eng 2011;38 (4):252-63.
[7]
Hidallana-Gamage HD, Thambiratnam DP, Perera NJ. Failure analysis of laminated glass panels subjected to blast loads. Eng Fail Anal 2014;36:1429.
[8]
Peng Y, Yang JK, Deck C, Willinger R. Finite element modeling of crash test
behavior
for
windshield
laminated
glass,
Int
J
Impact
Eng
2013;57:27-35. [9]
Xu J, Li YB. Study of damage in windshield glazing subject to impact by a pedestrian’s head, Proc Inst Mech Eng D 2009;223(D1):77-84.
[10]
Zang MY, Lei Z, Wang SF. Investigation of impact fracture behavior of automobile laminated glass by 3D discrete element method. Comput Mech
2007;41 (1):73-83. [11]
Chen SH, Zang MY, Wang D, Yoshimura S, Yamada T. Numerical analysis of impact failure of automotive laminated glass: A review. Compos Pt B-Eng 2017; 122:47-60.
[12]
Xu J, Li YB, Chen X, Yan Y, Ge DY, Zhu MY, et al. Characteristics of windshield cracking upon low-speed impact: numerical simulation based on the extended finite element method, Comput Mater Sci 2010;48(3):582-88.
[13]
Chen SH, Zang MY, Xu W. A three-dimensional computational framework for
impact fracture analysis of automotive laminated glass. Comput Methods
Appl Mech Eng 2015; 294:72-99. [14]
Silling SA. Reformulation of elasticity theory for discontinuities and longrange forces. J Mech Phys Solids 2000;48(1):175-209.
[15]
Silling SA, Askari E. A meshfree method based on the peridynamic model of solid mechanics. Comput Struct 2005;83(17-18):1526-35.
[16]
Breitenfeld MS, Geubelle PH, Weckner O, Silling SA. Non-ordinary statebased peridynamic analysis of stationary crack problems. Comput Method Appl Mech Eng 2014; 272:233-250.
[17]
Silling SA. Dynamic fracture modeling with a meshfree peridynamic code. In: Proceedings Second MIT Conference on Computational Fluid and Solid Mechanics, 17-20 June 2003.P. 641-44.
[18]
Kilic B, Madenci E. Prediction of crack paths in a quenched glass plate by using peridynamic theory. Int J Fract 2009;156(2):165-77.
[19]
Ha YD, Bobaru F. Studies of dynamic crack propagation and crack branching with Peridynamics. Int J Fract 2010;162(1-2):229-44.
[20]
Hu YL, Madenci E. Peridynamics for fatigue life and residual strength prediction of composite laminates. Compos Struct 2017; 160:169-84.
[21]
Huang D, Zhang Q, Qiao PZ. Damage and progressive failure of concrete structures using non-local peridynamic modeling. Sci China Technol Sci; 2011;54(3):591-6.
[22]
Xu J, Askari A, Weckner O, Silling SA. Peridynamic analysis of impact
damage in composite laminates. J Aerospace Eng 2008;21(3):187-94. [23]
Askari A, Nelson K, Weckner O, Xu JF, Silling SA. Hail impact characteristics of a hybrid material by advanced analysis techniques and testing. J Aerospace Eng 2011;24(2):210-7.
[24]
Azizi M, bin Mohd Ihsan A, Mohamed N. The peridynamic model of viscoelastic creep and recovery. Multidiscipline Modeling in Materials and Structures 2015; 11(4):579-97.
[25]
Mitchell JA. A non-local, ordinary-state-based viscoelasticity model for peridynamics. Sandia National Lab Report, 2011.
[26]
Delorme R, Tabiai I, Lebel L, Lévesque M, Generalization of the ordinary state-based peridynamic model for isotropic linear viscoelasticity. Mech Time-Depend Mater 2017; 21:549-75.
[27]
Dorduncu M, Barut A, Madenci E. Peridynamic truss element for viscoelastic
deformation,
57th
AIAA/ASCE/AHS/ASC
Structures,
Structural Dynamics, and Materials Conference, AIAA SciTech. [28]
Madenci
E,
Oterkus
S.
Ordinary
state-based
peridynamics
for
thermoviscoelastic deformation. Eng Fract Mech 2017; 175:31-45. [29]
Weckner O, Mohamed N. Viscoelastic material models in peridynamics, Appl Math Comput 2013; 219:6039-43
[30]
Kilic B, Madenci E. Structural stability and failure analysis using peridynamic theory. Int J Non-Linear Mech 2009; 44(8):845-54.
[31]
Chen SH, Zang MY, Wang D, Zheng ZM, Zhao CL. Finite element modelling
of impact damage in polyvinyl butyral laminated glass. Compos
Struct 2016; 138:1-11. [32]
Zang MY, Lei Z. An approach to combining 3D discrete and finite element methods based on penalty function method. Comput Mech 2010;46(4):60919.
S0263-8223(19)33108-3 https://doi.org/10.1016/j.compstruct.2019.111722 COST 111722
To appear in:
Composite Structures
Received Date: Revised Date: Accepted Date:
18 August 2019 2 October 2019 21 November 2019
Please cite this article as: Wu, L., Wang, L., Huang, D., Xu, Y., An ordinary state-based peridynamic modeling for dynamic fracture of laminated glass under low-velocity impact, Composite Structures (2019), doi: https://doi.org/ 10.1016/j.compstruct.2019.111722
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.
© 2019 Published by Elsevier Ltd.
An ordinary state-based peridynamic modeling for dynamic fracture of laminated glass under low-velocity impact Liwei Wu1, Lei Wang1, Dan Huang1 *, Yepeng Xu1 College of Mechanics and Materials, Hohai University, Nanjing, 211100, China.
Abstract: Laminated glass, as a typical kind of sandwich-structure composite, is widely used in various fields such as safety, vehicle and transportation engineering. Laminated glass mainly suffers from dynamic/impact loadings, and glass-ply cracking is the main failure pattern of laminated glass. This paper presents a non-local ordinary state-based peridynamic modeling and numerical approach for simulating the dynamic failure process of sandwiched laminated glass under low-velocity impact loading. The mechanical behavior of the PVB layer is simulated by reformulating classical visco-elastic model under the framework of ordinary state-based peridynamic theory, and the adhesion between the glass and PVB interlayer is described by using a penalty-based method. Fracture patterns of a laminated glass plate under drop-weight loading were investigated, and the numerical results compared well with experimental observations. Furthermore, a series of numerical simulations were conducted to in-depth analyze the effect of the thickness of the PVB interlayer and the glass layers on the fracture mode, initial locations and crack propagation speed of the laminated glass when subjected to impact. Keywords: Laminated glass; impact failure; peridynamics; non-local; visco-elastic
*Corresponding
Author: Dan Huang. Email: [email protected]
1. Introduction As a type of widely used sandwich composite structure, laminated glass not only owns the capacity of shielding for wind, rain, heat and sound, but also plays a significant role as a safety component in infrastructure as well as military fields. Typical laminated glass usually consists of two glass plies adhered by a soft interlayer, preferably polyvinyl butyral (PVB). Laminated glass suffers from low-velocity impact in various cases, such as impact of a small and low-velocity stone thrown by the wheel in automotive applications, and impact of wind-borne debris in architectural applications. Effectively declaring the failure mechanism of laminated glass when subjected to low-velocity impact is of important value for improving the resistance as well as safety of laminated glass. This motivates the present work to study the characteristics of dynamic fracture of the laminated glass under low-velocity impact. In the past decades, numerous experimental and numerical studies have been conducted on the mechanical behavior of laminated glass, particularly on the failure mechanism and fracture modeling. For example, Muralidhar et al. [1] and Xu et al. [2] experimentally investigated the debonding of glass-polymer interface under quasi-static and dynamic loading. Chen et al. [3, 4] experimentally studied the cracking behavior of the glass layer in a laminated glass panel. As to the modeling for failure of laminated glass, a series of numerical methods have been proposed in which the most commonly used is the element deletion or erosion contact technologies [5-7] based on finite element method, and incorporated continuum damage mechanics (CDM) method [8, 9]. Zang et al. [10] developed a three-dimensional discrete element
model to simulate the dynamic fracture of laminated glass, and Zhang and Lei [11] proposed a combined discrete element/finite element method (DEM/FEMs) based on a penalty function method to improve the calculation efficiency. Furthermore, Xu et al. [12] simulated performed a parametric study on the fracture of laminated glass plate subjected to low-velocity impact with the extended finite element method (XFEM), and Chen et al. [13] analyzed the crack propagation in the glass layer of a laminated glass plate under drop-weight impact by using the extrinsic cohesive zone model (CZM). In recent years, another alternative titled peridynamics (PD) [14] exhibits significant potential in handling fracture-concerned problems, in particular, simulating impact failure. The recently proposed non-local peridynamic theory [14] is a reformulation of the classical continuum mechanics particularly for handling problems with discontinuities and long-range forces. As opposed to conventional continuum mechanics theory, an integral-typed equation of motion without spatial derivatives is employed to describe the mechanical behavior of solids. The primary advantage of the peridynamic theory lies in its validity in the presence of discontinuities and its ability to model discontinuities naturally and spontaneously without any extra technical remedies employed in conventional continuous methods [15]. Owing to its unique advantages in naturally incorporating discontinuities within a single continuum framework, the peridynamic theory has been successfully applied to study various fracture problems, including quasi-static and dynamic fracture [16, 17], multi-crack propagation and crack branching [18, 19], and fatigue failure[20], etc. With respect to simulating
impact problems, Huang et al. [21] investigated the damage initiation, evolution and crack propagation in a concrete model subjected to out-of-plane impact loading. Xu et al. [22] employed the peridynamic model to simulate the delamination and matrix damage process in composite laminates under low-velocity impact. For hybrid composite, Askari et al. [23] investigated high- and low-energy hail impacts against a hybrid composite to explore the effect of introducing ply-level hybridization. Up to date, the peridyanmic theory and approach has not yet been employed to analyze the dynamic fracture of laminated glass composites. As an important part of laminated glass, the PVB interlayer is a type of strongly rate-, time- and temperature-dependent polymer material. In common studies, PVB is usually described by using elastic, hyper-elastic or visco-elastic models. With peridynamic modeling, Azizi’s work [24] shown that the bond-based peridynamic model was able to predict linear viscoelastic creep behavior including recovery by using the Burger’s viscoelastic equation. Under the framework of the ordinary state-based peridynamic theory, Mitchell [25] has successfully reformulated a viscoelastic model to investigate the relaxation response in solids. On the above basis, Delorme [26] proposed a generalization of the peridynamic viscoelastic model by taking into account the viscous effects both for spherical and deviatoric behaviors. Dorduncu [27] presented a peridynamic model with “truss” element for viscoelastic structures under mechanical and thermal loads, where the truss element enabled the peridynamic interactions among the finite element nodes. And in the study of Madenci [28], the peridynamic constitutive relations for viscoelastic deformation in
terms of Prony series under mechanical and thermal loads were presented. In addition, Weckner [29] applied the Fourier- and Laplace- transforms to local and nonlocal viscoelasticity, and derived integral-representation formulas using a Green’s function approach. Based on the previously available work in literature, the present work aims at developing a non-local model and numerical approach to describe the dynamic fracture of laminated glass under low-velocity impact such that inherits the advantages of peridynamics in modeling discontinuous problems. In particular, an ordinary state-based peridynamic visco-elastic model is developed to characterize the relaxation response of the PVB interlayer, while the ordinary state-based peridynamic elastic-brittle model is employed to describe the mechanical behavior of glass sheets. Furthermore, the adhesion between glass and PVB is of vital importance to the mechanical performance of laminated glass, and which is described by using a penalty-based method in the present work. The remainder of the paper is organized as follows: Section 2 briefly introduces the fundamentals of peridynamic theory, the viscoelastic constitutive modeling, failure and damage criterion, and numerical implementation. The modeling of numerical examples is presented in Section 3. After validating the proposed model and approach through examples, the effect of thickness of PVB interlayer and glass layers on the failure of laminated glass is discussed in Section 4. Finally, a few concluding remarks are summarized in Section 5.
2. Methodology 2.1. Fundamentals of Peridynamic theory In a peridynamic solid as shown in Fig. 1, the integral-typed equation of motion for any material point x in the reference configuration at time t t 0 is given as
x, t f ξ, dVx b x, t u
(1)
Hx
where and b denote the mass density and external applied body force density respectively, u x,t is the displacement field. f denotes the interacting force vector, and dVx denotes an infinitesimal volume of material point x . In the peridynamic theory, the motion of the body is analyzed by considering the interaction of a peridynamic material point x with any material point x , within a horizon H x defined as H x 0 x x . The connection between material points is referred to as the bond, and is denoted by ξ , which is defined as ξ x x . is the relative displacement of two interacting points.
u(x,t) x
y(x,t) ξ﹢η
ξ x'
u(x',t)
y(x',t)
Fig. 1. Kinematics of material points in the reference and current configuration.
In the state-based peridynamic theory, as a generalization of the bond-based peridynamic theory, the deformed bond is given by the deformation vector state Y , which is expressed as
Y x, t ξ y x ξ, t y x, t
(3)
The Equation of motion in the state-based peridynamic theory can be written as
x, t u
Hx
T x, t
x x T x, t x x dVx b x, t
(4)
where T x,t is the force vector state that represents the relationship between material points at time t . Generally, there exists three types of peridynamic models, namely bond-based, ordinary state-based and non-ordinary state-based peridynamic models, according to the description of interactive forces between material points, as shown in Fig. 2. f' x' f
ξ
x
Bond based peridynamics '
T' T' T
x'
T
x'
ξ
ξ x
x
Ordinary state-based peridynamics
Non-ordinary state-based peridynamics
Fig. 2. Schematics of bond based, ordinary and non-ordinary state-based peridyamics.
2.2 Reformulation of ordinary state-based viscoelastic model in pedidynamics In the state-based peridynamic theory, the non-local constitutive model provides the relationship between the force vector state and the deformation vector state of a material point. The general form of a constitutive model is written as
ˆ Y, T=T
(5)
where Tˆ is bounded and Riemann-integrable on H x , and where denotes all variables other than the current deformation vector state that T depend on for some particular material.
Let M denotes the unit state valued function, which is the deformed direction vector state. It can be defined as
M Y Dir Y Thus M Y x x
Y Y
(6)
is a unit vector that points from the deformed position of x
toward the deformed position of x . For the ordinary modeling, as shown in Fig. 2, there exists a scalar state t such that
T t M
(7)
where t is termed as the scalar force state. If material is elastic and there exists a differentiable scalar valued function such that
ˆ Y T=T
(8)
The extension scalar state e , the weighted volume m , and a scalar valued function termed as the dilatation ˆ are defined respectively as follows
e y x,
y Y,
x X
m x x
ˆ e
(9) (10)
3 x e m
(11)
where is influence function. Using the weighted volume m and a scalar valued function ˆ , the isotropic part and the deviatoric part of e are defined respectively by
ei
x 3
,
ed e ei ,
ˆ e
(12)
d To characterize the viscoelastic behavior of solids, the deviatoric extension state e
de db ( i ) d de db (i ) can be divided into elastic e and back extension e , e e e (shown as
Fig. 3). Using the above decomposition for the deviatoric extension state, the elastic energy is expressed as
W ,e ,e d
db i
k 2 d e e d i e d e dbi e d e dbi 2 2 2
(13)
The scalar volumetric and deviatoric force states t i and t d are derived by
t tˆ e ▽W e ,
ti
W W = , e i e i
td
W e d
(14)
And the scalar force state t can be expressed as 3 p (15) x i e d i e dbi m where p k is the peridynamic pressure and k is the bulk modulus. The elastic t ti td =
parameter and relaxation moduli i are defined as 15 (16) i m where is a positive constant, and is shear modulus of material. In order to
avoid unbounded creep of this model, it is assumed that 0 , i 0 . In addition, for simplification we define i = , where 0 1.
Fig. 3. Standard linear solid in viscoelasticity model.
The evolution equation for the scalar deviatoric back extension plays a significant
role in the viscoelasticity model, and which can be expressed based upon the following assumption t
where t
d i
i e
d i
db i
i e d e
db i
(17)
is the scalar deviatoric force state. The evolution equation for the
standard linear solid is given by db i e
b where i
1
b i
e
d
e
db i
(18)
i is a time-concerned constant. i
The deviatoric force state t d t d e d i e d i e
db i
can be expressed as
e d i e d e
db i
e d i e
db i
(19)
2.3 Damage and failure criterion For the failure criterion of bond in peridynamics, Silling and Askari [15] defined the stretch of bond as s
Y x x X x x e x x X x x X x x
(20)
When the stretch of bond beyond a predefined limit, the bond is broken. After bond failure, there is no force sustainable between two material points. Once the bond fails, there is no healing for the bond forever, making the model history-dependent. Defining x, ,t is a history-dependent scalar-valued function mapping the breakage of the bond and is expressed as
1 s t , s0 , 0 t t 0 else
x, , t
where s0 denotes the critical stretch or failure of the bond, s0
(21) 5G0 , where G0 is 9kδ
the energy release rate [15]. The bond is active when 1 , and it is broken when
0 . Fig. 4 shows the degradation of bond force. During the dynamic failure process of glass, due to the high pressure, high strain rate, and large deformation loading condition, compression damage takes an essential role. In this study, to take the different behavior and strength of glass when subjected to impact loading, we take tensile and compressive failure into account, as proposed in [30]. The critical bond stretch is divided into tensile and compressive critical stretch as s0fc c / E , ft s0 t / E ,
s0
(22)
s0
where c , t is the uniaxial compressive strength and tensile strength, respectively.
s 0 represents the compressive state and s 0 represents the tensile state.
f μ=0 s0fc
μ=1
μ=0 s0ft
s
Fig. 4. Relationship between bond stretch and pairwise force function.
Damage in the materials is expressed locally by the ratio of the number of broken bonds to the number of the total bonds in the field as
D x, t 1
Hx
x, , t dVx
Hx
dVx
(23)
where D x, t is a function of position and time. Note that 0 D 1 , with 0
representing original with no damage, and 1 representing complete disconnection of a material point from all other points connected in its horizon [15]. 2.4 Numerical implementation A standard, single-point-based spatial discretization of Eq. (3) is given by in T xin , t n x np xin T x np , t n xin x np V p b xin , t n u p
where n denotes the number of time steps, V p = x
3
(24)
is the involved volume of
in are the acceleration of point xi at time t n . material point x p . u
It is necessary to store the velocities and displacements for each material point at the current time-step and update them at the next time-step, employing explicit forward and backward difference techniques, which are given as
u n 1 un 1t u n
(25)
u n 1 u n 1t u n
(26)
where t is time step size. To obtain a stable numerical result, the time step t should satisfy [10], t
δ c
where δ denotes the size of horizon, c
(27)
2 /
represents the dilatational
wave speed, and are the Lame’s elastic constants of the material. The adhesion between glass plies and the PVB interlayer is of significant importance to the mechanical performance of laminated glass. If laminated glass exhibits a higher adhesion level, glass plies are adhered closer to the PVB film. On the contrary, the glass shards are more likely to splash from the PVB film during impact if the adhesion level is excessively low. In this study, the adhesion between the glass plies and the PVB film is described by taking use of a penalty-based method,
which resists the relative displacements between glass plies and the PVB interlayer by adding penalty springs (as shown in Fig. 5). The penalty force f g is calculated as follows [11] f g k u g u p
(28)
where k is a penalty factor, u g and u p are the displacements of the glass node and its counterpart node of the PVB respectively.
crack
penalty springs
Fig. 5. Schematic diagram of the penalty-based method for modeling of adhesion.
3. Case studies and results In this section, a numerical example was firstly conducted on a viscoelastic plate under bending loading to verify the deformation prediction of the ordinary state-based viscoelastic model. Then, two numerical cases were established to validate the capability of the proposed model and approach for simulating impact failure of laminated glass composites, and the influence of thickness of glass and PVB interlayer on the failure of the laminated glass was investigated further. For validation, a drop-weight experiment of a laminated glass plate was presented to investigate the dynamic fracture behavior, and the results were compared with experimental results in literature [13]. Furthermore, the impact failure behavior of a laminated glass beam was analyzed and the numerical results were compared with available experimental results and results from finite element analysis in literature [31].
3.1 Viscoelastic calculation: response of rectangular plate subjected to bending Considering a 3D rectangular plate with length L 0.2 m , width W 0.1 m , and thickness h 0.003 m (see Fig. 6). The viscoelastic constitutive model described in Section 2.2 was employed, with material property constants as described in [25]. The plate is fixed along the left end, and is subjected to vertically uniform stress along the right end (external loads acting on the rightest three layers of particles). The loading history in the creep-recovery test is divided into three stages, i.e. 0-t0, t0-t1, and after t1 ( t0 0.01 s , t1 0.02 s , t2 0.03 s ), as shown in Fig. 7. The load linearly increases to 0 =20 MPa , then keeps unchanged from t0 to t1, and decreases abruptly to zero at t1. In peridynamic simulation, the grid spacing size is 0.001 m and the horizon size is three times of the grid spacing size, and the time step size is 9.0×10-7 s after convergence and stability test.
σ
x 0.2
Fig. 6. The geometric model of the Plate.
0.1
y
σ σ0
t0
t1
t2
t
Fig. 7. Loading condition in the creep-recovery test.
Fig. 8 and Fig. 9 show the vertical displacement results from the ordinary state-based viscoelastic model and the numerical results by using finite element method. Firstly, the vertical displacement increases linearly with the increase of load. When external load remains constant, the vertical displacement has a slow and nonlinear increase, as referred to the creep process. Lastly, the displacement gradually recovers to zero after the unloading time. The peridynamic predictions successfully capture the elastic response as well as the creep and recovery responses, and quantitatively agree well with FEM results.
2.4x10-6
peridynamic results FEM results
-6
Vertical displacement (m)
2.1x10
1.8x10-6 1.5x10-6 1.2x10-6 9.0x10-7 6.0x10-7 3.0x10-7 0.0 0.000
0.005
0.010
0.020
0.015
0.025
0.030
Time (s) Fig. 8. Comparison of the vertically downward displacement of the central point at the right end of the plate versus time by using Peridynamic analysis and FEM.
(a) FEM results
(b) Peridynamic results
Fig. 9. Comparison of vertical displacement contours at t=0.002s by using Peridynamic analysis (right) and FEM (left).
3.2 Impact damage and crack propagation in a laminated glass plate The geometric model for the drop-weight experiment is shown in Fig. 10, which consists of an impactor, a laminated glass plate and supporter. The configuration of the specimen is 100 mm×40 mm, and the thicknesses of top-layer glass, PVB, and bottom-layer glass are 2 mm, 0.76 mm, and 3 mm, respectively. The initial velocity of
the impactor is 6.78 m/s, and the mass is 35.5 g. The impactor is simplified as a sphere in numerical simulation. For peridynamic modeling, the grid spacing size is 0.3 mm, and the horizon size is three times of the grid spacing size. This model includes 1,195,259 material points totally. For dynamic loading, the time step is set to be 8.0×10-8 s considering both calculation efficiency, convergence and stability. The critical bond stretch of the glass is 0.0001, and the PVB interlayer keeps visco-elastic deformation during the impact.
(2+0.76+3) mm
v=6.78 m/s m=35.5 g
40 mm
10 mm
200 mm
Fig. 10. The geometric model of the laminated glass plate.
The impactor in Fig. 10 is considered as a rigid body, and the glass plies are brittle linear elastic. The PVB interlayer is modeled by using the reformulated ordinary b state-based peridynamic viscoelastic constitutive model, where material constants i
and
are 0.2 ms and 0.5 respectively [25]. Additionally, the penalty factor
employed for adhesion modeling is k EPVB / (1 2 PVB ) [32] in this study, where
EPVB and PVB are Young’s modulus and Poisson’s ratio of the PVB. Other material parameters for this specimen are listed in Table 1.
Table 1 Material properties for laminated glass plate [31]. Material
Glass
PVB
Impactor
Support
Density (kg/m3)
2500
1000
—
1200
Young's modulus (GPa) 74.0
0.0158
210.0
0.042
Poisson’s ratio
0.49
0.27
0.49
0.20
As shown in Fig. 11, compressive waves are initially generated by the impactor in the contact area of the top-layer glass and subsequently spread toward all the directions of the top glass layer. The propagation of compressive waves is restricted due to the viscoelasticity property of the PVB interlayer. Hence, the compressive waves are reflected to be tensile waves when they shift vertically down and arrive at the boundary between the top-layer glass and PVB interlayer. This leads to the fracture of the glass layers. Cracks propagate gradually when the accumulated reflected tensile waves reach the fracture criterion and form radial cracks in the top-layer glass. Furthermore, the compressive waves arrive at the bottom-layer glass after passing through the PVB interlayer. Radial cracks will form and propagate in a manner similar to the cracks in the top-layer glass. In addition, the compressive waves are reflected to be tensile waves when they spread horizontally to the free boundary and eventually form the circular cracks in glass layers.
compressive waves
v
top-layer glass PVB
tensile waves
bottomlayer glass
Fig. 11. The schematic diagram of the impact waves spreading in the laminated glass plate.
The process of the dynamic fracture for the laminated glass plate in peridynamic simulation is shown in Fig. 12. As shown in Fig. 12 (a), when the compressive waves arrive at the bottom-layer glass and are reflect as tensile waves, the radial cracks of glass begin to propagate in different directions due to the tensile waves at t = 24 μs. Significant formation of circular cracks is observed when the compressive waves reach the free boundary (Fig. 12 (d)) at t =112 μs. After that, more and more radial cracks and circular cracks initiate and propagate.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Fig. 12. Dynamic fracture processes of a laminated glass plate subjected to impact loading (a): 24 μs (b): 56 μs (c): 80 μs (d): 112 μs (e): 144 μs (f): 168 μs (g): 200 μs (h): 232 μs.
Fig. 13 shows the corresponding experimental observations and finite element analysis results on fracture pattern of the laminated glass plate in literature [13]. Totally 22 radial cracks exist in Fig. 13 (b), and most of them propagate to the specimen boundary. From the peridynamic simulation results shown in Fig. 12 (h), we can see that all 22 radial cracks can be found, and most of them propagate close to the boundary. By comparing the numbers and final morphology of radial and circular cracks, it indicates that the fracture patterns of the laminated glass plate in peridynamic simulation keep in good agreement with the experimental results [13]. As shown in Fig. 13 (a), in finite element analysis, only some primary radial cracks and circular cracks of the laminated glass plate occur [13]. In addition, the impact force from the simulation results of peridynamic and finite element analysis are
compared with that measured in the experiment, as shown in Fig. 14. Comparing the numerical results by using two different methods, it indicates that the proposed peridynamic model and approach shows its capacity and advantages in handling impact failure of laminated glass both qualitatively and quantitatively.
(a) Finite element analysis results [13]
(b) Experimental results [13] Fig. 13. Fracture of the laminated glass plate in experiments and by using finite element analysis [13]. 2500
peridynamic simulation experiment finite element simulation
Force (N)
2000
1500
1000
500
0 0.0
0.1
0.2
0.3
0.4
0.5
Time (ms)
Fig. 14. The impact force-time curve of the impactor for different results.
3.3 Effects of the thicknesses of PVB interlayer/glass layer on impact failure of laminated glass To further demonstrate the capability of the proposed peridynamic model and approach in simulating impact failure of laminated glass, in the second example, we simulate the fracture process of a laminated glass beam under impact loading and compare the results with experimental observations and finite element analysis with cohesive zone model in literature [31]. And then, an in-depth investigation was conducted to analyze the effect of thickness of the PVB interlayer and glass layer on glass-ply cracking, since available researches [13, 31] have shown that the PVB thickness may take great effects on damage and crack propagation from top-layer to bottom-layer glass. The model is shown in Fig. 15, which includes an impactor, a laminated glass beam and supporter. In peridynamic simulation, the model is discretized to 23608 material points with a grid spacing of 1.5 mm. The time step is 4.0×10-8 s, and the horizon size is 4.5 mm. Velocity of the impactor is 3.13 m/s. Other material parameters keep the same as listed in Table 1.
v=3.13 m/s m=0.1 kg
(10+4+10) mm
10 mm
10 mm
10 mm
200 mm
Fig. 15. The geometric model of the laminated glass beam.
The dynamic fracture process of the laminated glass beam by using the proposed peridynamic approach and stress contours and crack propagation results from finite element analysis on the same model in reference [31] are shown in Fig. 16. The results show that due to the accumulation of the compressive waves in top-layer glass, the crack occurs first from the bottom of the top-layer glass, because the accumulated tensile stress there reaches the tensile strength of glass firstly. Then, the crack propagates upwardly until reaching the boundary of the top-layer glass. In peridynmaic simulation, the first crack occurs at 52 μs, and penetrates through the top-layer glass at 110 μs. While in finite element analysis, the crack initiates at about 59 μs propagates to the top surface at 133μs. In peridynamic simulation, same as experimental observations, the compressive waves will pass through the PVB film, and arrive at the bottom-layer glass. The crack in the bottom-layer glass occurs at about 175 μs, and eventually penetrates through the bottom-layer glass at 280 μs. In the present low-velocity impact loading cases, the principal crack pattern in the laminated glass beam is vertical crack in both top and
bottom glasses, and delamination between glass and PVB does not occur. Fig. 17 shows the corresponding experimental observations on crack initiation and propagation in the top-layer glass of the laminated glass beam with a time interval of 20 μs. It evidently declares that the crack initiates from the bottom side of the top-layer glass and propagates upwards to the top surface of the beam. By comparison, we can find again that the crack morphology, initial locations and propagation obtained by using the proposed peridynamic approach keep good agreement with the results of corresponding experiment and FEM results.
(a1)
(a2)
(b1)
(b2)
(c1)
(c2)
Fig. 16. Dynamic fracture process of a laminated glass beam subjected to impact loading in peridynmaic simulation (left), and finite element simulation (right). (a1). 52 μs, (b1). 60 μs, (c1). 205 μs; (a2). 55 μs, (b2). 65 μs, (c2). 234 μs.
Fig. 17. Experimental observations: crack initiation and propagation in the top-layer glass (with time interval of 20 μs) [31].
After qualitative and quantitative verification, five cases with different PVB interlayer thicknesses: 2 mm, 3 mm, 4mm, 5 mm, and 6 mm were performed to investigate the effect of the thickness of the PVB interlayer on impact failure of the laminated glass. The final fracture modes of the specimen are shown in Fig. 18. The final crack patterns of five cases are almost identical irrespective of the thickness of the PVB interlayer, thereby indicating that the variation of the thicknesses of the PVB interlayer takes little effect on the fracture mode of the glass under low-velocity impact.
(a)
(b)
(c)
(d)
(e) Fig. 18. Damage contours in the laminated glass beam with different PVB interlayer thickness (a): 2 mm, (b): 3 mm, (c): 4 mm, (d): 5 mm, (e): 6 mm.
To further investigate the effect of the thickness of the PVB interlayer, the plots of time variation versus crack initiation and growth in the bottom-layer glass with different PVB interlayer thicknesses are shown in Fig. 19. It indicates that the PVB thickness significantly influences the crack initiation time and has a slight influence on crack propagation velocity in the bottom-layer glass. The crack in the bottom-layer glass with a thinner PVB layer will initiate much earlier, demonstrating that the speed with which compressive waves pass through the PVB interlayer decreases when the PVB interlayer thickness increases, which results in the differences of crack initiation time in bottom-layer glass. However, in bottom-layer glass, the influence on the crack propagation velocity is much slighter. Fig. 19 shows that before the crack length reaches 8 mm, the crack propagation velocity with different PVB interlayer thicknesses is almost identical. After that, the crack propagates a little faster with
thinner PVB interlayer. It is worth to note that a special phenomenon can be found in Fig. 19. When the crack length reaches about 7 mm, the speed of crack propagation will decrease, no matter the thickness of the PVB interlayer changes.
360
2 mm 3 mm 4 mm 5 mm 6 mm
330 300
Time (s)
270 240 210 180 150 120 0
2
4
6
8
10
12
Crack length (mm)
Fig. 19. Plots of time versus crack initiation and growth with different PVB interlayer thickness.
Finally, five cases are conducted to analyze the effect of the thickness of glass layer. For the sake of simplicity, the thickness of the bottom-layer glass keeps 10 mm and the thickness of the top-layer glass changes from 6 mm, 8 mm, 10 mm, 12 mm to 14 mm. The thickness of the PVB interlayer keeps 4 mm, and the impactor velocity is 3.13 m/s. The final fracture of the laminated glass beam and the time variation versus crack length are plotted in Fig. 20 and Fig. 21, respectively.
(a)
(b)
(c)
(d)
(e) Fig. 20. Damage contours in the laminated glass beam with different top-layer glass thickness, (a): 6 mm, (b): 8 mm, (c): 10 mm, (d): 12 mm, (e): 14 mm.
360
6 mm 8 mm 10 mm 12 mm 14 mm
330 300
Time (s)
270 240 210 180 150 120 90 0
2
4
6
8
10
12
Crack length(mm)
Fig. 21. Plots of time variation versus crack length with different top-layer glass thickness.
As shown in Fig. 20, the final fracture modes for the former four cases are nearly identical to those previously. In the last case, the fracture of the bottom-layer glass is different from that of others. Because of the ultra-high value of the thickness of the top-layer glass (14mm), the crack in bottom-layer glass stops propagating at a length of approximately 3.3 mm. With respect to other four cases, the thickness of the top-layer glass affects the crack initiation time obviously in Fig. 21. The crack in the bottom-layer glass with a thinner top-layer glass initiates much earlier, which is similar to the influence of the thickness of the PVB interlayer. However, the crack propagation speed keeps almost unchanged when the thickness of the top-layer glass increases. Similarly, just like the observations mentioned above, the crack propagation speed will significantly decrease when the crack length reaches about 7 mm in call cases.
4. Summary
In this study, the dynamic fracture behavior of laminated glass composites under low-velocity impact is modeled under the framework of the state-based peridynamic theory for the first time. The relaxation response of the PVB interlayer is described by using a peridynamic reformulation of viscoelasticity. A penalty-based method is employed to characterize the adhesion between the glass plies and PVB interlayer, which is of key importance to the mechanical performance of laminated glass. The proposed model and numerical approach were validated through simulating the failure of a laminated glass plate under drop-weight impact loading, and comparing the
numerical results with available experimental observations. After model validation, the effects of the PVB interlayer thickness and glass layer thickness on glass-ply cracking were investigated through the simulation of a laminated glass beam. We compare the crack morphology, initial locations and propagation sequence with experiment results and finite element model results. In various cases, the glass layer keeps the same failure mode. Nevertheless, both the PVB thickness and glass thickness take effects on the crack initiation time in the bottom-layer glass. The crack in bottom-layer glass with thinner PVB interlayer or thinner top-layer glass will initiate significantly earlier. Furthermore, the speed of the crack propagation keeps nearly unchanged when the PVB thicknesses and glass thicknesses increase. In all simulation results for laminated glass with different PVB interlayer thicknesses or glass layer thicknesses, it shows that the crack propagation speed will significantly decrease when the crack length reaches some critical value (e.g. 7 mm under the working cases in this study). Further studies will focus on the mechanism of the crack growth.
Acknowledgements: The authors acknowledge the support from the Fundamental Research Funds for the Central Universities in China (No. 2017B13014, 2019B66014), the Research & Innovation Project of Jiangsu Province (No. SJKY19_0419), the National Natural Science Foundation of China (No. 51679077, 11932006), and the National Key R&D Program of China (No. 2018YFC0406703).
References [1]
Muralidhar S, Jagota A, Bennison SJ, Saigal S. Mechanical behaviour in tension of cracked glass bridged by an elastomeric ligament. Acta Mater 2000, 48(18-19):4577-88.
[2]
Xu J, Li YB, Liu BH, Zhu MY, Ge DY. Experimental study on mechanical behavior of PVB laminated glass under quasi-static and dynamic loadings. Compos B 2011;42(2):302-8.
[3]
Chen JJ, Xu J, Liu BH, Yao XF, Li YB. Quantity effect of radial cracks on the
cracking propagation behavior and the crack morphology. PLOS ONE
2014;9(7): e98196. [4]
[5]
Chen JJ, Xu J, Yao XF, Xu XQ, Liu BH, Li YB. Different driving mechanisms
of in-plane cracking on two brittle layers of laminated glass. Int
J Impact Eng
2014; 69:80-5.
Zhang XH, Hao H, Ma G. Laboratory test and numerical simulation of laminated glass window vulnerability to debris impact. Int J Impact Eng 2013; 55:49-62.
[6]
Pyttel T, Liebertz H, Cai J. Failure criterion for laminated glass under impact loading and its application in fnite element simulation, Int J Impact Eng 2011;38 (4):252-63.
[7]
Hidallana-Gamage HD, Thambiratnam DP, Perera NJ. Failure analysis of laminated glass panels subjected to blast loads. Eng Fail Anal 2014;36:1429.
[8]
Peng Y, Yang JK, Deck C, Willinger R. Finite element modeling of crash test
behavior
for
windshield
laminated
glass,
Int
J
Impact
Eng
2013;57:27-35. [9]
Xu J, Li YB. Study of damage in windshield glazing subject to impact by a pedestrian’s head, Proc Inst Mech Eng D 2009;223(D1):77-84.
[10]
Zang MY, Lei Z, Wang SF. Investigation of impact fracture behavior of automobile laminated glass by 3D discrete element method. Comput Mech
2007;41 (1):73-83. [11]
Chen SH, Zang MY, Wang D, Yoshimura S, Yamada T. Numerical analysis of impact failure of automotive laminated glass: A review. Compos Pt B-Eng 2017; 122:47-60.
[12]
Xu J, Li YB, Chen X, Yan Y, Ge DY, Zhu MY, et al. Characteristics of windshield cracking upon low-speed impact: numerical simulation based on the extended finite element method, Comput Mater Sci 2010;48(3):582-88.
[13]
Chen SH, Zang MY, Xu W. A three-dimensional computational framework for
impact fracture analysis of automotive laminated glass. Comput Methods
Appl Mech Eng 2015; 294:72-99. [14]
Silling SA. Reformulation of elasticity theory for discontinuities and longrange forces. J Mech Phys Solids 2000;48(1):175-209.
[15]
Silling SA, Askari E. A meshfree method based on the peridynamic model of solid mechanics. Comput Struct 2005;83(17-18):1526-35.
[16]
Breitenfeld MS, Geubelle PH, Weckner O, Silling SA. Non-ordinary statebased peridynamic analysis of stationary crack problems. Comput Method Appl Mech Eng 2014; 272:233-250.
[17]
Silling SA. Dynamic fracture modeling with a meshfree peridynamic code. In: Proceedings Second MIT Conference on Computational Fluid and Solid Mechanics, 17-20 June 2003.P. 641-44.
[18]
Kilic B, Madenci E. Prediction of crack paths in a quenched glass plate by using peridynamic theory. Int J Fract 2009;156(2):165-77.
[19]
Ha YD, Bobaru F. Studies of dynamic crack propagation and crack branching with Peridynamics. Int J Fract 2010;162(1-2):229-44.
[20]
Hu YL, Madenci E. Peridynamics for fatigue life and residual strength prediction of composite laminates. Compos Struct 2017; 160:169-84.
[21]
Huang D, Zhang Q, Qiao PZ. Damage and progressive failure of concrete structures using non-local peridynamic modeling. Sci China Technol Sci; 2011;54(3):591-6.
[22]
Xu J, Askari A, Weckner O, Silling SA. Peridynamic analysis of impact
damage in composite laminates. J Aerospace Eng 2008;21(3):187-94. [23]
Askari A, Nelson K, Weckner O, Xu JF, Silling SA. Hail impact characteristics of a hybrid material by advanced analysis techniques and testing. J Aerospace Eng 2011;24(2):210-7.
[24]
Azizi M, bin Mohd Ihsan A, Mohamed N. The peridynamic model of viscoelastic creep and recovery. Multidiscipline Modeling in Materials and Structures 2015; 11(4):579-97.
[25]
Mitchell JA. A non-local, ordinary-state-based viscoelasticity model for peridynamics. Sandia National Lab Report, 2011.
[26]
Delorme R, Tabiai I, Lebel L, Lévesque M, Generalization of the ordinary state-based peridynamic model for isotropic linear viscoelasticity. Mech Time-Depend Mater 2017; 21:549-75.
[27]
Dorduncu M, Barut A, Madenci E. Peridynamic truss element for viscoelastic
deformation,
57th
AIAA/ASCE/AHS/ASC
Structures,
Structural Dynamics, and Materials Conference, AIAA SciTech. [28]
Madenci
E,
Oterkus
S.
Ordinary
state-based
peridynamics
for
thermoviscoelastic deformation. Eng Fract Mech 2017; 175:31-45. [29]
Weckner O, Mohamed N. Viscoelastic material models in peridynamics, Appl Math Comput 2013; 219:6039-43
[30]
Kilic B, Madenci E. Structural stability and failure analysis using peridynamic theory. Int J Non-Linear Mech 2009; 44(8):845-54.
[31]
Chen SH, Zang MY, Wang D, Zheng ZM, Zhao CL. Finite element modelling
of impact damage in polyvinyl butyral laminated glass. Compos
Struct 2016; 138:1-11. [32]
Zang MY, Lei Z. An approach to combining 3D discrete and finite element methods based on penalty function method. Comput Mech 2010;46(4):60919.