Modelling impact fracture and fragmentation of laminated glass using the combined finite-discrete element method

International Journal of Impact Engineering 112 (2018) 15–29

Contents lists available at ScienceDirect

International Journal of Impact Engineering journal homepage: www.elsevier.com/locate/ijimpeng

Modelling impact fracture and fragmentation of laminated glass using the combined finite-discrete element method

MARK

Xudong Chena, Andrew H C Chanb,* a b

School of Civil Engineering, Suzhou University of Science and Technology, Suzhou 215011, China School of Engineering and ICT, University of Tasmania, Hobart 7001, Australia

A R T I C L E I N F O

A B S T R A C T

Keywords: Fracture Fragmentation Laminated glass Impact Combined finite-discrete element method

Fracture and fragmentation responses of laminated glass under hard body impact are modelled with the combined finite-discrete element method (FEM-DEM). The method is essentially a discrete element method with finite element mesh be embedded in, yielding more accurate contact forces as well as mass and energy conservation. Failure models of glass, interlayer and glass-interlayer interface are discussed respectively and proven to be reliable in simulating the rupture of laminated glass. Numerical examples are presented and validated with results from different sources, and the advantage of the FEM-DEM modelling on the impact failure of laminated glass over its parent FEM and DEM is demonstrated. The influences of Young's modulus of interlayer are further discussed, showing that a stiff resin can effectively improve the deformation performance at little expense of energy absorption. In general, modelling the impact failure of laminated glass using the combined finite-discrete element method is successful.

1. Introduction Laminated glass, which is typically composed of glass plies and interlayer(s), is increasingly popular for structural purposes in recent decades [1–2]. Consequently, its fracture and fragmentation responses under impact action attract growing attentions from both academics and practising engineers. Earlier investigations on the performance of laminated glass elements under bending or lateral pressure were conducted experimentally [3–4]. Recently, Zhang et al. [5] presented laboratory test results on the vulnerability of laminated window glass, and concluded that the interlayer thickness plays a crucial role in improving the penetration resistance. Marcon et al. [6] studied the windshield laminated glass experimentally and identified the damage process from the initial micro-crack to catastrophic macroscopic failure. Since high-speed camera is essential for capturing the transient fracture responses in laminated glass, direct experiments are expensive to perform. In addition, tests could yield some random results and be influenced by unpredictable factors. As an alternative, a variety of novel computational approaches are developed to solve complex fracture and fragmentation. Rabczuk and Belytschko proposed two and three-dimensional cracking particles method [7–9] within the elementfree Galerkin (EFG) framework. The method is capable of handling crack branching and fragmentation, and results were validated with numerical and experimental data. Meanwhile, peridynamics was

*

employed to model the rupture of composite laminates [10] under shock-type loading. Further, dual-horizon peridynamics (DH-PD) was developed by Ren et al. [11–12] and used to examine the crack propagation in Kalthoff–Winkler plate as well as crack branching problems.. As for the impact fracture and fragmentation of laminated glass, the finite element method (FEM) is commonly used. In Flocker and Dharani [13–14], linear elastic constitutive relationship was employed to evaluate the fracture of glass plies subject to low velocity missiles based on the maximum tensile stress criteria using an explicit finite element code DYNA2D. They assumed that the interlayer was made out of Polyvinyl Butyral (PVB) and a Hertzian cone was obtained. The laminated glass was assumed to be perfectly bonded without any potential debonding or slipping during impact and only the glass ply in contact with the projectile can be damaged. Du Bois et al. [15] applied two coincided elements (a Belytschko–Tsai shell element for glass and a membrane element for PVB interlayer) in their laminated glass impact investigation. A piecewise linear plastic fracture criterion that allows for a small plastic strain is used for glass failure, and was employed to predict the behaviour of windscreen under a sphere impact. However, the use of coincidental elements assumes that both outer and inner glass plies fail at the same time, which may be unrealistic. Pyttel et al. [16] modelled the crack initiation and propagation in windshield laminated glass subjects to rigid ball impact based on the critical energy threshold and

Corresponding author. E-mail addresses: [email protected] (X. Chen), [email protected] (A.H.C. Chan).

http://dx.doi.org/10.1016/j.ijimpeng.2017.10.007 Received 2 April 2017; Received in revised form 1 October 2017; Accepted 12 October 2017 Available online 14 October 2017 0734-743X/ © 2017 Elsevier Ltd. All rights reserved.

International Journal of Impact Engineering 112 (2018) 15–29

X. Chen, A.H.C. Chan

local Rankine (maximum stress) criteria using the commercial finite element solver PAM-CRASH. Recently, Chen et al. [17] developed a linear contact algorithm LC-Grid and cohesive model, and coupled them with FEM modelling on impact fracture of laminated glass. It is understood that although the FEM is capable of predicting the crack initiation and propagation in the laminated glass with reasonable accuracy, possible fragmentation and interactions between fragments are difficult to be simulated due to the natural limitations, i.e. the continuity requirement of the FEM and this cannot be readily resolved by using enhanced shape functions such as XFEM. For the impact fracture and fragmentation analysis of laminated glass, the discrete element method (DEM) can better accommodate the inherent discontinua mechanism. Oda et al. [18] employed the DEM to model the impact behaviour of laminated glass. Later, the approach was extended to laminated glass of bi-layer type [19]. Fragmentation of laminated glass subjects to ball impact was investigated by Zang et al. [20] using rigid spherical discrete elements and a Mohr-Coulomb failure criterion. In the work of Gao [21], spherical discrete elements were also used and qualitative analysis on fracture of laminated glass is performed. Xu and Zang [22] extended the sub-domain concept and attempted to employ the spherical discrete elements for potential fracture region and cubical finite elements for non-fracture region with four-point combined DE/FE algorithm to evaluate the interactions on the DE-FE interfaces. An extrinsic cohesive fracture model is employed for the brittle fracture of glass, while the interlayer is assumed not to fracture during the impact process. Though spherical discrete elements are capable of simulating the post-damage behaviour, crack patterns are usually of low precision and only some non-recognisable cracks can be obtained. Generally, the DEM is good at simulating the fragments and fragmentations of laminated glass under impact, however, precision in predicting crack initiation and fracture patterns is usually unsatisfactory. To obtain accurate crack initiation and propagation, realistic fracture patterns and practical fragmentation predictions, a novel computational method, i.e. the combined finite-discrete element method (FEM-DEM) [23] is employed in this paper. It overcomes the limitations of the FEM and DEM and best suits for transient computation of discontinua where contact is involved. In the analysis, the structure is fully discretised into a number of individual discrete elements, and each discrete element interacts with those that are in contact with it. Within a discrete element, a finite element mesh is formulated, leading to a more accurate estimate of contact forces, stress distribution and deformation. Fracture models are implemented in the joint interfaces between two adjacent elements. By using the FEM-DEM with an appropriate fracture model, the fracture process of laminated glass under impact can be better represented than its parent FEM and DEM. Unlike some FEM approaches, fragments in the FEM-DEM simulation generated by the impact also take part in the following contact interactions, making the mass and energy comply with the law of conservation. The combined FEM-DEM was developed by Munjiza and co-workers in 1990s [24–25]. A Munjiza-No Binary Search (NBS) algorithm [26] was proposed by Munjiza and Andrews for contact detection. Later, a combined single and smeared crack model [27] was developed by Munjiza et al. and implemented into the FEM-DEM program “Y” [28]. Issues of mesh sensitivity were discussed and addressed by Munjiza and John [29]. These principles and computational points of the FEM-DEM were systematically summarised in a monograph [23]. The FEM-DEM is still developing and a recent paper [30] proposed a shell FEM-DEM element for the fracture and fragmentation of thin shells. Furthermore, a parallelization solution based on space decomposition is also introduced by Lukas et al. [31]. A newly published book [32] provided ‘down to earth’ illustration of the framework for developing large-strain based nonlinear material laws and the large-strain large-displacement based FEM-DEM. Though there are some applications of the FEM-DEM on fracture considerations, most of them are emphasised on the static or dynamic behaviour of brittle or quasi-brittle civil engineering materials like rock falling [33–34], concrete breakage [35–36] and masonry

