Simulating the breakage of glass under hard body impact using the combined finite-discrete element method

Computers and Structures 177 (2016) 56–68

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Simulating the breakage of glass under hard body impact using the combined finite-discrete element method Xudong Chen a,1, Andrew H.C. Chan b,⇑,1, Jian Yang c,d a

School of Naval Architecture and Civil Engineering, Jiangsu University of Science and Technology, Zhangjiagang 215600, China School of Engineering and ICT, University of Tasmania, Hobart 7001, Australia c School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao Tong University, Shanghai 200240, China d School of Civil Engineering, University of Birmingham, Birmingham B15 2TT, UK b

a r t i c l e

i n f o

Article history: Received 21 November 2015 Accepted 31 August 2016

Keywords: Glass Breakage Fracture Impact Combined FEM-DEM

a b s t r a c t The fracture behaviour of glass under hard body impact is simulated using the combined finite-discrete element method with a finite element incorporated into the discrete elements. A Mode I fracture model is employed and proven to be successful in modelling the crack initiation and tracking the crack propagation in glass. Cone, flexural and punching damages are obtained thus qualitatively verified the glass fracture model. A parametric study on tensile strength, fracture energy and impact velocity further examines the breakage behaviour of glass, demonstrating the applicability of this computational technique for the analysis of glass impact mechanism. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction The time of glass making by human can be traced back to around 10,000 BC in Egypt [1]. Around 1960, Sir Alastair Pilkington introduced a revolutionary manufacturing method for float glass, which is still in use today and covers 90% float glass production. However, due to the brittle nature of glass, cracking may occur under impact and often results in a catastrophic collapse, bringing risks to the users. Consequently, the fracture and breakage of glass under impact is a topic of continuing research. Preston [2] demonstrated that the manner of breakage in annealed glass is directly related to the appearance of fracture surfaces and this view has been supported by Murgatroyd [3] and Shand [4]. Some of the general fracture features in glass were described by many other researchers and can be found in Refs. [5–7]. In building construction, glass is commonly used as balconies, facades and other structural members. The stress and safety for ordinary glass liable to human impact were discussed by Toakley [8] in order to minimise possible human injuries. Minor [9] examined the responses of building glazing subjects to windborne debris impact, with the consideration of both small and large missiles. ⇑ Corresponding author. E-mail addresses: [email protected] (X. Chen), [email protected] (A.H.C. Chan), [email protected] (J. Yang). 1 Formerly School of Civil Engineering, University of Birmingham, Birmingham B15 2TT, UK. http://dx.doi.org/10.1016/j.compstruc.2016.08.010 0045-7949/Ó 2016 Elsevier Ltd. All rights reserved.

The finite element method (FEM) and the discrete element method (DEM) have been used for the numerical simulation of the breakage of glass. Based on different theories, appropriate fracture models were developed and incorporated into the FEM or DEM to investigate the fracture behaviour of glass, e.g., the continuum damage mechanics (CDM) model [10], the fracture mechanics approach [11] and the two-parameter Weibull distribution [12]. Recently, Pyttel et al. [13] presented a fracture criterion for laminated glass and implemented it into an explicit finite element solver. They suggested that the crack initiation is based on the critical energy threshold while propagation is related to ‘‘local Rankine” (maximum stress) behaviour. In general, the FEM has been employed in crack initiation prediction with reasonable success, however, it is difficult for the FEM to model the post-damage behaviour, which is a major concern in glass design. Compare with the FEM, the DEM is less frequently used in glass fracture analysis. Oda et al. [14] employed the DEM to model the impact behaviour of laminated glass. Later, the method was extended to double layer glass [15]. Impact behaviours of both single and laminated glass subjects to ball impact were investigated by Zang et al. [16] using spherical discrete elements. In [16], penetration as well as fragmentation was predicted; however, the fracture patterns obtained are not realistic. Their recent research [17,18] attempted to employ discrete elements for potential fracture regions while finite elements in the non-fracture regions, but realistic crack patterns are still not readily attained. In general,

57

X. Chen et al. / Computers and Structures 177 (2016) 56–68

the DEM is capable of simulating the post-damage behaviour, including the kinetic movement and interactions of fragments, but with low precision in predicting crack initiation and fracture patterns. In order to overcome the respective limitations of FEM and DEM in glass impact simulations, a recently introduced numerical method, i.e., the combined finite-discrete element method (FEMDEM) [19,20], is employed in this paper. In this method, structures are first discretized into a number of discrete elements which are connected to each other by interfaces that incorporate a fracture model. Each discrete element interacts with those that are in contact with it. Within a discrete element, finite element mesh can be formulated, resulting into a more accurate estimate of contact forces between elements and the deformation of the structure, thus enabling the cracks initiate and propagate more accurately and fragments move in a more rational manner. The FEM-DEM was pioneered and developed by Munjiza and co-workers during 1990s. A newly published book [21] provides ‘down to earth’ presentation of the framework for developing large-strain based nonlinear material laws and the large-strain and large-displacement based FEM-DEM. With several applications of the FEM-DEM on fracture consideration, more attention is given to the geologic behaviour of rock [22] or the damage of concrete. For the breakage of glass, less attention was given previously. A recent paper [23] presented the fracture and fragmentation of 3D glass thin shells under impact using shell FEM-DEM elements, but results were not compared with other numerical or experimental data and no parametric study was given. The present study emphasises on the application of the FEM-DEM to simulate glass breakage under hard body impact. For the results reported, a Mode I fracture model for glass modified from a combined single and smeared crack model [24] is discussed in details and implemented into the FEM-DEM computer program Y [19]. Both 2D and 3D glass impact simulations from the FEM-DEM are validated with numerical or experimental results, showing that the FEM-DEM is applicable, reliable and robust for the fracture mechanism analysis of glass under impact loading. A parametric study on tensile strength, fracture energy and impact velocity is performed and comparative benchmarks are provided for future research.