failure [37–39]. Little attention is given to the impact fracture and fragmentation of glass and laminated glass. The present authors [40] successfully simulated the impact fracture of monolithic glass using the combined FEM-DEM, with the acquisition of a variety of damage modes, e.g. flexural crack, Hertzian cone crack, punching failure and etc. Impact fracture and fragmentation process of laminated glass is highly nonlinear and discontinuous, thus the discontinua-based FEMDEM approach is appropriate. As was indicated by Rabczuk and Belytschko [7], “for problems with many cracks and less predictability, such as fragmentation, simple methods such as interelement separation may be appropriate.” Since many cracks occur in laminated glass within the impact duration and inter-element separation is employed in the FEM-DEM, the FEM-DEM is suitable for the impact fracture and fragmentation modelling of laminated glass. Some simple computational issues on the FEM-DEM modelling of laminated glass were addressed in [41–42]. This paper is a substantial extension of [40] and emphasised on the modelling of impact fracture and fragmentation of laminated glass. Following the theoretical background of FEM-DEM, the fracture criteria of glass, interlayer and glassinterlayer interface are discussed respectively. Fracture and fragmentation of laminated glass under hard body impact is simulated using the FEM-DEM program. Results are proven to be reliable by comparing with those from existing literature. A further discussion is presented on the effect of the Young's modulus of interlayer, quantitatively revealing the role of interlayer in the impact fracture and fragmentation mechanism of laminated glass. 2. The FEM-DEM The combined finite-discrete element method (FEM-DEM) aims at solving large scale transient dynamics where contact and fracture are involved. It considers the solid a combination of both continua and discontinua. Details of the FEM-DEM are found in [23], and some key aspects of the FEM-DEM are addressed in this section. 2.1. Motion of elements In the FEM-DEM, the stress and deformation are formulated by the standard FEM approach, while the motion and contact forces of elements are considered by the DEM. The translational and rotational motions of a single discrete element i are expressed in terms of Newton's second law as

mi Ii

d2 ri = Fi dt 2

d ωi = Ti dt

(1) (2)

where mi is the mass of discrete element i; ri is the position; Ii is the moment of inertia; ωi is the angular velocity; Fi and Ti are net external force and torque, respectively. Accordingly, velocity and position of an arbitrary discrete element can be determined explicitly at each time step. 2.2. Contact algorithm Contact algorithm, which can be classified as contact detection and contact interaction, is the central of the FEM-DEM. Contact detection is a particular research domain of DEMs and beyond the scope of this paper. The Munjiza-No Binary Search (NBS) algorithm is employed in the current FEM-DEM so that computational time is linearly proportional to the number of elements involved in the simulation. Full details on the Munjiza-NBS contact detection algorithm can be referred to [23,26].Once contact detection is completed, contact interaction follows where contact forces are defined according to interaction law. In the FEM-DEM, distributed contact force is adopted for two discrete 16

International Journal of Impact Engineering 112 (2018) 15–29

X. Chen, A.H.C. Chan

S Pt,Pc

df

Fig. 1. Contact force in the 2D FEM-DEM.

elements in contact. For the pair of discrete elements in contact, one is denoted as contactor, and the other is target. In 2D, contact forces are evaluated from the overlapping area S. Referring to Fig. 1, penetration of any elemental area dA of the contactor into the target results in an infinitesimal contact force

d f = −d ft + d f c

(a)

(b)

(c)

(d)

(3)

where the subscripts t and c represent the target and contactor, respectively. In Eq. (3), dft and dfc can be written in the form of

d ft = −Ep gradϕc (Pc) dA

(4)

d f c = −Ep gradϕt (P) t dA

(5)

where Pc and Pt are the points sharing the same coordinate on S; ϕc and ϕt are potentials; Ep is a contact penalty parameter, which is usually a large value; grad is the gradient. With the integration of elemental area dA over the whole S, contact force f is obtained as

f = Ep

dA ∫S [grad ϕc (Pc) − grad ϕt (P)] t

Fig. 2. Joint interfaces in the FEM-DEM: (a) 2D initial; (b) 2D deformed; (c) 3D initial; (d) 3D deformed.

(6)

By replacing the overlapping area S with overlapping volume V, Eq. (6) is applicable to three-dimensional problems. The contact penalty parameter Ep plays an important role in evaluating contact forces and controlling displacement errors. Ep cannot be too small, otherwise noticeable element penetration could occur, which is not realistic [23]. A large contact penalty guarantees effective contact. Referring to [28], Ep is usually 2–100 times greater than λ (λ is the Eν first Lame constant, λ = (1 + ν)(1 − 2ν) where E is the Young's modulus and ν is the Poisson's ratio). In glass impact simulations, the value of Ep depends on glass and projectile material properties.

2.4. Multiple-material modelling The FEM-DEM is capable of modelling multiple materials within one entity, providing great advantage and convenience in laminated glass simulation. Should there are three different materials A, B and C within one entity (Fig. 3), joint properties exists in the interfaces between different materials. The joint properties can be automatically identified by the FEM-DEM, and all that users need to do is to appropriately define them. Table 1 shows the material and joint properties definitions recognised by the FEM-DEM. In Table 1, PA, PB and PC are material properties of A, B and C, respectively. Data in the third column are joint properties for the interfaces between different materials. For example, JPAB is the joint property of the interface between material A and B. Similar definitions apply to JPAC and JPBC. Users only need to define JPAB, JPAC and JPBC, and the FEM-DEM program will detect the interfaces by itself. The capability of modelling multiple materials within one entity makes the modelling of laminated glass possible. Fig. 4 shows a typical three-layer laminated glass. It is composed of two glass plies (one is called the outside ply and the other is named the inside ply) and one interlayer. The thickness of the outside glass ply, inside glass ply and the interlayer are ho, hi and hinter, respectively. Accordingly, there are two interfaces between both the outside and inside glass plies and the interlayer. By simply defining the joint properties, laminated glass can

2.3. Joint interface In the current FEM-DEM, only three-node constant strain triangular elements and four-node constant strain tetrahedral elements are available, saving considerable computational time for contact detection and interaction. Since contact algorithm in the FEM-DEM is quite time consuming, the employment of the simplest elements enable the improvement of computational efficiency [23]. Although the implemented elements are of the simplest type, triangle or tetrahedral elements can be applied to any conceivable boundary with complex curvatures, and satisfactory results can still be achieved by increasing the number of elements in desired area. Cracks are assumed to occur coincide with element edges so that re-meshing can be avoided. Between each pair of discrete elements in contact, linked interfaces (Fig. 2) are employed to evaluate the material from continue to discontinua. In Fig. 2, both initial and deformed interfaces in 2D and 3D are presented, where the dotted lines in Fig. 2(b) and (d) represent the initial position of the interface. Interfaces can break according to the deformed separation δ between the two adjacent element boundaries, and the change of δ follows a suitable fracture law based on the nature of examined material. Fracture models are implemented in these interfaces. Different materials can be modelled separately by incorporating appropriate fracture criteria. Though elements can only break along their edges, possible mesh bias can be reduced to the minimum if unstructured and small enough elements are employed [29].

A C B

Fig. 3. Multiple materials within one entity.

17

International Journal of Impact Engineering 112 (2018) 15–29

X. Chen, A.H.C. Chan

1.0

Material property

Material property

Joint property

PA PA PB

PB PC PC

JPAB JPAC JPBC

Outside glass ply Interlayer

Normalised bonding stress - z

Table 1 Material and joint properties.

ho

0.8 a=1.2, b=-1.0, c=1.0

0.6 0.4 0.2

hinter 0.0

Inside glass ply

0.0

hi

Fig. 4. Portion of a typical three-layer laminated glass.

⎧⎡2 δ − δ ⎪⎣ ⎪⎢ p σn = f z ⎨ t ⎪2 δ f ⎪ δp t ⎩

δ δp

t

0 ≤ δ ≤ δp δp < δ ≤ δc δ<0

(8)

(9)

In Eq. (9), D is a fracture damage index within the interval of [0, 1]. If δp < δ ≤ δc, D = (δ − δp)/(δc − δp) ; when δ > δc, D = 1, leading to the total fracture and free movement of two adjacent elements. δp is calculated from the elastic modulus and the tensile strength while δc is obtained from the fracture energy. It is verified in [40] that for the fracture of glass, a = 1.2, b = −1.0 and c = 1.0, and a normalised glass softening curve is plotted in Fig. 5. A full introduction of the glass fracture model can be referred to [40].