The combined FEM-DEM is a numerical method developed to analyse mechanics of solids while considering it as a combination of both continua and discontinua. In this method, the deformability and stress of each single discrete element is described by a standard continuum formulation (FEM) while the contact and motion of elements are considered by discontinua equations (DEM). For a single discrete element i, translational and rotational motions are controlled by the net external force Fi and torque Ti respectively according to Newton’s second law. 2

Ii

d

dt

r 2 i

¼ Fi

d xi ¼ T i dt

dA

Pt,Pc dF Fig. 1. Contact force in the FEM-DEM.

dF ¼ Ep ½graduc ðPc Þ  gradut ðPt ÞdA

ð3Þ

where Pc and Pt are the points share the same coordinate on S belonging to the contactor and the target, respectively; uc and ut are potentials for the contactor and the target, respectively; Ep is a contact penalty parameter, and grad represents for the gradient. By integrating over the area S, the 2D contact force is obtained in Eq. (4).

Z F ¼ Ep S

½graduc ðPc Þ  gradut ðPt ÞdA

ð4Þ

The 3D contact force can be simply extended by replacing the overlapping area S in Eq. (4) with an overlapping volume V. In the current FEM-DEM computer program, 3-node constant strain triangles and 4-node constant strain tetrahedrons are implemented. Fracture criterion is implemented between each pair of discrete elements in contact using linked interfaces (Fig. 2). For the 2D, line interfaces are employed. Similarly, triangular interfaces are used in the 3D. Interfaces can break according to the separation d between the two adjacent element boundaries. The change of d follows a suitable fracture law based on Mode I breakage. In the present model, cracks are assumed to occur along element edges [25], and possible mesh dependency can be reduced to the minimum if unstructured and small enough elements are used [26]. Details of the facture law employed in this work are given in the next section.

3. Fracture model of glass

2. The combined FEM-DEM

mi

S

ð1Þ

ð2Þ

In Eqs. (1) and (2), mi is the mass of discrete element i; ri is its position; Ii is the moment of inertia and xi is the angular velocity. Explicit numerical integration determines the velocity and position of a particular discrete element. In 2D, contact forces between elements are evaluated by the overlapping area S (Fig. 1). The penetration of any elemental area dA of the contactor into the target results in an infinitesimal contact force dF [19]

For the fracture of brittle materials, e.g., glass, Mode I tension breakage is most relevant [27]. There are also evidences [28] indicating that the dominant mechanism of brittle fracture in compression is still Mode I cracking. Consequently, a Mode I fracture model is adopted for glass. The cracking model used in this work is similar to Hillerborg’s [29]. The cohesive model originally proposed by Hillerborg is widely used in the fracture analysis of quasi-brittle materials, e.g., concrete or rock. From the view point of computational mechanics, quasi-brittle (concrete) and brittle (glass) are relative concepts relying on their extents of brittleness i.e., the magnitude of energy dissipation during softening. Theoretically, brittle material e.g., glass also has a plastic zone, albeit, considerably smaller. And the application of this model to glass, which is more brittle than concrete is computationally acceptable. When implementing the computational model, the plastic zone which is represented by a strain softening curve should be considered though it is very small. The advantages of using cohesive models with a strain softening curve are: (i) Brittle behaviour can be depicted should a small fracture energy is assigned. (ii) Computational stability can be obtained as a sudden drop of stress is avoided.

58

X. Chen et al. / Computers and Structures 177 (2016) 56–68

σ ft

B

A

(a)

(b)

O

Fig. 2. Interface joint elements in the FEM-DEM where dotted lines in the middle represent the initial positions: (a) 2D and (b) 3D.

Furthermore, Holloway [5] suggested ‘‘the fracture process in glass is not ideally brittle” from experimental measurements, which also provides theoretical grounds of using cohesive models. In Hillerborg’s model, a crack will develop when the bonding stress r reaches the tensile strength ft, which is the onset of the damage. Afterwards, the bonding stress is assumed to decrease gradually as the crack width d increases. When d reaches a critical value dc, the stress drops to zero and the crack is completely open. This process is depicted in a descending curve, with the area under the curve equals to the fracture energy as is given in Eq. (5).

Z

dc

G¼ 0

ð5Þ

rdd

In this work, standard continuum formulation (FEM) is used before the material reaches its tensile strength. Afterwards, the smeared crack approach is applied, capturing the damage process through a weak constitutive relation of the material, enabling the cracks being smeared out over the continuum. Since no preexisted crack is assumed, crack initiation is detected directly without any bias along the interface between discrete elements. A typical strain hardening and softening curve is given in Fig. 3. The curve can be divided into two parts: zone A is the elastic region and zone B represents strain softening. When the bonding stress reaches ft, the crack initiates with a crack tip process zone where stress declines. As the bonding stress r drops to zero, a maximum separation will occur. At this point, two adjacent discrete elements in contact are considered completely dissociated. Separation d between adjacent element edges is defined in [19]. It is noted that dP is the onset separation, which equals to the separation when the bonded stress reaches the tensile strength ft. dc is the ultimate separation when the bonded stress decreases to zero. The calculation of bonded tensile stress r is classified in condition of: (i) d < 0; (ii) 0 6 d 6 dp ; and (iii) d > dp. The complete relations between the bonding stress r and the separation d are given in Eq. (6) from Ref. [19].

rn ¼

8  2  > d d > ft  2 > d dp p < f tz > > > :2 d f dp t

0 6 d 6 dp dp < d 6 dc

ð6Þ

d<0

where z is a heuristic parameter.

  a þ b  1 Dðaþbc=ððaþbÞð1abÞÞÞ ½að1  DÞ þ bð1  DÞc  e z¼ 1 aþb

ð7Þ

In Eq. (7), a, b and c are constant parameters; D is an independent variable indicating the fracture damage index and D 2 [0, 1]. For dp < d 6 dc , D ¼ ðd  dp Þ=ðdc  dp Þ; for d > dc, D = 1, which suggests that no bonding exists between two adjacent elements, leading to the total fracture and free movement. It can be mathematically argued that Eq. (7) is C1 continuous on [0, 1] with respect to D,

δp

δc

δ