3.2. Interlayer Most interlayer materials commonly used are PVB-based. Though the PVB material has been employed in the laminated windshield as interlayer for some decades, research on its mechanical properties is still limited and not well documented. However, it is widely held that PVB is largely linear viscoelastic when it is loaded slowly, while it behaves elasto-plastic or even brittle under high-speed tension [44]. It is considered nearly incompressible and the mechanical properties are temperature and strain rate sensitive. For the sake of simplicity, the interlayer material is implemented without the consideration of thermal effect and loading rate in this study. And since the investigated impact velocities are not very low and loading time is very short, viscous effect can be ignored and the interlayer can be considered as pure elastic before it reaches the tensile strength. The softening behaviour can be represented by a relatively flat descending curve. Hence, the cohesive crack model in Section 3.1 is extended for the interlayer. As the interlayer also exhibit some ductility, a large fracture energy is usually assigned so that the curve shape in Fig. 6 is considered to be reasonable and acceptable. According to Tvergaard and Hutchinson [45], the shape of the cohesive strain softening curve plays a less significant role than the cohesive fracture energy (Gf) and the maximum cohesive traction (ft) when simulating the fracture of materials with ductility. Consequently, a normalised curve (Fig. 6) with a = 5.5, b = −5.0 and c = 1.0 (a, b and c are as defined in Section 3.1) is assumed to represent the softening behaviour of interlayer when the bonding stress exceeds the tensile strength of the material.

3.1. Glass There are evidences [43] showing that the dominant mechanism of brittle fracture in glass is Mode I cracking, therefore a cohesive crack model in [40] is adopted here for glass. It considers a crack develops when the bonding stress σ reaches the tensile strength ft, which is the onset of the damage. After this point, σ is assumed to decrease gradually as the crack width δ increases. When δ reaches a critical value δc, σ drops to zero and the crack is assumed to be completely open. In the FEM-DEM, standard FEM approach is used before the material reaches its tensile strength. The fracture developing process is captured through a weak constitutive relation of the material, enabling the cracks being smeared out over the continuum. This stage is named with softening. The softening is depicted through a descending σ-δ curve, where the area under the curve equals to the fracture energy Gf as is in Eq. (7).

σ dδ

1.0

2

( ) ⎤⎦⎥ f

a + b − 1 D (a + bc /((a + b)(1 − a − b))) ⎤ z = ⎡1 − e [a (1 − D) + b (1 − D)c ] a+b ⎣ ⎦

This section briefly discusses the fracture models used in the impact fracture and fragmentation analysis of laminated glass. The rupture of glass, interlayer and the glass-interlayer interface are introduced respectively, suggesting the damage of glass, interlayer and their interface are all considered in the FEM-DEM modelling. It should be noted that potential failures of interlayer and interface are predefined, however, these failure criteria may or may not be used in FEM-DEM simulations due to the specification of the examined example.

δc

0.8

where z is a heuristic parameter

3. Fracture models

∫0

0.4 0.6 Normalised crack index - D

Fig. 5. Normalised bonding stress z = σ / ft versus the crack index D for glass.

be automatically assembled and genuinely modelled. Apart from the above aspects, damping is also accounted for in the FEM-DEM to better accommodate energy dissipation. According to [23], any energy dissipation in contact can be attributed to friction or plastic straining of surface asperities, and the latter is approximated by a viscous damping model. The physical meaning of this contact damping in the FEM-DEM is plastic deformation and/or breaking of surface asperities. According to [28], a viscous damping with the value of 2h Eρ is used, where h is the characteristic size of the smallest element, E is the Young's modulus and ρ is the material density. Besides contact damping, material damping due to deformation is naturally covered by discrete elements in the FEM-DEM when they are discretised into finite elements [23].

Gf =

0.2

(7)

Since the assumption of no pre-existing crack is held, crack initiation and propagation are determined directly by the actual stress and deformation distributions within the material. Note that δp is the onset separation, which equals to the value when σ = ft. δc is the critical separation when σ reduces to zero. The complete relation between σ and δ is given in Eq. (8) [27]. 18

International Journal of Impact Engineering 112 (2018) 15–29

X. Chen, A.H.C. Chan

Normalised bonding stress - z

1.0

interlayer interfacial behaviour can be characterised. 4. Numerical examples

0.8

In this section, numerical examples on the impact fracture and fragmentation of laminated glass are presented. In Section 4.1, elastic impact behaviour without any fracture in laminated glass body is studied; in Section 4.2, an example with laminated glass beam under blunt impact producing a cone type crack is investigated and verified; in Section 4.3, fragmentation and post-damage failure of laminated glass is examined and compared with results from DEM; in Section 4.4, a laminated glass beam subjects to a rectangular projectile is simulated, and results are compared with high speed camera images. These examples are such chosen that they cover the range from elastic to postfailure, and are well compared with data from both computational and experimental sources.

0.6

0.4

0.2

a=5.5, b=-5.0, c=1.0

0.0 0.0

0.2

0.4 0.6 Normalised crack index - D

0.8

1.0

Fig. 6. Normalised bonding stress z = σ / ft versus the crack index D for interlayer.

4.1. Elastic impact

glass element

Fig. 8(a) shows a 1 m long 22 mm high laminated glass beam with both ends clamped in fixed channels. ho = hi = 10 mm and hinter = 2 mm. The centre of the beam subjects to the impact of a 10 mm radius steel ball. A triangular mesh with 9180 elements is employed for the laminated glass beam. Elements with a characteristic size of 1 mm are used in the impact effective area (Fig. 8(b)), while coarse mesh (characteristic size of 2 mm) is adopted in the far field taking advantage of the Saint–Venant's principle. This mesh is defined as mesh A as comparison will be made with results from other mesh later. The boundary of the ball is also meshed with small elements to make sure a reasonably circular shape can be guaranteed. Point A, which situates at the bottom middle of the inside glass ply, is selected. Point B locates in the bottom middle of the interlayer and will be used in a later section. Time step used for this analysis is 0.2 ns. To simulate the elastic behaviour, very low impact velocity v of the projectile is applied as v = 0.2, 0.3 and 0.4 m/s. It is assumed that the interlayer is perfect elastic and the glass-interlayer interface is strong enough to resist potential debonding. For typical glass, the density ρ = 2500 kg/m3, Young's modulus E = 70 GPa and Poisson's ratio ν= 0.2. The tensile strength of glass varies according to different types of source. The present authors [40] suggest that 30 MPa is a reasonable value for the tensile strength of glass, which was supported by a variety of literature like Sedlacek et al. [47], Overend and Zammit [48] and Ballarini et al. [49]. For the fracture energy of glass, a value of 4.0 J/m2 is widely accepted. The interlayer is assumed to be PVB material. Detailed material properties used in the FEM-DEM simulation are tabulated in Table 2. From Table 2, λglass = 1.944 × 1010 Pa and λprojectile = 1.154 × 1011 Pa. The contact penalty Ep is 1.0 × 1012 Pa, which is 8.67 times of the max{λglass, λprojectile} and satisfies the requirement in [28]. Same impact problems are solved with the ABAQUS package with

interface

interlayer element Fig. 7. Interfaces between glass and interlayer elements.

3.3. Glass-interlayer interface Fig. 7 shows the interface between glass and interlayer elements. Physically, the failure behaviour of glass-interlayer interface under impact is complicated, and there is no specific value for the adhesive properties in international or national standards as well as existing literature. To characterise the adhesive properties of the glass-interlayer interface and investigate the potential debonding, the bilinear descending softening relationship reported in Fig. 5 is employed. The main reasons of using that model are: i) Similar cohesive zone fracture models have been employed in existing literature [46] for glass-interlayer interface; ii) In a single impact, the laminated glass mainly experiences Mode I failure so that the Mode I-based fracture model in Section 3.1 is a reasonable assumption; iii) Different values of ft and Gf can be assigned to the glass-interlayer interface, making the model versatile for future calibration with potential experimental data. In general, with different combinations of ft and Gf, the glass-

Fig. 8. Configurations of the elastic impact problem: (a) overall view; (b) detailed mesh in the impact effective area.

19

International Journal of Impact Engineering 112 (2018) 15–29

X. Chen, A.H.C. Chan