Fig. 3. A typical strain hardening and softening curve.

which avoids any possible stress interruption. dp is calculated from the elastic modulus and the tensile strength while dc is obtained from the fracture energy. The above mentioned cohesive model can be applied to the fracture of glass, as similar models was employed for the fracture of brittle materials, like ceramics or glass. Camacho and Ortiz [30] developed a Lagrangian finite element method for fracture and fragmentation in brittle materials, with a linear cohesive relation to propagate multiple cracks along arbitrary paths. In Chen et al. [31], a cohesive element was adaptively inserted into the common surface of two adjacent glass elements, which follows an irreversible Camacho–Ortiz cohesive law with linear softening. Holmquist and Johnson [32] also considered that damage softening is necessary as it allows material to soften gradually in lieu of a sudden failure, and linear softening curves are employed in their work. The shape of the descending curve largely depends on the nature of glass. For typical glass, the density is approximately 2500 kg/m3, with Young’s modulus E = 70 GPa and Poisson’s ratio m = 0.2. The tensile strength of glass varies amongst different types of glass. In Poleshko and Goralik [33], a 33.32 MPa tensile strength was measured for Glass MKR-1; Ledbetter et al. [34] suggested that the strength of annealed glass is in the range of 7–28 MPa; Sedlacek et al. [35] indicated that the tensile strength of float glass is in the range of 20–100 MPa; Overend and Zammit [36] also determined the reference tensile strength of glass in the range of 24.3–34.6 MPa with 95% reliable probability based on their computational algorithm. Based on the above references, a tensile strength of 30 MPa assumed in this study is well within the reasonable range. In regard to the fracture energy, a value of 4.0 N/m is reasonable according to [5,37,38]. These parameters result in a 1 lm magnitude critical separation and a very steep strain softening curve. In computational term, according to the power softening law (Eq. (8)) given by Foote et al. [39], an exponential curve with a large exponent n is suitable for the fracture of brittle materials like glass.

 n d ¼ 1 dc ft

r

ð8Þ

As a comparison, Rama Chandra Murthy et al. [40] also recommends a bilinear curve for the softening of brittle materials (Fig. 4), which has a similar shape as the curve used in this study. In the present work, a normalised strain softening curve with parameters a = 1.2, b = 1.0 and c = 1.0 in Eq. (7) is plotted in Fig. 5. By using this decay curve, the bilinear property in [40] is represented in a smooth way. This softening curve has not been used before for the fracture of glass and is a good attempt with the FEM-DEM. When typical glass material parameters (e.g., ft = 30 MPa and Gf = 4.0 N/m) are employed in this curve, a very steep descending of strain softening can be represented with a tiny critical separation distance dc, which guarantees that glass experiences

59

X. Chen et al. / Computers and Structures 177 (2016) 56–68

σ ft

4.1. Validation for two-dimensional glass beams 4.1.1. Elastic behaviour Consider a beam of 2 m long and 20 mm thick with both ends clamped in rigid channels, and subjects to the impact of a 25 mm radius steel ball in the middle. Triangular elements with irregular orientation is meshed into the glass beam so that potential mesh dependency is reduced to the minimum. The total number of elements for the glass beam is 26,976. In this analysis, precise deformation of individual discrete elements is not important as the element size is very small therefore only one finite element is used for each discrete element to minimise the computational cost. Fine meshes with 0.5 mm characteristic size elements are assigned in the central impact effective area to obtain a realistic fracture pattern, while coarse meshes in the far field are set based on the Saint-Venant’s principle. The boundary of the projectile is also meshed with fine elements so that a circular shape and blunt impact can be guaranteed. Time step for this problem is 0.2 ns. This mesh is termed as mesh A so that comparisons can be made in the following sections. Configurations of this problem is shown in Fig. 6, where point A is a point in the middle at the bottom of the beam, and point B is a point on the projectile in contact with the beam. Table 1 lists the material properties of the glass beam, projectile and supports. The projectile hits the glass beam at a velocity of 0.5 m/s and is considered not generating fracture during the simulation. Result from the FEM-DEM program Y (Fig. 7) shows that the beam is intact except some insignificant and very local damage near the impact point, which will not affect the global behaviour of the beam. This suggests that the whole beam can be considered within its elasticity as a whole. The same problem is solved by the FEM program ABAQUS with the same mesh in the FEM-DEM analysis and CPS3 elements. In this study, the deformation of the beam and the velocity of the projectile are examined to show the elastic behaviour of glass and the contact effect exerted onto the projectile. Results obtained from both the FEM-DEM and ABAQUS are given in Figs. 8 and 9.

σ1

O

δc

δ1

δ

Fig. 4. A bilinear strain softening curve adapted from [40].

Normalised bonding stress- z

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

0.6 0.4 0.2 0.0 0.0

0.2

0.4

0.6

0.8

1.0

Normalised crack opening - D Fig. 5. Normalised bonding stress z = r/ft versus the normalised crack opening D = (d  dp)/(dc  dp).

brittle fracture. Results on the fracture of glass in the further sections of this paper are all based on this descending curve. In order for the computational process to be stable, the energy dissipation in the softening zone cannot be too small. If a very fine mesh is used [26], the crack tip stress concentration can be modelled more realistically.

Table 1 Material properties used in the beam analysis.

4. Validation

Density (kg/m3) Young’s modulus (GPa) Poisson’s ratio Tensile strength (MPa) Fracture energy (N/m)

This section provides validation of the computational model for the breakage of glass subjects to impact with the glass fracture model and the FEM-DEM code Y. Examples in both 2D and 3D are included and corresponding dynamic behaviours are investigated.

Glass

Projectile

Support

2500 70 0.2 30 4.0

7800 200 0.3 – –

7800 200 0.3 – –

(a)

Point B

(b)

Point A

Fig. 6. Configurations of the impact problem: (a) overall view and (b) details of mesh A in the impact effective area.

60

X. Chen et al. / Computers and Structures 177 (2016) 56–68

Time (ms)

Damage

0.002 Fig. 7. Fracture response of the beam under impact of

v = 0.5 m/s at t = 0.6 ms.

Y-Displacement (mm)

0.00

0.01

-0.05 FEM-DEM

-0.10

ABAQUS

0.02

-0.15 -0.20 -0.25 0.0

0.04 0.2

0.4

0.6

0.8

1.0

Time (ms) Fig. 8. The vertical displacement curves of point A.

0.1

Y-Velocity (m/s)

-0.1

-0.2

0.2 FEM-DEM

-0.3

Fig. 10. Crack development in the glass beam with time.

ABAQUS

-0.4

-0.5 0.0

0.2

0.4

0.6

0.8

1.0

Time (ms) Fig. 9. The vertical velocity curves of point B.

From Fig. 8, the vertical displacement curves of point A on the glass beam from the FEM-DEM and ABAQUS are almost identical, with relative errors less than 5% in general. The results from the FEM-DEM shows some material delay at the beginning of the simulation, showing that the FEM-DEM is more realistic in dynamic impact simulation than the FEM. Fig. 9 shows the change of vertical velocity of point B on the projectile in contact with the beam from the FEM-DEM and ABAQUS. It can be observed that the two curves are in good agreement. The jumping tendency in the ABAQUS simulation can be attributed to the multiple contacts between the projectile and the target. Since the ABAQUS and the FEM-DEM employ different algorithms for contact problems, such jumping tendency is not observed in the FEM-DEM simulation, whose curve is smooth and shows advantage over the ABAQUS contact simulation. Comparisons between the results from both the FEM-DEM and ABAQUS demonstrates that the combined FEM-DEM together with the glass fracture model is well capable of simulating the responses of glass beams under a low impact velocity when the material is within its elasticity.

4.1.2. Fracture behaviour For the beam investigated in Section 4.1.1, fracture will occur should the impact velocity exceeds a certain level. Fig. 10 shows

Fig. 11. A cone crack obtained from an impact experiment (after Knight et al. [41]).

the fracture development of the same beam investigated in the previous section under an impact velocity of 5.85 m/s using the FEM-DEM. Referring to Fig. 10, the damage initiation, crack propagation and the formation of a Hertzian type cone crack can be observed. Almost immediately after the impact, some local damage is found on the top glass surface (t = 0.002 ms). Later on, horizontal cracks occur (t = 0.01 ms) and form a new fracture surface. With the time elapsing, a Hertzian type cone crack is formed gradually and a central flexural crack is observed throughout the thickness (from t = 0.02 ms to 0.2 ms). Meanwhile, small fragments are found around the contact surface. The projectile is still in contact with the glass beam and has not punched through yet at the end of this simulation. The formation of a cone type crack is well recorded in literatures (e.g., Fig. 11) when brittle glass subjects to blunt impact. In this example, a cone crack is formed in a very short duration by about t = 0.1 ms, which can be considered as ‘‘instantaneous”. These phenomena agree well with the observation of Hertz [42],

61

X. Chen et al. / Computers and Structures 177 (2016) 56–68

Fig. 12. Beam with mesh B: (a) mesh configuration in the impact effective area and (b) cracks and fragments at t = 0.2 ms.

4.1.3. Convergence In this section, mesh B (Fig. 12(a)) is employed for the same beam investigated in Sections 4.1.1 and 4.1.2 to study the convergence of the FEM-DEM and glass fracture model. Elements with characteristic size of 1 mm are meshed into the impact effective area, and the total number of elements for the glass beam is 9952. All other parameters are kept unaltered. Impact results for mesh B is shown in Fig. 12(b). It can be observed that a cone type crack highly similar to the one in Section 4.1.2 is obtained. Further investigations in Fig. 13 shows the change of the total kinetic energy curves for both mesh A and B. Despite some minor differences, the total kinetic energy curves have converged and are almost identical when most cracks are formed, e.g., after t = 0.2 ms. 4.1.4. Post-damage behaviour Zang et al. [16] analysed the impact fracture behaviour of monolithic glass using the spherical discrete element model. In this verification, a 2D FEM-DEM analysis is conducted to examine the post-damage behaviour. The beam, with its two ends clamped, is 101 mm  4.76 mm and meshed with 6080 elements with a characteristic size of 0.4036 mm. Time step is 0.1 ns. The impactor is 6 mm in radius with a mass of 1 kg, and hits the glass at a velocity of 30 m/s. Material properties used in the analysis are given in Table 2. As it is well known that the DEM naturally simulates the fragments and fragmentations. The emphasis of the comparison between the FEM-DEM and the DEM is to investigate the postdamage behaviours of glass, e.g., failure patterns, shape and movement of fragments, relative positions of projectiles, etc. Fig. 14 shows the profile of the damage process from both the FEM-DEM and the DEM in the same scale. It can be concluded that the transient responses in Fig. 14(a) and (b) are similar. In both simulations, projectiles punch the glass and their relative positions reach very good agreement, suggesting that the velocity field of the projectile in both the FEM-DEM and the DEM are quite similar, which further extends the conclusion obtained in Fig. 9 from elastic to post-damage stage. Fragmentations obtained using the FEM-DEM is similar with those from the DEM. Differences between the damage behaviours of the FEM-DEM and the DEM can still be found. Some pieces of large fragments are observed in the FEM-DEM result, while most of the fragments are linked individual discrete elements in the DEM simulation. Since triangular elements are used to mesh the structure fully in the FEM-DEM, a better structural simulation is

100%

Normalised total kinetic energy

showing that the FEM-DEM and the glass fracture model are applicable and accurate for glass fracture modelling.

98% Mesh A

96%

Mesh B

94% 92% 90% 0

0.15

0.3

0.45

0.6

Time (ms) Fig. 13. Normalised total kinetic energy curves for mesh A and B.

Table 2 Material properties used in post-damage analysis.

Density (kg/m3) Young’s modulus (GPa) Poisson’s ratio Tensile strength (MPa) Shear strength (MPa) Fracture energy (N/m)

Glass

Projectile

2500 74.1 0.2 34.6 17.9 4.0

1.105  106 200.0 0.3 – – –