velocities retain constant value until around t = 0.6 ms; stage III is from t = 0.6 ms till around t = 0.9 ms, where velocities are still constant but either slower than their counterpart in stage II (v0 = 0.2 and 0.3 m/s cases) or getting positive as the projectiles are bounced back (v0 = 0.4 m/s case); stage IV is from t = 0.9 ms until the end of the simulation where all the projectiles are bounced back. The transit time bands are marked in Fig. 10. Should Figs. 9 and 10 are compared together, one can find that the stages in Fig. 10 largely correspond to different stretches of the deformation curves in Fig. 9. The slope of the deformation curve changes from t = 0.1 ms, which roughly corresponds to the first transit band; another major change on the deformation curve occurs around t = 0.6 ms, as the deformation curvatures become non-regular after that time point. Since there is no contact between the projectile and the laminated glass as the ball is bounced back for v0 = 0.4 m/s case, the laminated glass enters a free vibration phase leading to the deformation oscillation. Though in the v0 = 0.2 and 0.3 m/s cases, the residual velocity of the projectile is still negative and there is contact between the projectile and laminated glass, the actual velocity is as low as −0.0145 and −0.0137 m/s, which is only 7.25% and 4.57% of the initial impact velocities and almost can be neglected. The projectile can be considered “adhered” to the laminated glass and the whole system enters the free vibration phase, resulting in noticeable oscillation in Fig. 9. And after t = 0.9 ms, residual velocities of the projectiles for v0 = 0.2 and 0.3 m/s cases also get positive and the laminated glass beams enter a complete free vibration phase. It should be noted that residual velocities in v0 = 0.2 and 0.3 m/s cases are higher than that in the v0 = 0.4 m/s case, as the projectiles in the former two cases are pushed by the large amplitudes of the beam oscillation. It is worth mentioning that due to the existence of interlayer and the bouncing back of the projectiles, deformation of laminated glass is sophisticated and highly non-linear even when the impact velocity is low, and these characteristics have been well captured by the FEM-DEM.

Table 2 Material properties used in the elastic impact analysis.

Density (kg/m3) Young's modulus (GPa) Poisson's ratio Tensile strength (MPa) Fracture energy (J/m2)

Glass

PVB

Projectile/support

2500.0 70.0 0.2 30.0 4.0

100.0 0.1 0.4 –– ––

7800.0 200.0 0.3 –– ––

Y-Displacement (mm)

0 -0.01 -0.02 -0.03 -0.04 -0.05 0

0.2

0.4

0.6

0.8

1

Time (ms) 0.2m/s FEM-DEM

0.2m/s ABAQUS

0.3m/s FEM-DEM

0.3m/s ABAQUS

0.4m/s FEM-DEM

0.4m/s ABAQUS

Fig. 9. Vertical deformations of point A under different low impact velocities.

the same mesh used in the FEM-DEM and CPS3 element. In ABAQUS analysis, materials are assumed not able to fracture and only elastic impact is considered. Fig. 9 shows the vertical displacement of point A along with time under different impact velocities. It can be observed that the deformation obtained from the FEM-DEM agrees very well with that from ABAQUS. Both curves at the same impact velocity exhibit highly similar trends. When v = 0.2 m/s, displacement curves from the FEM-DEM and ABAQUS reach the closest agreement. With the increase of impact velocity, difference between the FEM-DEM and ABAQUS become slightly obvious but are still within acceptable tolerance (< 5%) except at some individual data points. The small difference can be attributed to the different approaches used to simulate the problem. A further investigation on the vertical velocity of the projectile in Fig. 10 reveals some possible relationship between the projectile and the target on the elastic deformation of the laminated glass. The velocity data is obtained from the FEM-DEM. Referring to Fig. 10, the velocity curves of projectiles can be divided into four stages regardless of its specific initial value. Stage I is from the onset of the impact to t = 0.06–0.08 ms, where the velocities of the projectile decrease dramatically from their initial value to less than 0.05 m/s; in stage II, the

4.2. Convergence In this section, initial impact velocity is raised to v = 20 m/s so that fracture occurs without the perforation. Convergence study on the residual velocity of projectile, total kinetic energy, load-deflection characterisation and crack propagation speed are examined. A different mesh (mesh B as is shown in Fig. 11) is used. Comparing with mesh A, 1 mm characteristic size elements are distributed over a larger area and a total number of 13,080 elements are employed. All other parameters are kept same as those in Section 4.1. The laminated glass fractures quite quickly and simulation up to t = 0.2 ms is enough as most cracks occur by that time. To evaluate the difference between results from mesh A and B, a relative error η is defined as

0.1

η=

Y-velocity (m/s)

0

-0.2 0.2m/s 0.3m/s 0.4m/s

-0.4 0

0.2

0.4 0.6 Time (ms)

0.8

(10)

where ΩA and ΩB are results from mesh A and B, respectively. Fig. 12(a) shows the residual velocity curves of projectile with different mesh, and error analysis in Fig. 12(b) shows that η is smaller than 5%. It can be concluded that residual velocity of projectile is converged against different mesh. Assume the out-of-plane thickness is unit one. Total kinetic energy curves are plotted in Fig. 13(a). It can be observed that both curves are very close to each other. Analysis on η is given in Fig. 13(b), showing that a notably less than 3% difference is achieved and energy convergence with different mesh is demonstrated. According to Newton's third law of motion, the magnitude of impact force onto the laminated glass equals to the force that the projectile subjects to. According to the impulse theorem (FΔt = mΔv), since the time interval Δt is quite small between each time step, impact force (or load) can be approximated. Fig. 14 schematically shows the load-

-0.1

-0.3

ΩA − ΩB × 100% ΩA

1

Fig. 10. Vertical velocities of the projectile under different impacts obtained from the FEM-DEM.

20

International Journal of Impact Engineering 112 (2018) 15–29

X. Chen, A.H.C. Chan

Fig. 11. Configuration of the impact effective area with mesh B.

600

Mesh A

-5

Total Kinetic energy (J)

Residual velocity (m/s)

0

Mesh B

-10 -15 -20

400 300 Mesh A

200

Mesh B

100 0

-25 0

0.05

0.1

0.15

0

0.2

0.05

0.1

Time (ms)

Time (ms)

(a)

(a)

0.15

0.2

0.15

0.2

5.0%

Relative error η

5.0%

2.5%

Relative error η

500

0.0%

2.5%

0.0%

-2.5%

-2.5% -5.0% 0

-5.0% 0

0.05

0.1

0.15

0.05

0.1

Time (ms)

0.2

Time (ms)

(b)

(b)

Fig. 13. Total kinetic energy and error with different mesh: (a) kinetic energy curves; (b) error η.

Fig. 12. Residual velocity of projectile and error analysis: (a) residual velocity; (b) error η.

500

deflection curves for mesh A and B, respectively. The largest value of impact force occurs at the very beginning, and then decrease along with the increase of deflection with oscillation. Obviously, load-deflection curves are in good agreement, they exhibit highly similar development trends and converge well. As there are many cracks within the laminated glass when it subjects to impact, it is not realistic to evaluate crack propagation speed of all cracks. A dominant flexural crack locates near the centre of the inside glass ply is chosen to be examined. Table 3 and Fig. 15 present the crack propagation details. It is noted that the crack length is same for mesh A and B as crack propagates over the same number of element edges in the simulation. The propagation speed can be considered as reasonably converged for a DEM based simulation tool. Convergence on different temporal discretisation can simply be satisfied in the FEM-DEM as long as the time step Δt is smaller than a critical value, and relevant examples can be found in [23,41]. Referring to [23],

Δt ∝ h

ρ E

Mesh A

Impact force (KN)

400

Mesh B

300 200 100 0 0

0.5

1

1.5

2

Deflection (mm) Fig. 14. Load-deflection curves for mesh A and B.

where h is the characteristic size of element. Once a critical time step is determined from Eq. (11), any smaller time step guarantees a stable and converged simulation.

(11) 21

International Journal of Impact Engineering 112 (2018) 15–29

X. Chen, A.H.C. Chan

Table 3 Crack propagation of dominant flexural cracks in inside glass with different mesh.

Initiation time (μs) Completion time (μs) Crack length (m) Propagation speed (m/s)

Mesh A

Mesh A

Mesh B

20 33 0.011 846.2

21 35 0.011 785.7

Table 4 Material properties used in the simulation.

Young's modulus (GPa) Poisson's ratio Fracture energy (J/m2) Tensile strength (MPa) Density (kg/m3)

Glass

PVB

Steel

70 0.25 3.9 30.0 2500

0.1 0.4 – – 1100

200 0.29 – – 7800