obtained. Triangular elements introduce higher rigidity in the FEM-DEM, while fragmentation is dispersed to a wider range in the DEM analysis due to softer connections between particles. Despite these differences, the effectiveness of the post-damage simulation of the FEM-DEM is well validated.

4.1.5. Laminated glass In [18], an impact test on a laminated glass beam is described and high speed camera records are provided. The laminated glass beam is 200 mm  24 mm  10 mm in dimension. Thicknesses of the upper and lower glass layers are both 10 mm. The PVB interlayer is 4 mm thick. The projectile is 4 mm  4 mm  10 mm with the mass of 1 kg, and hits the beam in the middle at a velocity of 3.13 m/s. The beam is constrained by four fixed supporters, and the geometry of the model is shown in Fig. 15(a). For the laminated glass investigated, the interlayer is assumed to be perfect elastic and interfaces between the glass and the interlayer are considered to be well bonded. Consequently, the fracture of the laminated glass can be entirely attributed to the failure of glass, which is within the theme of this paper. A FEM-DEM simulation is conducted to compare the transient fracture responses with

62

X. Chen et al. / Computers and Structures 177 (2016) 56–68

Fig. 14. Impact fracture behaviour sequences-from top to bottom at t = 0, 0.1, 0.2, 0.3, 0.4, 0.5 ms: (a) FEM-DEM and (b) DEM (after Zang et al. [16]).

(a) Geometry (after [18])

(b) Mesh configuration of the FEM-DEM model Fig. 15. Geometry and mesh configurations of the laminated glass beam.

high speed camera images given in [18]. The triangular mesh distribution of the FEM-DEM modelling is shown in Fig. 15(b) with 10,106 elements of a characteristic size of 1 mm for the laminated glass body. Time step is 0.5 ns. Material properties used in the analysis are taken from [18] and tabulated in Table 3.

Fig. 16 shows the diagrams of two sets of failure sequences from both the FEM-DEM simulation and the high speed camera images in [18], respectively. Each sequence is ordered by the number at its left top starts with t = 0 s (sequence number 1), and the time interval between two consecutive fracture sequences is 100 ls. From the FEM-DEM simulation (Fig. 16(a)), one can find:

63

X. Chen et al. / Computers and Structures 177 (2016) 56–68 Table 3 Material properties for the laminated glass beam analysis.

3

Density (kg/m ) Young’s modulus (GPa) Poisson’s ratio Tensile strength (MPa) Fracture energy (N/m)

Glass

Interlayer

Projectile

Supporter

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 – –

(a) FEM/DEM simulation results

(b) High speed camera records (after Ref. [18])

Fig. 16. Results from the FEM-DEM simulation and the experimental images.

(i) Glass starts to fail just beneath the projectile in the upper layer and cracks propagate along the thickness of the upper layer at t = 100–200 ls (sequence number 2 and 3). (ii) A dominant flexural crack is formed in the lower glass layer slightly left to the projectile which can be observed at t = 300 ls (sequence number 4). (iii) The flexural crack in the lower layer propagates towards the interlayer, and is captured at t = 400 ls (sequence number 5). (iv) The final damage pattern at t = 500 ls (sequence number 6) is given, with small fragments near the contact surface. By comparing Fig. 16(a) and (b), conclusions can be drawn that the results from the FEM-DEM simulation agree very well with the high speed camera records. Cracks in the upper glass layer are observed both in the FEM-DEM simulation and the experiment. For the fracture of lower glass layer, dominant flexural cracks are formed at similar locations and time points. The comparison between the simulation and the high speed camera records further demonstrates that the FEM-DEM together with the fracture model employed in this paper is capable of simulating the failure of glass, and is capable of yielding correct results.

Fig. 17. Mesh configurations of the plate impact problem.

Table 4 Material properties used in the glass plate analysis.

Density (kg/m3) Young’s modulus (GPa) Poisson’s ratio Tensile strength (MPa) Fracture energy (N/m)

Glass

Projectile

2500 70.0 0.23 30.0 4.4

7800 200.0 0.3 – –

4.2. Validation for three-dimensional glass plate analysis Pauw [43] performed ABAQUS simulations and experiments for impact breakage of a circular glass plate with the radius R = 50 mm and thickness h = 4 mm. The plate is clamped around the entire circumference. A 33 mm high hemisphere cylinder (bullet) with the radius r = 10 mm hits the centre of the glass plate at a velocity of 1.98 m/s. In the 3D FEM-DEM modelling, 4358 four-node constant strain tetrahedral elements with the characteristic size of 4 mm are used to mesh the glass plate. Time step is 5 ns. Mesh configurations for this problem are shown in Fig. 17.

Material properties used in the analysis are tabulated in Table 4. In this section, the final damage pattern is investigated as it is the reflection of stress wave distribution. For a particular impact problem, correct simulations should yield correct fracture patterns on a macroscopic level. The fracture patterns obtained from both the FEM-DEM and [43] at t = 10 ms are given in Fig. 18. It can be observed that damages are mainly localised around the centre of the plate and the projectile starts to penetrate the glass plate. The simulation from the FEM-DEM agrees well with that

64

X. Chen et al. / Computers and Structures 177 (2016) 56–68

Fig. 18. Damage patterns of the plate at t = 10 ms: (a) FEM-DEM and (b) ABAQUS (after Pauw [43]).

from ABAQUS. In [43], although the centre of the glass plate is meshed with finer elements, damage is still localised and no crack propagates to the edge of the plate. Current results from the FEMDEM simulation is also in good agreement with the experimental damage pattern given in Fig. 19, where a hole is formed in the centre of the specimen. 4.3. Summary

5. 2D parametric study In Section 4, the validity of the use of the combined FEM-DEM for modelling two and three-dimensional glass impact breakage is established. This section reports the parametric study performed to check the stability of the results obtained from the FEM-DEM and provides comparative benchmarks for future research.

Fig. 19. Experimental damage pattern of the glass plate (after [43]).

2.5

2D

2.0

3D

Time (s)

Numerical simulations with both 2D and 3D examples are performed with the FEM-DEM computer code Y and a Mode I glass fracture model. In the 2D analysis, elastic and fracture behaviours of glass beams subject to impact are investigated. Elastic analysis in Section 4.1.1 shows that the FEM-DEM produces as good results as that from the FEM when material is within the elastic range. Both the time evolution of the deformation of the beam and the velocity of the projectile demonstrates that they are in good agreement with ABAQUS results. A highly similar cone crack is obtained in Section 4.1.3 when the mesh is changed. Apart from that, the total kinetic energy curves from different meshes converge well. Postdamage simulation in Section 4.1.4 shows similar fracture and fragmentation as that from the DEM, leaving the fragments in free kinetic movement. An impact simulation on a laminated glass beam in Section 4.1.5 demonstrates that the FEM-DEM together with the fracture model yield excellent agreement with high speed camera records. In the 3D, simulation of a clamped circular glass plate subjects to bullet impact is investigated and qualitatively validated with results from ABAQUS and experiment. Comparisons between the results of the FEM-DEM and ABAQUS shows that a central punching is observed, no matter whether finer mesh is employed at the centre or not. In summary, the intrinsic discrete formulation makes the FEMDEM simulation of discontinua (e.g., fracture and fragmentation) a natural result. The implemented Mode I fracture model is applicable and accurate enough to describe crack initiation and propagation in the glass body.

1.5 1.0 0.5 0.0 0

5000

10000

15000

20000

25000

30000

Number of elements Fig. 20. The computational time needed for each step versus number of elements.

Currently, the serial FEM-DEM code takes a relatively long time to run should there are many elements involved in the computation. Fig. 20 shows the computational efficiency for both 2D and 3D problems of the FEM-DEM program. The computations were performed on an Intel core 2 Duo 2.66 GHz desktop with 2 GB memory and Windows XP operating system. For practical purpose, computational time increases almost linearly with the number of elements, and the computational time per element for 3D is about 60 times as that in 2D for each step. To improve the computational efficiency, a parametric study is performed for 2D and conclusions can be similarly extended to 3D problems as direct investigation on 3D is computationally

65

X. Chen et al. / Computers and Structures 177 (2016) 56–68

expensive. The study is based on the glass beam investigated in Section 4.1.1 by varying the tensile strength ft, fracture energy Gf and impact velocity v. Simulations are terminated at t = 0.6 ms as most cracks, fragments and fragmentations are formed by that time.

σ

σ

f2 f1

f1 f3

5.1. Tensile strength According to Elices et al. [44], tensile strength ft and fracture energy Gf play important roles in the cohesive glass fracture model discussed in Section 3. Taking a linear softening curve as an illustration, Fig. 21 shows that both parameters affect the shape of the softening curve substantially. The solid line in Fig. 21 represents the original softening curve with the tensile strength f1 and the critical separation dc1. As shown in Fig. 21(a), with the fracture energy be held constant, a higher tensile strength raises the starting tensile strength of the softening curve from (O, f1) to (O, f2) and decreases the critical separation from dc1 to dc2, resulting in an even steeper softening curve, and vice versa. Fig. 22 shows the fracture patterns of the beam of different tensile strengths with fracture energy Gf = 4 N/m and impact velocity v = 5.85 m/s at t = 0.6 ms. It can be observed from the figure that for a lower tensile strength (e.g., ft = 20 MPa), flexural cracks are obtained and the material behaves in a more ductile manner. On the other hand, the higher the ft, the more brittle the beam is, with no observation of any obvious flexural failure but a cone type crack instead. However, it can also be observed that less local damage occurs due to a higher tensile strength.

O

δ c1 δ c3

O

(a)

δ c1

δ c4

(b)

Fig. 21. Influences of material parameters on softening curves: (a) tensile strength and (b) fracture energy.

ft (MPa)

Damage

20

50

5.2. Fracture energy The influences of fracture energy Gf on the breakage of glass are simpler. From Fig. 21(b), a higher Gf produces a larger critical separation dc4 if the tensile strength ft is left unaltered, pushing the stress going a flatter softening path during the damage evolution process. Fig. 23 lists the fractures and fragmentations of the beam over different fracture energy Gf with the tensile strength ft = 30 MPa and the impact velocity v = 5.85 m/s at t = 0.6 ms. It can be understood from Fig. 23 that generally materials with small Gf exhibit brittleness while large Gf increases the ductility. When Gf is low (e.g., 3–5 N/m, which is a reasonable range for float glass according to a variety of literatures [37,38]), a cone type crack is usually obtained. Once Gf becomes larger, a central throughthickness crack appears. The disappearance of a cone crack can be attributed to the increase of the fracture energy, which makes some element boundaries more resistant to separate and only a dominant central flexural crack is observed. The critical fracture energy which shifts the failure mode from cone to flexural crack is 5–6 N/m.

δ c2

70

80 Fig. 22. Breakages of the glass beam with Gf = 4 N/m and v = 5.85 m/s at t = 0.6 ms.

impact velocity v reaches a certain level, e.g., 5 m/s in this study. For moderately high impact velocities (5–10 m/s), a classical cone type crack is found. The higher the impact velocity, the more severe the fragmentation will be. Cone type cracks are still recognisable when the impact velocity is below 20 m/s. Any impact with v > 20 m/s results in crushing, leaving the impact effective area failing quickly as the amount of input energy completely overwhelms the resistance of glass.

5.3. Impact velocity

5.4. Discussions

Impact velocity exhibits a direct influence on the breakage behaviour of glass. Simply referring to the equation of kinetic energy Ek = Mv2/2, where M is the mass of the projectile and v is the impact velocity, a slight change in the impact velocity changes the input energy dramatically. In this sub-section, velocities below 50 m/s are investigated as this is the basic wind speed in wind zone one (110–120 mph) according to ASTME1996. Fig. 24 schematically shows the glass breakage under a variety of impact velocities ranging from 0.5 m/s to 50 m/s with ft = 30 MPa and Gf = 4 N/m at t = 0.6 ms. For very low impact velocities (e.g., v 6 1 m/s), damages are restricted to the local area around the contact point. With the increase of v, central flexural cracks are observed. A cone fracture does not appear until the