Mesh B fracture of outside glass ply was assumed in the FEM simulation, possible fractures in both outside and inside glass plies are considered in the FEM-DEM analysis to make the simulation more realistic. Material properties used in the FEM-DEM simulation are tabulated in Table 4. Contact penalty is Ep = 1.0 × 1012 Pa. In [14], a non-reflective boundary was used for the laminated glass along the direction that perpendicular to the impact. Since there is currently no such boundary in the FEM-DEM, the beam is modelled with enough length so that crack-forming process in the middle will not be influenced by the stress waves that are reflected from the two ends of the laminated glass. Referring to the classic stress wave theory from Kolsky [51], the velocity of dilatational wave propagation in a solid is

Initiation

c=

Completion

λ + 2μ ρ

(12) Eν (1 + ν )(1 − 2ν )

E

where λ and μ are Lame constants with λ = and μ = 2(1 + ν) . Obviously, the velocity of stress wave propagation in glass is much faster than that in the PVB interlayer. From Eq. (12), velocity of dilatational wave in glass is c = 5796.55 m/s. In the FEM-DEM modelling, the length of the laminated glass is set as L = 120 mm so that fracture behaviour within the impact effective area will not be influenced by reflected stress waves for 15–20 μs. Based on the above calculation, a 120 mm long 12.7 mm high laminated glass beam clamped in two rigid channels is modelled by the FEM-DEM, and ho, hi and hinter are same as those in Ref. [14]. Configurations of the FEM-DEM model is shown in Fig. 17. The FEM-DEM model is of 42,225 elements in total. 36,582 elements are meshed within the glass while 5292 elements are meshed within the PVB interlayer. Elements with a characteristic size of 0.2 mm are employed in the middle impact effective area while a coarse mesh is defined in the far filed. Time step is 0.01 ns. Fig. 18 shows the cracks of the laminated glass from both the FEMDEM and FEM at t = 7.6 μs. In [14], only the computational algorithm for cone damage was considered, and hence Flocker and Dharani obtained an ideal symmetric cone crack with the propagation angle of 30°. However, such ideal cone crack is not possible in the FEM-DEM simulation since: (1) mesh is not symmetric; and (2) other crack types can also be captured by the FEM-DEM. Referring to Fig. 18(a), a recognisable cone crack is obtained in the FEM-DEM simulation, where the propagation angles are 26.5° and 29° at both sides. The obtained propagation angles in FEM-DEM agree very well with that of FEM, and the majority damage occurs within the cone area. No damage in the inside glass ply is identified from the FEM-DEM simulation, demonstrating that the assumption “rupture only occurs in the outside glass ply” in the FEM simulation is correct. For dynamic impact, Knight et al. [52] pointed out that the cone angle varies with the impact velocity, and proposed a simple formulation for evaluating the cone angle α if initial impact velocity v0 is lower than 50 m/s

Fig. 15. Initial and completion crack patterns of the dominant flexural crack in inside glass.

4.3. Cone fracture Flocker and Dharani [14] presented analysis on fracture of laminated glass where explicit finite element code DYNA2D was employed. They proposed a method which allows traditional wave propagation approach to model small hard projectile impact on architectural laminated glass. In [14], a three-layer laminated window with two soda-lime glass plies and one PVB interlayer was considered. Fig. 16 illustrates the laminated glass unit for the analysis. The thicknesses for the outside and inside glass plies are ho = 4.76 mm and hi = 6.35 mm, respectively. Thickness of the PVB interlayer is 1.59 mm. The laminated glass subjects to a normal impact of a 3.97 mm radius steel ball at a velocity of v0 = 35.8 m/s. Flocker and Dharani assumed that the glass and interlayer are perfectly bonded. In addition, steel ball and the inside glass ply are modelled with linear elastic materials, and the PVB interlayer is modelled as a linear elastic solid. These assumptions suggest the fracture of this laminated window glass can be attributed to the failure of the outside glass ply only. Since no specific PVB properties are given in [14], standard PVB material properties are used in the FEM-DEM simulation from existing literature, e.g. Vallabhan [50]. Though only the

α=

v0 (90 − αst ) + αst 290

(13)

where αst is the cone angle for static indentation. According to Cook and Pharr [53], Hertzian indentation produces a cone with αst ≈ 22°. When v0 = 35.8 m/s, the cone angle should be α = 30. 4 °. With the consideration of relatively random crack distributions and uncertainty in

Fig. 16. Configurations of the lamented glass unit (after Flocker and Dharani [14]).

22

International Journal of Impact Engineering 112 (2018) 15–29

X. Chen, A.H.C. Chan

Fig. 17. Geometry and mesh configuration of the laminated glass in FEM-DEM simulation: (a) overall view; (b) mesh in the impact effective area.

FEM-DEM modelling, cone angles of 26.5°/29° are quite accurate and acceptable. Further, according to the observation of Chaudhri and Kurkjian [54], fine debris can be found ejecting around the contact point during impact, which cannot be observed in the FEM simulation but are available in the FEM-DEM result as small flying fragments are achieved in Fig. 18(a). This is an intrinsic merit of the FEM-DEM over the FEM on simulating the movement of fragments. In addition, traces of median and lateral cracks, which occur very soon after the start of impact, are identified by the FEM-DEM, suggesting genuine dynamic impact process is simulated. Table 5 summarises the major features from the results simulated by the FEM-DEM and FEM. It can be concluded that result from the FEMDEM reaches a good agreement with that from the FEM on cone angle, and its damage pattern is closer to the reality. Median and lateral cracks are found in the FEM-DEM simulation with fine debris around the impact point. In both the FEM-DEM and FEM simulations, the inside glass ply is intact until the end of the simulation. It is argued that the laminated glass may experience further deformation and flexural cracks in the inside glass ply may be observed with the increase of simulation time.

Table 5 Summary on the comparison of results between FEM-DEM and FEM.

Cone angle Median & lateral crack Fragments Inside glass damage

FEM-DEM

FEM

26.5°/29° Yes Yes No

30° No No No

4.4. Post-damage fragmentation Zang et al. [20] presented a study on the fracture and fragmentation behaviour of automobile laminated glass subjects to rigid impact by using spherical DEM. They investigated a clamped laminated glass with two glass plies and one PVB interlayer. In this section, 2D FEM-DEM analysis is performed to examine the post-damage fragmentation of the laminated glass, i.e. failure mode, shape and movement of fragments as well as relative positions of the projectile. The beam is 100 mm × 4.76 mm, with 0.76 mm thick for PVB and 2 mm thick for each glass ply. The laminated glass subjects to the impact of a 12 mm diameter ball with the mass of 1 kg at the velocity of

Fig. 18. Crack patterns of the laminated glass from FEM-DEM and FEM at t = 7.6 μs: (a) FEM-DEM; (b) FEM (after [14]).

23

International Journal of Impact Engineering 112 (2018) 15–29

X. Chen, A.H.C. Chan

Fig. 19. The configuration of the laminated glass in FEM-DEM simulation.

laminated glass reach a very good agreement. Conclusions can be drawn that in both the FEM-DEM and DEM modelling, the residual velocity fields of projectiles are highly similar, suggesting in the postdamage stage the movement of the projectile in FEM-DEM agrees very well with that from traditional DEM. Furthermore, the final longitudinal damaged area in both FEM-DEM and DEM simulations are almost identical. Despite these similarities, differences between the fracture and fragmentation behaviours of the FEM-DEM and traditional DEM are obvious and can be observed as follows:

Table 6 Material and interfacial properties for the FEM-DEM simulation.

3

Density (kg/m ) Young's modulus (GPa) Poisson's ratio Tensile strength (MPa) Shear strength (MPa) Fracture energy (J/m2)

Glass

PVB

Projectile

2500 74.1 0.2 34.6 17.9 4.0

101 0.006 0.42 18.62 17.9 20.0

1.105 × 106 200 0.3 – – –

30 m/s. In the FEM-DEM simulation, the laminated glass is modelled with 7202 elements (6148 elements for glass and 1054 elements for interlayer) with a characteristic size of 0.38 mm (Fig. 19). Time step is 0.1 ns. Material and interfacial properties used in this analysis are given in Table 6. The contact penalty Ep is 1.0 × 1012 Pa. To enable possible fracture of interlayer, a fracture energy of 20.0 J/m2 is assumed for the PVB interlayer. Potential damage for the glass-interlayer interface is also considered and a fracture energy of 100 J/m2 is used. Fig. 20 shows the entire fracture and fragmentation process (t = 0 to 0.5 ms) from both the FEM-DEM and DEM simulations. It can be concluded that the transient responses in Fig. 20(a) and (b) are quantitatively similar. In both simulations, the projectile penetrates the laminated glass, and relative positions between the projectile and