5.4.1. Index of brittleness It is understood from Sections 5.1 and 5.2 that for glass with a given thickness h, the higher the tensile strength ft (or the smaller the fracture energy Gf), the more brittle the material is. Define D to be the characteristic size of the smallest element, based on this conclusion, a parameter j is proposed for evaluating the material brittleness, where

j/

ft D Gf

ð9Þ

Should the impact velocity is relatively high, cone and flexural cracks are usually observed. The shift from cone to flexural cracks indicates that the material becomes ductile. Further investigation

66

X. Chen et al. / Computers and Structures 177 (2016) 56–68

Gf (N/m)

Table 5 Estimate of threshold tensile strengths between cone and flexural cracks.

Damage

3

Gf (N/m)

Transit ft (MPa)

j (106 D)

3 4 5 6

<20 20–30 20–30 >30

<6.67 5–7.5 4–6 >5

on fracture behaviours gives the estimate threshold tensile strengths between cone and bending cracks for different fracture energy Gf in Table 5. It can be observed from Table 5 that in order to keep j a constant level (should in the range of (5–6)  106 D), a higher ft requires a higher Gf to compensate. Thus the expression given in Eq. (9) is verified and can be used as a quick check for brittleness.

5

6

8 Fig. 23. Breakage of the glass beam with ft = 30 MPa and

v = 5.85 m/s at t = 0.6 ms.

5.4.2. Breakage regimes The trend of breakage patterns under different impact velocities can be determined. Very localised damage occurs at the surface around where the projectile and the target in contact when the impact velocity is very low. Then central flexural fracture, which can also be termed as ‘bending’, is followed by cone type cracks. Higher impact velocities introduce crushing, making any regular crack difficult to be observed. Consequently, four breakage regimes are obtained from Fig. 24, they are: (1) local minor fracture, (2) central bending crack, (3) cone type crack and (4) punching. From

Fig. 24. Breakage of the glass beam with ft = 30 MPa and Gf = 4 N/m at t = 0.6 ms from top to bottom: (a)

v = 0.5, 1, 2, 3, 4 m/s and (b) v = 5, 10, 20, 30, 50 m/s.

X. Chen et al. / Computers and Structures 177 (2016) 56–68

3200

Impact kinetic energy (J)

Minor-Bending Bending-Cone

Punching

2400

Cone-Punching

1600

Cone

800

0 0

0.02

0.04

0.06

0.08

0.1

Fracture resistance (N)

(a) Four damage regimes in one sacle

Impact kinetic energy (J)

80

Minor-Bending 60

Bending-Cone

Cone

40

Bending 20

Minor 0

0

0.02

0.04

0.06

0.08

0.1

Fracture resistance (N)

(b) Minor and bending damage regimes only Fig. 25. Breakage regimes and thresholds.

the energy point of view, a diagram indicating the shift of regimes is given in Fig. 25(a), where minor and bending damage regimes are so small that details are shown separately in Fig. 25(b). Note that by using the product of fracture energy and thickness Gfh (N) as the horizontal axis, and the input kinetic energy Mv2/2 (J) as the vertical axis, the relationship between the glass fracture resistance, input energy and breakage scenario can be established. In Fig. 25, the area between two adjacent curves is the regime represents a particular type of breakage (e.g., bending, cone or punching). With the increase of Gfh, the input energy Mv2/2 needs for producing a certain type of damage increases accordingly. Among these regimes, local minor fracture occupies the smallest area, which is almost unnoticeable in Fig. 25(a), and bending cracks also exist in a narrow band; a great amount of area is given to the cone type cracks, while there is no upper limit for punching. 6. Conclusions The combined finite-discrete element method is employed in this paper for glass breakage analysis under hard body impact. A Mode I cohesive crack model similar to Hillerborg’s is used for the fracture analysis of glass. Results from the FEM-DEM together with the employed glass fracture model are validated and found to be in good agreement with existing computational and experimental data both in 2D and 3D. (i) Comparison with the ABAQUS simulation in Section 4.1.1 shows that the FEM-DEM is capable of simulating the responses of a glass beam under low impact velocity when the material is within the elasticity.

67

(ii) A Hertzian type cone crack, which is widely known in existing literatures, is simulated by the FEM-DEM in Section 4.1.2 when the impact velocity is increased. (iii) Highly similar cone cracks and kinetic energy curves are obtained in Section 4.1.3 when applying different meshes to the FEM-DEM simulations, demonstrating that the results obtained has converged. (iv) Comparison with the DEM simulation in Section 4.1.4 shows the effectiveness of the post-damage simulation of the FEMDEM, and its advantage in achieving realistic crack patterns and fragments over the particle-based DEM. (v) Comparison with the high speed camera records in Section 4.1.5 shows that dominant flexural cracks are formed at similar locations and time points, which further demonstrates that the FEM-DEM together with the fracture model employed in this paper is capable of yielding correct impact simulation results. (vi) Comparison with the 3D ABAQUS simulation and experiment in Section 4.2 shows similar damages are obtained in the FEM-DEM simulation. The projectile penetrates the plate with obvious material removal, leaving a hole in the centre. Therefore, conclusions can be reached that the FEM-DEM together with the glass fracture model is applicable, reliable and robust to glass impact analysis with the considerations of fracture and fragmentations. A parametric study is performed by varying the tensile strength, fracture energy and impact velocity. Tensile strength is a starting point for the fracture commencement, and fracture energy exhibits the control over the damage process. The ratio j of tensile strength over fracture energy can be used as a quick index to evaluate the material brittleness. Performance of glass under impact can be greatly improved by raising the tensile strength as well as the fracture energy. Different glass breakage regimes and thresholds are established based on the impact-resistance diagram, providing comparative benchmarks for future research. Though glass is investigated in this work, the combined FEMDEM and employed fracture criterion can be further developed to model the breakage behaviour of other brittle or quasi-brittle materials. Acknowledgements The present work was financially supported by the China Scholarship Council (No. 2009621028), the University of Birmingham (PGTA–Postgraduate Teaching Assistantship), Jiangsu University of Science and Technology (No. 1732921402) and Jiangsu ShuangChuang Program. The combined FEM-DEM code Y and technical supports from Professor Antonio Munjiza, Queen Mary University of London are highly appreciated. References [1] Axinte E. Glasses as engineering materials: a review. Mater Des 2011;32 (4):1717–32. [2] Preston FW. A study of the rupture of glass. J Soc Glass Technol 1926;10:234–69. [3] Murgatroyd JB. The significance of surface marks on fractured glass. Soc Glass Technol 1942;26:153–5. [4] Shand EB. Experimental study of the fracture of glass. I. The fracture process. J Am Ceram Soc 1954;37(2):52–9. [5] Holloway DG. The fracture of glass. Phys Educ 1968;3:317–22. [6] Clarke DR, Faber KT. Fracture of ceramics and glasses. J Phys Chem Solids 1987;48(11):1115–57. [7] Ward JB. Analysis of glass fracture. Mater Des 1987;8(2):100–3. [8] Toakley AR. Stresses and safety levels for glass liable to human impact. Build Environ 1977;12:87–95. [9] Minor JE. Windborne debris and the building envelope. J Wind Eng Ind Aerodyn 1994;53:207–27.