i) Damages of interlayer are different. Before t = 0.1 ms, the PVB interlayer in the FEM-DEM simulation has ruptured, while it experiences quite large deformation in the DEM simulation and was almost intact before t = 0.2 ms. Since ball connections are used in the DEM modelling, the deformation of interlayer is much larger than that in the FEM-DEM, where triangular elements are relatively rigid and cannot undergo such significant deformation. ii) Types of fragments are different. In the FEM-DEM, besides numerous elemental-sized fragments, there were several large pieces of laminated glass fragments (with both glass and PVB elements). However, in the DEM simulation most fragments are linked-particles. iii) In general, it can be claimed that the results from the FEM-DEM is Fig. 20. Impact fragmentation sequences from top to bottom at t = 0, 0.1, 0.2, 0.3, 0.4 and 0.5 ms: (a) FEM-DEM; (b) DEM (after Zang et al. [20]).

24

International Journal of Impact Engineering 112 (2018) 15–29

X. Chen, A.H.C. Chan

Fig. 21. Configurations of the laminated glass beam: (a) geometry (after Xu and Zang [22]); (b) FEM-DEM mesh.

projectile in the outside ply and propagate through its thickness until t = 200 μs (sequence number 3). A dominant flexural crack which can be observed at t = 300 μs (sequence number 4) is formed in the inside glass ply and propagates towards the interlayer. Some small fragments are identified near the contact surface in the FEM-DEM simulation. Further investigation shows that the formation time for the dominant flexural crack (0.011 m length) in the inside glass is 1.2 × 10-5s, and the propagation speed is approximately 916.7 m/s. Fig. 24 shows the evolution of energy along with time. At the beginning, the initial input energy, total kinetic energy and the kinetic energy of projectile are the same. With the time elapsing, both total kinetic energy and the kinetic energy of projectile drop. Ideally, zone I between curves of initial input energy and total kinetic energy equals to the energy consumed for fracture and elastic strain; while zone II between curves of the total kinetic energy and kinetic energy of projectile relates to the oscillation of laminated glass and free movement of fragments. Based on the above comparisons, conclusion can be drawn that results from the FEM-DEM simulation agree very well with the data from experiments. With the increase of recording time interval, fracture in the inside glass ply (e.g. the dominant flexural crack) is observed at a similar location and time as is captured in the experiment.

closer to the reality and a better structural simulation is achieved.

4.5. Fracture under impact of a rectangular projectile Xu and Zang [22] conducted a number of impact tests on laminated glass beams with high speed camera images recorded. The laminated glass beam which is constrained by four fixed supporters (Fig. 21(a)), is 200 mm × 24 mm × 10 mm. ho = hi = 10 mm and hinter = 4 mm. The projectile is 4 mm × 4 mm × 10 mm with a mass of 1 kg, and it hits the middle of the beam at a velocity of 3.13 m/s. An FEM-DEM simulation is conducted and the triangular mesh is shown in Fig. 21(b). There are 10,106 elements with a characteristic size of 1 mm for the laminated glass body. Time step is 0.5 ns. Ref. [22] suggests the interlayer is perfect elastic and interfaces between the glass and the interlayer are well bonded, so that rupture only occurs within the glass plies. Material properties for the analysis are tabulated in Table 7, and contact penalty Ep = 1.0 × 1012 Pa is used. The corresponding fracture history of the laminated glass from the FEM-DEM simulation is presented and compared with high speed camera records. In Fig. 22, time interval between two consecutive images is 20 μs. Only crack evolution in the outside glass ply is observed while the inside glass ply is kept intact up to t = 140 μs for sequence number 8. Though the simulated crack type in the outside glass ply is not exactly the same as that in the test records, the trends and occurrence time match well. The cracks become more and more obvious from t = 60 μs (sequence number 4). In general, results from the FEM-DEM and experiment agree well for this 20 μs-interval comparison. Fig. 23 shows the failure sequences from the FEM-DEM modelling and another experiment with 100 μs interval between each image. Referring to Fig. 23(a), cracks start to develop just beneath the

4.6. 3D fracture In this sub-section, impact fracture of a 3D laminated glass example inspired by Timmel et al. [55] is presented. In [55], a 1500 mm × 1500 mm square laminated glass plate is clamped along its four edges by a rigid frame. The total thickness of the plate is 5.38 mm (5 mm for glass ply and 0.38 mm for the interlayer). A rigid spherical projectile with 300 mm in diameter and 70 kg in mass hits the centre of the plate at a velocity of 10 m/s. In [55], glass plies are modelled with coinciding shell elements and the PVB interlayer is modelled as a hyperelastic membrane. Additionally, no rupture in interlayer is considered. Since only one shell layer is employed for the two glass plies, it is held in [55] that both glass plies fail at the same time. Some reasonable simplifications have to be made in the FEM-DEM simulation for the sake of computational efficiency. The thickness of interlayer is set to be 5 mm. It is well attached to and beneath the 5 mm glass sheet, thus simultaneous fracture of both outside and inside glass

Table 7 Material properties for the impact fracture analysis of laminated glass beam.

Density (kg/m3) Young's modulus (GPa) Poisson's ratio Tensile strength (MPa) Fracture energy (J/m2)

Glass

PVB

Projectile

Support

2500.0 74.09 0.2 100.0 10.0

100.0 0.05 0.48 –– –

– 210.0 0.269 – –

2400.0 5.0 0.4 – –

25

International Journal of Impact Engineering 112 (2018) 15–29

X. Chen, A.H.C. Chan

Fig. 22. Comparison between the FEM-DEM simulation and the experiment with a 20 μs interval: (a) FEM-DEM results; (b) high speed camera records.

5

plies is achieved. Configuration of the laminated glass plate in FEMDEM simulation is schematically shown in Fig. 25. Material properties of glass and PVB used for the FEM-DEM analysis are same as those in Table 4, while Poisson's ratio of glass is 0.23 in accordance with [55]. There are 1818 elements and the time step is 1.0 ns. The contact penalty is 1.0 × 1012 Pa. Fracture pattern at t = 1.6 ms simulated by the FEM-DEM is given in Fig. 26. It can be observed that cracks occur along both radial and circumferential directions, and a typical spider-shaped crack is obtained. This spider crack qualitatively validates the effectiveness of the FEM-DEM approach, as similar cracking patterns are commonly found in 3D laminated glass impact literature, e.g. [5,17]. Comparing with results from previous 2D simulations, 3D modelling is capable of representing the impact cracking behaviour of laminated glass from a perspective view. Spider cracks, which are not available in 2D due to the natural limitations, are observed in 3D simulation. Albeit a coarse mesh has to be employed due to low computational efficiency in 3D [40], FEM-DEM simulation on the impact fracture of laminated glass plate is still successful and beneficial. And on the other hand, fine mesh can be employed in 2D modelling, and detailed cracking behaviour along the thickness can be examined.

I

Energy (J)

4.5 II 4 Initial input energy Total kinetic energy

3.5

Projectile kinetic energy

3 0

0.1

0.2

0.3

0.4

0.5

Time (ms) Fig. 24. Evolution of energy along with time.

4.7. Summary Numerical examples on fracture and fragmentation simulated by the FEM-DEM program are presented and validated with existing computational and experimental results. In Section 4.1, very low impact velocities are applied to the projectile and the elastic impact behaviour of laminated glass is investigated and verified with that from ABAQUS; in Section 4.2, convergence study on the residual velocity of projectile, total kinetic energy, load-deflection and crack propagation speed against change of mesh is performed; in Section 4.3, a cone type crack

Fig. 25. Configuration of the laminated glass plate in the FEM-DEM simulation.

Fig. 23. Comparison between the FEM-DEM simulation and the experiment with a 100 μs interval: (a) FEM-DEM results; (b) high speed camera records.

26

International Journal of Impact Engineering 112 (2018) 15–29

X. Chen, A.H.C. Chan

Fig. 26. Crack pattern of the laminated glass plate by the FEM-DEM simulation.

three-dimensional impact fracture example is presented and a radialcircumferential crack is obtained. Crack propagation velocity in glass is estimated around 1000 m/s from the FEM-DEM simulations. It was pointed out that the maximum crack speed cannot exceed the Rayleigh wave speed [7,56–58]. The obtained values of crack propagation speed in this paper are well less than the corresponding Rayleigh surface wave speed in glass, which is about 3100 m/s given in [14]. Since these examples are quite representative and verified with results from widely acknowledged methods, it can be concluded that the FEM-DEM together with the fracture models in Section 3 are capable of yielding good results on the fracture and fragmentation behaviour of laminated glass under hard body impact. It is worth mentioning that the employed fracture model is rateindependent, which may not be seamlessly applicable to high velocity impact where high strain rate is involved. However, rate-independent fracture models are acceptable for low velocity impact on glass. Holmquist and Johnson [59] gave the relation between strength f and the strain rate ɛ˙ in glass as

f = f0 (1 + C ln ɛ˙)

where f0 is the strength at ɛ˙ = 1.0 and C is a dimensionless constant (C = 0.003 according to [60]). Impact velocities investigated in this paper correspond to strain rates much smaller or around the level of 102 s−1 in glass ply according to Behr et al. [61], showing that glass is not that sensitive to strain rate. In practice, rate-independent approaches also have been successfully employed in existing impact fracture considerations of laminated glass [17,22]. For higher impact velocities at the magnitude of 103m/s, more sophisticated fracture models e.g. JHB model [59] with the consideration of high strain rate (> 104 s−1), high pressure effect and friction-heat transfer are necessary for the better representation of the physics.

Fig. 27. Failure responses with different Young's modulus for the interlayer.

0%

Normalised displacement

(14)

5. Further simulations for different Young's modulus of interlayer

-10%

The interlayer in the examined examples in Section 4 are all PVBbased. PVB interlayers are of good elasticity and capable of undergoing large deformation. In recent years, interlayers made from a stiff resin, e.g. SentryGlas®Plus (SGP), are employed in structural laminated glass applications [62–63]. One of the important differences between PVB and SGP is the value of the Young's modulus Einter. For PVB, Einter is around 100 MPa while it is about 1 GPa for SGP. According to Hooke's law, for a given stress σ, a higher modulus E requires a smaller strain ɛ to achieve equilibrium. In other words, the interlayer will be stiffer with the increase of its Young's modulus. In this section, further simulations are conducted based on the 1 m long 22 mm high laminated glass beam examined in Section 4.1. To enable the potential interfacial debonding between glass plies and interlayers, a 200 J/m2 fracture energy is defined for the interface joint. All other material properties are same as those in Section 4.1. Influences of Einter on impact failure responses of laminated glass are investigated using a range of values from 0.1 to 2.0 GPa, which is reasonable for commonly used interlayer materials. Selected failure modes of the laminated glass beam under v = 18 m/ s impact at t = 1 ms with different Einter are shown in Fig. 27. The

0.1GPa

-20%

0.3GPa 0.6GPa

-30%

1.2GPa 2.0GPa

-40% 0

0.2

0.4

0.6

0.8

1

Time (ms) Fig. 28. Normalised vertical deformation (Δ/h) against time for different Einter.

in a laminated glass beam is obtained and compared with the FEM results; in Section 4.4, impact fragmentation of laminated glass is examined, providing good similarity with results from DEM; in Section 4.5, impact response under a rectangular projectile is studied by the FEM-DEM, fracture sequences are compared with high speed camera images and very good agreement is reached; in Section 4.6 a 27

International Journal of Impact Engineering 112 (2018) 15–29

X. Chen, A.H.C. Chan

Nomalised total kinetic energy

100%

failure of other brittle or quasi-brittle multi-layered materials should it is further developed.

80%

Acknowledgements The authors appreciate the financial support from China Scholarship Council (No. 2009621028), the University of Birmingham (Postgraduate Teaching Assistantship) and Jiangsu Shuang–Chuang Program for this research. Supports on the FEM-DEM from Professor Antonio Munjiza are gratefully acknowledged.

60% 40% 0.1GPa 20%

1.0GPa

References

0% 0

0.2

0.4 0.6 Time (ms)

0.8

[1] Behr RA, Minor JE, Norville HS. Structural behaviour of architectural laminated glass. J Struct Eng 1993;119(1):202–22. [2] Ledbetter LR, Walker AR, Keiller AP. Structural use of glass. J Archit Eng 2006;12(3):137–49. [3] Hooper JA. On the bending of architectural laminated glass. Int J Mech Sci 1973;15:309–23. [4] Behr RA, Minor JE, Linden MP, Vallabhan CVG. Laminated glass units under uniform lateral pressure. J Struct Eng ASCE 1985;111(5):1037–50. [5] Zhang X, 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] Marcon B, Fouvry S, Sassy O, Guegan J, Daniel G. Fracture mechanics of impacted laminated glass subjected to various fatigue stressing conditions. Eng Fract Mech 2014;127:71–82. [7] Rabczuk T, Belytschko T. Cracking particles: a simplified meshfree method for arbitrary evolving cracks. Int J Numer Methods Eng 2004;61(13):2316–43. [8] Rabczuk T, Belytschko T. A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Comput Methods Appl Mech Eng 2007;196(29–30):2777–99. [9] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A simple and robust three-dimensional cracking-particle method without enrichment. Comput Methods Appl Mech Eng 2010;199(37–40):2437–55. [10] Diyaroglu C, Oterkus E, Madenci E, Rabczuk T, Siddiq A. Peridynamic modeling of composite laminates under explosive loading. Compos Struct 2016;144:14–23. [11] Ren H, Zhuang X, Cai Y, Rabczuk T. Dual-horizon peridynamics. Int J Numer Methods Eng 2016;108(12):1451–76. [12] Ren H, Zhuang X, Rabczuk T. Dual-horizon peridynamics: a stable solution to varying horizons. Comput Methods Appl Mech Eng 2017;318:762–82. [13] Flocker FW, Dharani LR. Stress in laminated glass subject to low velocity impact. Eng Struct 1997;19(10):851–6. [14] Flocker FW, Dharani LR. Modelling fracture in laminated architectural glass subject to low velocity impact. J Mater Sci 1997;32:2587–94. [15] Du Bois PA, Kolling S, Fassnacht W. Modelling of safety glass for crash simulation. Comput Mater Sci 2003;28:675–83. [16] Pyttel T, Liebertz H, Cai J. Failure criterion for laminated glass under impact loading and its application in finite element simulation. Int J Impact Eng 2011;38(4):252–63. [17] Chen S, Zang M, Xu W. A three-dimensional computational framework for impact fracture analysis of automotive laminated glass. Comput Methods Appl Mech Eng 2015;294:72–99. [18] Oda J, Zang MY, Mori T, Tohyama K. Simulation of dynamic fracture behavior of laminated glass by DEM. Trans 8th calculation dynamics symp JSME. 1995. p. 429–30. [19] Oda J, Zang MY. Analysis of impact fracture behaviour of laminated glass of bi-layer type using discrete element method. Key Eng Mater 1998;145–149:349–54. [20] 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. [21] Gao W, Zang M. The simulation of laminated glass beam impact problem by developing fracture model of spherical DEM. Eng Anal Boundary Elem 2014;42:2–7. [22] Xu W, Zang M. Four-point combined DE/FE algorithm for brittle fracture analysis of laminated glass. Int J Solids Struct 2014;51:1890–900. [23] Munjiza A. The combined finite-discrete element method. John Wiley and Sons; 2004. [24] Munjiza A. Discrete elements in transient dynamics of fractured media Ph.D. thesis Department of civil engineering, University of WalesSwansea; 1992. [25] Munjiza A, Owen DRJ, Bicanic N. A combined finite-discrete element method in transient dynamics of fracturing solids. Eng Comput 1995;12:145–74. [26] Munjiza A, Andrews KRF. NBS contact detection algorithm for bodies of similar size. Int J Numer Methods Eng 1998;43(1):131–49. [27] Munjiza A, Andrews KRF, White JK. Combined single and smeared crack model in combined finite-discrete element analysis. Int J Numer Methods Eng 1999;44(1):41–57. [28] Munjiza, A., 2014. Manual for the “Y” FEM/DEM computer program. [29] Munjiza A, John NWM. Mesh size sensitivity of the combined FEM/DEM fracture and fragmentation algorithms. Eng Fract Mech 2002;69:281–95. [30] Munjiza A, Lei Z, Divic V, Peros B. Fracture and fragmentation of thin shells using the combined finite-discrete element method. Int J Numer Methods Eng 2013;95:478–98. [31] Lukas T, Schiava D'Albano GG, Munjiza A. Space decomposition based parallelization solutions for the combined finite-discrete element method in 2D. J Rock