68

X. Chen et al. / Computers and Structures 177 (2016) 56–68

[10] Sun X, Khaleel MA, Davies RW. Modelling of stone-impact resistance of monolithic glass ply using continuum damage mechanics. Int J Damage Mech 2005;14:165–78. [11] Dharani LR, Yu J. Failure modes of glass panels subjected to soft missile impact. In: Brebbia CA, Varvani-Farahami A, editors. Damage and fracture, vol. VIII. WIT Press; 2004. p. 163–71. [12] Dharani LR, Wei J, Yu J, Minor JE, Behr RA, Kremer PA. Laminated architectural glass subjected to blast impact loading. Am Ceram Soc Bull 2005;84(1):42–4. [13] 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. [14] Oda J, Zang MY, Mori T, Tohyama K. Simulation of dynamic fracture behavior of laminated glass by DEM. In: Trans 8th Calculation Dynamics Symp JSME; 1995. p. 429–30. [15] Oda J, Zang MY. Analysis of impact fracture behaviour of laminated glass of bilayer type using discrete element method. Key Eng Mater 1998;145– 149:349–54. [16] 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. [17] 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. [18] 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. [19] Munjiza A. The combined finite-discrete element method. John Wiley and Sons; 2004. [20] Munjiza A. Discrete elements in transient dynamics of fractured media. PhD thesis. Swansea: Department of Civil Engineering, University of Wales; 1992. [21] Munjiza A, Rougier E, Knight EE. Large strain finite element method: a practical course. John Wiley and Sons; 2015. [22] Lisjak A, Grasselli G. Rock impact modelling using FEM/DEM. In: Munjiza A, editor. DEM5; 2010. p. 269–74. [23] Munjiza A, Lei Z, Divic V, Peros B. Fracture and fragmentation of thin shells using the combined finite-discrete element method. Int J Numer Meth Eng 2013;95:478–98. [24] Munjiza A, Andrews KRF, White JK. Combined single and smeared crack model in combined finite-discrete element analysis. Int J Numer Meth Eng 1999;44 (1):41–57. [25] Williams JR. Contact analysis of large numbers of interacting bodies using discrete modal methods for simulating material failure on the microscopic scale. Eng Comput 1988;3:197–209.

[26] Munjiza A, John NWM. Mesh size sensitivity of the combined FEM/DEM fracture and fragmentation algorithms. Eng Fract Mech 2002;69:281–95. [27] Anderson TL. Fracture mechanics: fundamentals and applications. 2nd ed. Boca Raton: CRC Press; 1991. [28] Wang EZ, Shrive NG. On the Griffith criteria for brittle fracture in compression. Eng Fract Mech 1993;46(1):15–26. [29] Hillerborg A, Modéer M, Petersson PE. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cem Concr Res 1976;6(6):773–81. [30] Camacho GT, Ortiz M. Computational modelling of impact damage in brittle materials. Int J Solids Struct 1996;33:2899–938. [31] 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. [32] 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. [33] Poleshko AP, Goralik ET. Tensile strength of technical ceramics and glass. Strength Mater 1981;13(3):344–7. [34] Ledbetter LR, Walker AR, Keiller AP. Structural use of glass. ASCE – J Archit Eng 2006;12(3):137–49. [35] Sedlacek G, Blank K, Gusgen J. Glass in structural engineering. Struct Eng 1995;73(2):17–22. [36] Overend M, Zammit K. A computer algorithm for determining the tensile strength of float glass. Eng Struct 2012;45:68–77. [37] Linger KR, Holloway DG. The fracture energy of glass. Philos Mag 1968;18 (156):1269–80. [38] Michio I, Kazuhiro U, Shinsuke T, Yasuo G, Mototsugu S. Work of fracture and crack healing in glass. J Am Ceram Soc 2006;68(12):704–6. [39] Foote RML, Mai Y-W, Cotterell B. Crack growth resistance curves in strainsoftening materials. J Mech Phys Solids 1986;34(6):593–607. [40] Rama Chandra Murthy A, Palani GS, Lyer NR. State-of-the-art review on fracture analysis of concrete structural components. Sadhana 2009;34 (2):345–67. [41] Knight CG, Swain MV, Chaudhri MM. Impact of small steel spheres on glass surfaces. J Mater Sci 1977;12(8):1573–86. [42] Hertz H. On the contact of elastic solids. London: Macmillan; 1896. [43] Pauw SD. Experimental and numerical study of impact on window glass fitted with safety window film. PhD Thesis. Ghent University; 2010. [44] Elices M, Guinea GV, Gomez J, Planas J. The cohesive zone model: advantages, limitations and challenges. Eng Fract Mech 2002;69:137–63.