1

Fig. 29. Normalised total kinetic energy (Tcurrent/T) with different Einter.

velocity is such chosen that both outside and inside glass plies close to the contact point are ruptured, and deformation of the beam is largely controlled by the interlayer. It can be observed from Fig. 27 that when Einter is small (e.g. Einter = 0.1 GPa), the majority fractured glass in both outside and inside plies adheres to the interlayer, with obvious bending of the interlayer. With the increase of Einter from 0.5 to 2.0 GPa, debonding of the glass-interlayer interface becomes noticeable, and in addition, less bending of the interlayer is obtained. The trend can be concluded that the larger the Einter, the more severe the debonding is. Fig. 28 shows the normalised vertical deformation Δ/h (Δ is the actual vertical deformation and h is the total thickness of the laminated glass) for point B (definition of point B can be found in Section 4.1) against time with different Einter. It can be concluded from Fig. 28 that a smaller Einter (e.g. Einter = 0.1 GPa) yields a larger deformation, indicating that the laminated glass is ‘softer’. With the increase of Einter, Δ/ h decreases and gets a much smaller value than that for Einter = 0.1 GPa. Fig. 29 presents the normalised total kinetic energy of the examined laminated glass system with different Einter under the 18 m/s impact. Notation Tcurrent and T are the current and initial total kinetic energies of the system, respectively. It can be found that laminated glass with different Einter exhibit similar energy reduction capacities, and the difference is less than 5%. This suggests that a stiff interlayer has a similar energy absorption capacity as a soft interlayer, where the small difference can be ignored. Based on these investigations, it is worthwhile to use a stiff resin like SGP in engineering applications if displacement and stiffness control are vital.

6. Conclusions The FEM-DEM is employed in this paper for the fracture and fragmentation analysis of laminated glass under hard body impact. Cohesive fracture models for the failure of glass, interlayer and glassinterlayer interface are discussed respectively. Numerical examples on the impact facture and fragmentation of laminated glass are simulated and validated. The FEM-DEM results are found to be in good agreement with data in existing literature, showing that the present method produces highly realistic and reliable failure patterns. And advantages of FEM-DEM modelling on the impact fracture and fragmentation behaviour of laminated glass over its parent FEM and DEM are demonstrated. Further discussions on the Young's modulus of interlayer Einter reveals that use of a stiff resin like SGP can effectively improve the deformation performance at little expense of energy absorption, providing guidance to laminated glass manufacturing and engineering applications. Though only laminated glass is investigated in this work, the FEMDEM and the employed fracture models can be used to simulate the 28

International Journal of Impact Engineering 112 (2018) 15–29

X. Chen, A.H.C. Chan

p. 513–20. [47] Sedlacek G, Blank K, Gusgen J. Glass in structural engineering. Struct Eng 1995;73(2):17–22. [48] Overend M, Zammit K. A computer algorithm for determining the tensile strength of float glass. Eng Struct 2012;45:68–77. [49] Ballarini R, Pisano G, Royer-Carfagni G. The lower bound for glass strength and its interpretation with generalized Weibull statistics for structural applications. J Eng Mech 2016;142(12):04016100. [50] Vallabhan CVG, Das Y, Ramasamudra M. Properties of PVB interlayer used in laminated glass. J Mater Civil Eng 1992;4(1):71–6. [51] Kolsky H. Stress waves in solids. Oxford: Clarendon Press; 1953. [52] Knight CG, Swain MV, Chaudhri MM. Impact of small steel spheres on glass surfaces. J Mater Sci 1977;12(8):1573–86. [53] Cook RF, Pharr GM. Direct observation and analysis of indentation cracking in glasses and ceramics. J Am Ceram Soc 1990;73(4):787–817. [54] Chaudhri MM, Kurkjian CR. Impact of small steel spheres on the surfaces of “normal” and “anomalous” glasses. J Am Ceram Soc 1986;69(5):404–10. [55] Timmel M, Kolling S, Osterrieder P, Du-Bois PA. A finite element model for impact simulation with laminated glass. Int J Impact Eng 2007;34:1465–78. [56] Freund LB. Crack propagation in an elastic solid subjected to general loading—I. constant rate of extension. J Mech Phys Solids 1972;20(3):129–40. [57] Freund LB. Crack propagation in an elastic solid subjected to general loading—II. non-uniform rate of extension. J Mech Phys Solids 1972;20(3):141–52. [58] Larcher M. Development of discrete cracks in concrete loaded by shock waves. Int J Impact Eng 2009;36(5):700–10. [59] Holmquist TJ, Johnson GR. A computational constitutive model for glass subjected to large strains, high strain rates and high pressures. J Appl Mech-ASME 2011;78(5):051003. [60] Holmquist TJ, Johnson GR, Lopatin CM, Grady DE, Hertel ESJ. High strain rate properties and constitutive modeling of glass. 15th international symposium on ballistics. 1995. p. T14929. [61] Behr RA, Kremer PA, Dharani LR, Ji FS, Kaiser ND. Dynamic strains in architectural laminated glass subjected to low velocity impacts from small projectiles. J Mater Sci 1999;34(23):5749–56. [62] Bennison SJ, Smith CA, Van Duser A, Jagota A. Structural performance of laminated glass made with a “stiff” interlayer. glass in buildings. In: Block V, editor. ASTM international: West Conshohocken, PA; 2002. p. 57–64. ASTM STP 1434. [63] Delincé D, Callewaert D, Belis J, Van Impe R. Post-breakage behaviour of laminated glass in structural applications. In: Veer BL, editor. Challenging glass conference on architectural and structural applications of glass. 2008. p. 459–66.

Mech Geotech Eng 2014;6:607–15. [32] Munjiza A, Rougier E, Knight EE. Large strain finite element method: a practical course. John Wiley and Sons; 2015. [33] Mahabadi OK, Lisjak A, Munjiza A, Grasselli G. Y-Geo: new combined finite-discrete element numerical code for geomechanical applications. Int J Geomech 2012;12(6):676–88. [34] Liu HY, Kang YM, Lin P. Hybrid finite-discrete element modeling of geomaterials fracture and fragment muck-piling. Int J Geotech Eng 2013;9(2). 1939787913Y.000. [35] Živaljić N, Nikolić Ž, Smoljanović H. Computational aspects of the combined finite–discrete element method in modelling of plane reinforced concrete structures. Eng Fract Mech 2014;131:669–86. [36] Guo L, Latham J-P, Xiang J. Numerical simulation of breakages of concrete armour units using a three-dimensional fracture model in the context of the combined finite-discrete element method. Comput Struct 2015;146:117–42. [37] Smoljanović H, Živaljić N, Nikolić Ž. A combined finite-discrete element analysis of dry stone masonry structures. Eng Struct 2013;52:89–100. [38] Smoljanović H, Nikolić Ž, Živaljić N. A combined finite–discrete numerical model for analysis of masonry structures. Eng Fract Mech 2015;136:1–14. [39] Smoljanović H, Nikolić Ž, Živaljić N. A finite-discrete element model for dry stone masonry structures strengthened with steel clamps and bolts. Eng Struct 2015;90:117–29. [40] Chen X, Chan AHC, Yang J. Simulating the breakage of glass under hard body impact using the combined finite-discrete element method. Comput Struct 2016;177:56–68. [41] Chen X. Investigation on the impact damage of glass using the combined finite/ discrete element method PhD thesis, UK: School of Civil Engineering, University of Birmingham; 2013 [42] Chen X, Chan AHC, Yang J. FEM/DEM modelling of hard body impact on the laminated glass. Appl Mech Mater 2014;553:786–91. [43] Anderson TL. Fracture mechanics: fundamentals and applications. 2nd ed Boca Raton: CRC Press; 1991. [44] Zhang X, 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. [45] Tvergaard V, Hutchinson JW. The relation between crack growth resistance and fracture process parameters in elastic–plastic solids. J Mech Phys Solids 1992;40(6):1377–97. [46] Pelfrene J, Van Dam S, Van Paepegem W. Numerical and experimental study of the peel test for assessment of the glass-PVB interface properties in laminated glass. In: Belis J, Louter C, Bos F, Lebet J-P, editors. Proceedings of challenging glass 4 & COST action TU0905 final conference. London, UK: Taylor & Francis Group; 2014.

29