Impact Fracture and Fragmentation of Glass via the 3D Combined Finite-Discrete Element Method ^†

Abstract

A driving technical concern for the automobile industry is their assurance that developed windshield products meet Federal safety standards. Besides conducting innumerable glass breakage experiments, product developers also have the option of utilizing numerical approaches that can provide further insight into glass impact breakage, fracture, and fragmentation. The combined finite-discrete element method (FDEM) is one such tool and was used in this study to investigate 3D impact glass fracture processes. To enable this analysis, a generalized traction-separation model, which defines the constitutive relationship between the traction and separation in FDEM cohesive zone models, was introduced. The mechanical responses of a laminated glass and a glass plate under impact were then analyzed. For laminated glass, an impact fracture process was investigated and results were compared against corresponding experiments. Correspondingly, two glass plate impact fracture patterns, i.e., concentric fractures and radial fractures, were simulated. The results show that for both cases, FDEM simulated fracture processes and fracture patterns are in good agreement with the experimental observations. The work demonstrates that FDEM is an effective tool for modeling of fracture and fragmentation in glass.

1. Introduction

It is a well-known fact that a great number of traffic accidents result in windshield breakage that unduly causes great harm to automobile to both passengers and pedestrians. Due to the consequences of these accidents, a forensics approach was developed wherein analysts can use windshield fracture patterns to reconstruct some aspects of accidents [1]. Needless to say, the study of automobile glass impact fracture mechanisms is of theoretical and practical importance to the automobile industry as it affects passenger protection, passive safety measures, and traffic accident reconstruction.

Current automobile windshield safety analysis standards heavily rely upon both experimental and analytical approaches to meet Federal guidelines for breakage standards. Due to computational mechanic’s advancements, numerical methods offer an alternative approach that is proving to be effective for the study of automobile glass fracture. Given the critical nature and need for computational methods that can effectively capture glass breakage many efforts have been undertaken. Bios et al. simulated the impact of a sphere into a glass plate and the behavior of a windscreen during a roof crash via the finite element method (FEM), where failed elements were deleted from the calculation after the strain was greater than a pre-set failure strain [2]. Timmel et al. presented a computational technique for the modeling of laminated glass using an explicit finite element code. In their work, different material models, element techniques, and the influence of the mesh were discussed [3]. Xu et al. used the extended finite element method (XFEM) to characterize radial crack and circumferential crack propagations of windshield cracking under low-speed impacts [4]. Pyttel et al. proposed a failure criterion for laminated glass under impact loading [5]. The criterion was implemented into a commercial finite element code and was validated by comparing against experiments. Peng et al. simulated the mechanical behavior of a windshield-laminated glass given the impact of a pedestrian’s head using an element deletion approach implemented in a commercial finite element code [6]. Different finite element models were tested and the numerical results were verified against experimental data. Chen et al. utilized a cohesive zone model based computational framework for the modeling of impact fracture on automotive laminated glass [7]. Xu et al. investigated radial multi-cracking phenomenon in laminated glass subject to dynamic loading using XFEM [8]. Lin et al. modeled automotive windshield impact fracture behavior via an intrinsic cohesive approach [9]. The numerical results were compared against the element deletion method and experimental observations. In all the above mentioned work, the simulations were conducted in the context of the finite element method (FEM).

In the meantime, computational mechanics of discontinuum approaches have also been extensively utilized to simulate the glass impact fracture processes. Oda et al. simulated the dynamic fracture behavior of laminated glass using a 2D discrete element method (DEM) implementation, where both the glass and polyvinyl butyral (PVB) film were divided into identical circular discrete elements [10,11]. Zang et al. investigated the impact fracture behavior of automobile glasses using a 3D particle-based DEM, wherein the advantage of laminated glass in passenger safety was thoroughly demonstrated [12]. In order to simulate the large deformation of the PVB layer in laminated glass, Lei and Zang developed a numerical framework that combined particle-based DEM and explicit FEM, where the mechanical response accounted for glass fracture by using DEM while the large deformation of the PVB was modeled using FEM [13]. Lei then studied the impact fracture mechanisms of automobile glass using the 3D combined finite-discrete element method (FDEM), where the community accepted combined single smeared crack model was extended to 3D for the modeling of mode I and II fractures [14]. Munjiza et al. proceeded by developing a model for fracture and fragmentation of multi-layered thin shells in the context of FDEM [14,15]. In their work, the impact fracture patterns for a thin flat glass shell and a thin spherical glass shell were simulated. Chen and Chan simulated the fracture and fragmentation responses of laminated glass under hard body impact using FDEM [16], where different fracture patterns (e.g., cone and flexural) were simulated in both 2D and 3D. Xu et al. proposed a 3D adaptive algorithm, which automatically converts distorted finite elements into spherical discrete elements, for the simulation of impact fracture of laminated glass [17]. Wang et al. compared four different numerical methods (i.e., FEM, XFEM, DEM, FDEM) for the fracture of brittle materials with specific reference to glass [18]. They concluded that FDEM yields the most satisfactory performance for the modeling of the dynamic fracture of materials.

In this work, FDEM was used to simulate dynamic impact fracture processes of a laminated glass beam and a glass plate in 3D. The enriched fracture details, such as the fracture processes’ time sequences of the laminated glass beam and the concentric and radial fractures of the glass plate, were compared against the experimental observations. The rest of the paper is organized as follows: a brief overview of FDEM is introduced in order to provide the reader with a general framework of the method; a generalized traction-separation model for modeling fracture and fragmentation is then introduced; a laminated glass impact fracture simulation is then presented; finally, FDEM glass plate fracture pattern phenomenology is discussed.

2. Overview of the Combined Finite-Discrete Element Method

FDEM is an effective tool for addressing a variety of physics problems formulated not in terms of the continuum assumption and differential equations, but in terms of a large number of discrete entities interacting with each other [19,20,21]. In FDEM the solid domains (called discrete elements) are discretized into finite elements, where finite rotations and finite displacements are assumed a priori. Through failure, fracture, and fragmentation, single domains represented by separate finite element meshes are transformed into a number of interacting domains. The finite element discretization of the solid domains is also conveniently used to discretize the contact between discrete elements. Utilizing this approach, discretized contact solutions can then be used for both contact detection and contact interaction [22].

The generalized governing equation of a FDEM system is

$$M\ddot{x}+C\dot{x}=f$$

(1)

where $M$ is the lumped mass matrix, $C$ is the damping matrix, $x$ is the displacement vector, and $f$ is the equivalent force acting on each node and includes all forces existing in the system such as the body forces, boundary tractions, forces due to material deformation as well as contact forces between solid domains and cohesion forces in the damaged areas [20]. Equation (1) is then integrated in time in order to obtain the transient evolution of the system. There are several time integration schemes that can be used for this purpose. In this work the central difference time integration scheme was adopted [23].

From an algorithmic point of view, FDEM includes: material deformation, contact detection, contact interaction, and continua-discontinua transition [20]. An in-depth description of the FDEM is outside of the scope of this paper. The interested reader is referred to the following seminal works for more detailed descriptions of the method [20,22,24].

2.1. Material Deformation

In FDEM, large strain finite element methods are used to simulate material deformation. In early versions of FDEM, only the elastic deformation of the solid material was taken into account. Hyper-elastic material constitutive laws, which defined stress as a function of finite strains (e.g., the Green-St. Venant strain), were implemented to simulate solid material deformation under finite rotations and finite displacements [20].

For more recent versions of FDEM, the solid deformation is calculated using a multiplicative decomposition-based formulation [24]. This approach naturally decomposes deformation into translation, rotation, plastic stretches, elastic stretches, volumetric stretches, shear stretches, etc. In essence, a total deformation description is obtained from the displacement field via a decomposition of the respective deformation functions; which, when derivation is applied, results in multiplication, thus the term multiplicative decomposition [24]. Of note, the multiplicative decomposition-based formulation has been applied to define material constitutive laws in different ways, from anisotropic elastic formulation for rock materials [25,26], to plastic formulation for metals [27] and for a generalized elastoplastic formulation for anisotropic solids [28].

In terms of mesh “element” technology, domain discretization has usually been conducted by implementing constant strain elements (i.e., constant strain triangles and tetrahedrons) in FDEM [20]. However, it is well-known that these lower order elements can experience numerical locking. This occurs when two or more physical mechanisms compete for the available degrees of freedom within each finite element, especially in incompressibility conditions [24]. Recently, in order to alleviate this deficiency, composite triangular and tetrahedral finite elements were developed with a selective integration scheme to properly relax the constraints on the volumetric term of the material deformation [25,26]. It is noted here that since the volumetric deformation is calculated separately from other stress-bearing mechanisms, it is relatively easy to construct a selective integration scheme in the composite elements that have a multiplicative decomposition-based formulation. Moreover, the composite finite elements, which are constructed with a group of low order sub-elements, are good choices for FDEM since only the information from the sub-elements is needed. The approach is especially convenient for cases where re-meshing is necessary (e.g., dynamical fracture propagation). In addition, the low order sub-elements enable robust contact interaction algorithms that maintain both a relatively high computational efficiency and accuracy as well.

2.2. Contact Detection

In the context of FDEM, contact enforcement is usually conducted in two steps, contact detection and contact interaction [20]. The goal of contact detection is to determine the contact relationship between each finite element and its adjacent elements, while contact interaction algorithms calculate the contact forces between the elements in contact. In typical FDEM applications, the simulation system can have thousands to millions of discrete elements with irregular shapes which can freely move and rotate in space until they are in contact with other elements. As a result, random collisions can occur between any free faces of the discrete elements at any moment. In order to correctly resolve contact relationships in such a complex multibody system, contact detection algorithms are required to have good robustness, be highly efficient and require low memory usage.

Most of the modern contact detection algorithms can be classified either as tree-based or grid-based algorithms [20]. In the tree-based algorithms, the position and size of each contact body (discrete element or finite element) is represented using a special tree structure, while grid-based algorithms build a data structure by mapping each contact body onto identical background grids (bins). Both groups of algorithms usually contain sort and search components: the data structure is first sorted via selected sorting algorithms and the final contact search is then conducted on the sorted data structure. Han et al. compared the performance of several selected tree-based and grid-based contact algorithms and showed that selected grid-based algorithms can be 100× times faster than tree-based algorithms [29].

The NBS algorithm is a linear contact detection algorithm that was first proposed by Munjiza for bodies of similar size [30]. Since NBS features linear sort and linear search algorithms, the total detection time is linearly proportional to the number of contact bodies in the system. Moreover, a special linked-list structure is used in NBS which guarantees that memory usage is nearly proportional to the number of contact bodies. Other grid-based algorithms with linear complexity include the MR and CGRID algorithms, of which both were built on top of the NBS algorithm. The memory usage in MR is exactly proportional to the number of contact bodies [31], while computational efficiency and memory usage in CGRID are insensitive to the size of the contact bodies [32].

For more recent versions of FDEM, a contact detection algorithm called MRCK was developed to improve contact detection performance [22]. In MRCK, the contact targets are sorted using the MR linear sort which takes advantage of the fact that no contact target can move more than the size of a single cell in a single time step. For a given contactor, a process called “contactor rendering” is used to detect all the cells that the contactor currently occupies. The rendering is done in conjunction with the sorted list of targets in such a manner that only the cells that have targets mapped to them are rendered. MRCK significantly speeds up the rendering process since a contactor with no contacts is not rendered at all.

2.3. Contact Interaction

From its inception, FDEM has used a discretized distributed potential contact force algorithm to resolve interaction between contact bodies (discrete elements and/or finite elements) [20,33]. In the earliest versions, the “element to element” contact interaction, such as “triangle to triangle” in 2D and “tetrahedron to tetrahedron” in 3D, was implemented [20]. This “element to element” approach exactly considers the geometry of both the contactor and the target elements and the integration of the contact forces distributed along the edges/faces of these contact bodies was done analytically. Since this approach integrated contact forces exactly, it was therefore quite time consuming.

In the latest version of FDEM, the contact interaction has been simplified by using “element to point” contact interaction kinematics [22], e.g., “triangle to point” in 2D and “tetrahedron to point” in 3D. For these cases, the target elements are discretized into a series of contact points distributed on the free boundary edges/faces of the discrete elements. Each target point is considered a Gauss integration point through which distributed interaction forces are integrated. This simplified “element to point” contact interaction is much faster than the original “element to element” approach as it integrates the distributed contact forces in an approximate form.

In the discretized distributed potential contact force algorithm, the amplitude of the contact force is usually a function of the contact potential, while the direction of the contact force is perpendicular to the contour of the contact potential field [34]. As a result, the contact forces calculated through this method rely on the evaluation of the contact potential field. In the earlier versions of FDEM, the potential field was defined locally, according to the geometry of each finite element, which introduced an artificial numerical non-smoothness in the contact force, i.e., both the jump in the amplitude and direction of the contact force could be observed when the contact point moved from one finite element to another [34]. In order to overcome the non-smoothness in this contact force, a smooth contact algorithm was recently developed [34]. In the new smooth contact algorithm, the contact potential is defined globally as the geometrical information of each discrete element is accounted for via nodal connectivity and existing discrete element boundaries, thereby yielding smooth potentials as well as contact forces [34].

In its original form, for contact processes the energy balance calculated from the discretized distributed potential contact force algorithm was always preserved. For the recently developed generalized contact interaction law, different mechanisms that dissipate the energy during the contact process were introduced in the context of the discretized distributed potential contact force algorithm [34].

2.4. Continua-Discontinua Transition

In FDEM, cohesive zone models are usually used to simulate the fracture and fragmentation in solids [20]. In the cohesive zone models, the mechanical response of the solid material is decomposed into a solid matrix part and an interfacial part that represents multi-scale processes taking place behind the fracture front (e.g., breaking of bonds between grains). In the solid matrix part, the finite element method along with constitutive laws (e.g., elastic, plastic, continuum damage models) are used to simulate the solid’s bulk deformation, while fracture and fragmentation are handled in the interfacial part via cohesive elements or cohesive points that are assumed to coincide with the finite element boundary [20]. The mechanical response of the cohesive elements/points are dominated by the so-called “traction-separation” constitutive law. To control when and how the cohesive elements/points are introduced into the system, different cohesive zone models have been developed and implemented in FDEM.

One group of cohesive zone models is called the intrinsic cohesive zone model [35]. In this model, the cohesive elements/points are inserted in advance between two adjacent elements to connect those two originally adjacent elements. In order to enforce the continuity of the material, “strain-hardening” with a sufficiently large initial slope must be introduced in the “traction-separation” curve, which reduces the stiffness of the system artificially. Another group of cohesive zone models are called the extrinsic cohesive zone model [36]. In the extrinsic models, the cohesive elements are dynamically inserted into the system according to the stress state. The extrinsic models avoid the artificial compliance seen in the intrinsic model approach, however they may end up being “time-discontinuous” at the point when the material is transitioning from continua to discontinua.

In the earlier versions of FDEM, the combined single smeared crack model [37,38], which belongs to the family of intrinsic cohesive zone models, was developed for the simulation of fracture and fragmentation. Originally, the combined single smeared crack model was aimed at mode I loaded cracks only and it was implemented in 2D [37]. In 2010, that crack model was extended to 3D for both mode I and mode II fractures and this extended 3D crack model was implemented in the open-source FDEM code Y-code [14,39]. Since then, with the open distribution of Y-code, this 3D crack model has been widely used in different research groups for different topics.

The latest version of FDEM features a unified cohesive zone model [40]. Similar to the traditional extrinsic cohesive zone model, the unified cohesive zone model dynamically inserts the cohesive element into the system based on the local stress state. However, the transition from continua to discontinua is smoothly achieved and as such, the unified cohesive zone model eliminates the “time-discontinuous” issue. Moreover, within the unified cohesive zone model the point of transition from continua to discontinua is controllable. Both the intrinsic cohesive zone model and enhanced extrinsic cohesive model can be retrieved from the unified model through the introduction of a threshold parameter [40]. Thus, the unified cohesive zone model, which unifies the intrinsic and extrinsic cohesive zone models, has all the advantages of both approaches while overcoming most of the disadvantage of existing cohesive zone models [40].

3. Generalized Traction-Separation Model for Fracture and Fragmentation

In the cohesive zone models, the bonding stresses along the damage zone or discrete crack are defined as functions of the local separation (displacement). As shown in Figure 1a, the local separation at any point on the surfaces of a crack can be divided into three components

$$\mathbf{\delta }={\delta }_{n}n+{\delta }_{t_1}{t}_1+{\delta }_{t_2}{t}_2$$

(2)

where $n$, $t_1$, and $t_2$ are the normal and tangential material axes which define a local coordinate systems moving and rotating with the material, while ${\delta }_n$, ${\delta }_{t_1}$, and ${\delta }_{t_2}$ are the normal and tangential separations at any point in the damage zone, respectively.

Accordingly, the traction vector $p$ is also divided into three components with respect to the material axis $n$, $t_1$, and $t_2$ (Figure 1b),

$$p=\sigma n+{\tau }_1{t}_1+{\tau }_2{t}_2$$

(3)

where $\sigma$, ${\tau }_1$, and ${\tau }_2$ are the normal and tangential stresses in the direction of $n$, $t_1$, and $t_2$.

One approach that can be used to construct a relationship between traction and separation is to establish a potential-based solution that defines the traction vector as a function of the interfacial potential gradient with respect to the components of the separation vector [35,36]. The interfacial potential is defined as a function of the components of the separation vector as

$$p=\nabla \phi \left({\delta }_n,{\delta }_{t_1},{\delta }_{t_2}\right)=\frac{\partial \phi }{\partial {\delta }_n}n+\frac{\partial \phi }{\partial {\delta }_{t_1}}{t}_1+\frac{\partial \phi }{\partial {\delta }_{t_2}}{t}_2$$

(4)

One can further assume that the interfacial potential is a function of an intermediate single variable $\delta$ which is a function of the separation vector, thus

$$p=\nabla \phi \left(\delta \right)=\left(\frac{\partial \delta }{\partial {\delta }_n}n+\frac{\partial \delta }{\partial {\delta }_{t_1}}{t}_1+\frac{\partial \delta }{\partial {\delta }_{t_2}}{t}_2\right)f\left(\delta \right)$$

(5)

where

$$\delta =\sqrt{{\delta }_n^2+\frac{{\beta }_1^2}{{\alpha }_1^2}{\delta }_{t_1}^2+\frac{{\beta }_2^2}{{\alpha }_2^2}{\delta }_{t_2}^2}$$

(6)

is defined as the weighted combination of the components of the separation vector, where the ${\alpha }_1$, ${\alpha }_2$, ${\beta }_1$, and ${\beta }_2$ are the material parameters which will be discussed later in this section. The derivatives of $\delta$ with respect to ${\delta }_n$, ${\delta }_{t_1}$, and ${\delta }_{t_2}$ are

$$\begin{array}{l}\frac{\partial \delta }{\partial {\delta }_n}=\frac{{\delta }_n}{\delta }\\ \frac{\partial \delta }{\partial {\delta }_{t_1}}=\frac{{\beta }_1^2}{{\alpha }_1^2}\frac{{\delta }_{t_1}}{\delta }\\ \frac{\partial \delta }{\partial {\delta }_{t_2}}=\frac{{\beta }_2^2}{{\alpha }_2^2}\frac{{\delta }_{t_2}}{\delta }\end{array}$$

(7)

While

$$f\left(\delta \right)=\frac{\mathrm{d}\phi \left(\delta \right)}{\mathrm{d}\delta }$$

(8)

defines the shape of the traction-separation curve and should be determined according to experimental data. Theoretically, $f\left(\delta \right)$ could be in any form, such as a linear function, bi-linear function, or a general nonlinear function, however different types of materials should use different traction-separation curves [41,42]. The discussion of the influence of the shape of the traction-separation curves on the numerical results is out of the scope of this work. The interested reader can find more information in [41,42].

The components of the traction vector are

$$\begin{array}{l}\sigma =\frac{{\delta }_n}{\delta }f\left(\delta \right)\\ {\tau }_1=\frac{{\beta }_1^2}{{\alpha }_1^2}\frac{{\delta }_{t_1}}{\delta }f\left(\delta \right)\\ {\tau }_2=\frac{{\beta }_2^2}{{\alpha }_2^2}\frac{{\delta }_{t_2}}{\delta }f\left(\delta \right)\end{array}$$

(9)

The energy release calculated from Equation (5) is

$$G={\displaystyle \int p\cdot \mathrm{d}\mathbf{\delta }}={\displaystyle \int \nabla \phi \left(\delta \right)\cdot \mathrm{d}\mathbf{\delta }}={\displaystyle \int f\left(\delta \right)\mathrm{d}\delta }$$

(10)

which indicates that the total energy release for a growing crack does not depend on the opening path, i.e., the fracture energy release rates calculated from Equation (5) for mode I and mode II are the same.

In order to introduce different fracture energy release rates for different modes, following the approach of Snozzi and Molinari [43], the tractions are modified to

$$\begin{array}{l}\sigma =\frac{{\delta }_n}{\delta }f\left(\delta \right)\\ {\tau }_1=\frac{{\beta }_1^2}{{\alpha }_1}\frac{{\delta }_{t_1}}{\delta }f\left(\delta \right)\\ {\tau }_2=\frac{{\beta }_2^2}{{\alpha }_2}\frac{{\delta }_{t_2}}{\delta }f\left(\delta \right)\end{array}$$

(11)

which needs four material parameters (${\alpha }_1$, ${\beta }_1$, ${\alpha }_2$, and ${\beta }_2$) and one function $f\left(\delta \right)$ to define the material’s traction-separation law. It is noted that the traction-separation law introduced in [43] can be retrieved from the proposed model when ${\alpha }_1={\alpha }_2$ and ${\beta }_1={\beta }_2$.

For pure tension, one has

$$\begin{array}{l}\delta ={\delta }_n\\ \sigma =f\left(\delta \right)\end{array}$$

(12)

which implies that the traction-separation curve can be determined through pure tension. Equation (12) yields

$$\begin{array}{l}{\delta }_{nc}={\delta }_c\\ {\sigma }_c=f\left(0\right)\\ G_{\mathrm{I}}={\displaystyle \int \sigma \mathrm{d}{\delta }_n}={\displaystyle \int f\left(\delta \right)\mathrm{d}\delta }\end{array}$$

(13)

where ${\delta }_c$ is a critical value at which point the function $f\left({\delta }_c\right)$ = 0, while ${\delta }_{nc}={\delta }_c$ is the maximum tensile separation at which point the tensile stress is $\sigma =0$. ${\sigma }_c$ is the tensile strength and $G_{\mathrm{I}}$ is the fracture energy release rate for mode I fracture.

For pure shear in direction $t_1$, one has

$$\begin{array}{l}\delta =\frac{{\beta }_1}{{\alpha }_1}{\delta }_{t_1}\\ {\tau }_1={\beta }_{1}f\left(\delta \right)\end{array}$$

(14)

which yields

$$\begin{array}{l}{\delta }_{t_{1}c}=\frac{{\alpha }_1}{{\beta }_1}{\delta }_c=\frac{{\alpha }_1}{{\beta }_1}{\delta }_{nc}\\ {\tau }_{1c}={\beta }_{1}f\left(0\right)={\beta }_1{\sigma }_c\\ G_{{\mathrm{II}}_1}={\displaystyle \int {\tau }_1\mathrm{d}{\delta }_{t_1}}={\alpha }_1{\displaystyle \int f\left(\delta \right)\mathrm{d}\delta }={\alpha }_1{G}_{\mathrm{I}}\end{array}$$

(15)

where ${\delta }_{t_{1}c}$ and ${\tau }_{1c}$ are the maximum shear separation and the shear strength in direction $t_1$, respectively, while $G_{{\mathrm{II}}_1}$ is the fracture energy release rate with respect to $t_1$. Equation (15) implies that the material parameter ${\alpha }_1$ is the ratio of $G_{{\mathrm{II}}_1}$ and $G_{\mathrm{I}}$, while the parameter ${\beta }_1$ is the ratio of shear strength ${\tau }_{1c}$ and tensile strength ${\sigma }_c$.

Similarly, for pure shear in direction $t_2$, one has

$$\begin{array}{l}\delta =\frac{{\beta }_2}{{\alpha }_2}{\delta }_{t_2}\\ {\tau }_2={\beta }_{2}f\left(\delta \right)\end{array}$$

(16)

which yields

$$\begin{array}{l}{\delta }_{t_{2}c}=\frac{{\alpha }_2}{{\beta }_2}{\delta }_c=\frac{{\alpha }_2}{{\beta }_2}{\delta }_{nc}\\ {\tau }_{2c}={\beta }_{2}f\left(0\right)={\beta }_2{\sigma }_c\\ G_{{\mathrm{II}}_2}={\displaystyle \int {\tau }_2\mathrm{d}{\delta }_{t_2}}={\alpha }_2{\displaystyle \int f\left(\delta \right)\mathrm{d}\delta }={\alpha }_2{G}_{\mathrm{I}}\end{array}$$

(17)

where ${\delta }_{t_{2}c}$ and ${\tau }_{2c}$ are the maximum shear separation and the shear strength in direction $t_2$, respectively, while $G_{{\mathrm{II}}_2}$ is the fracture energy release rate with respect to $t_2$. Equation (17) implies that the material parameter ${\alpha }_2$ is the ratio of $G_{{\mathrm{II}}_2}$ and $G_{\mathrm{I}}$, while the parameter ${\beta }_2$ is the ratio of shear strength ${\tau }_{2c}$ and tensile strength ${\sigma }_c$.

One can further introduce the damage as

$$d=\frac{\delta }{{\delta }_c}=\sqrt{\frac{{\delta }_n^2}{{\delta }_c^2}+\frac{{\beta }_1^2{\delta }_{t_1}^2}{{\alpha }_1^2{\delta }_c^2}+\frac{{\beta }_2^2{\delta }_{t_2}^2}{{\alpha }_2^2{\delta }_c^2}}$$

(18)

It is noted that the damage $d$ should be less than 1.0 since $d\ge 1.0$ means the material is completely damaged. To stabilize the damage representation equation, $d$ is therefore set to 1.0 as long as $d\ge 1.0$. Substituting Equations (13), (15), and (17) into Equation (18), one gets

$$d=\sqrt{d_n^2+d_{t_1}^2+d_{t_2}^2}=\sqrt{{\left(\frac{{\delta }_n}{{\delta }_{nc}}\right)}^2+{\left(\frac{{\delta }_{t_1}}{{\delta }_{t_{1}c}}\right)}^2+{\left(\frac{{\delta }_{t_2}}{{\delta }_{t_{2}c}}\right)}^2}$$

(19)

where

$$d_n=\frac{{\delta }_n}{{\delta }_{nc}};d_{t_1}=\frac{{\delta }_{t_1}}{{\delta }_{t_{1}c}};d_{t_2}=\frac{{\delta }_{t_2}}{{\delta }_{t_{2}c}}$$

(20)

are the damage components in $n$, $t_1$, and $t_2$, respectively.

The tractions can also be reformed as functions of damage

$$\begin{array}{l}\sigma =\frac{{\delta }_n}{\delta }f\left(\delta \right)=\frac{d_n}{d}{\sigma }_{c}z\left(d\right)\\ {\tau }_1=\frac{{\beta }_1^2}{{\alpha }_1}\frac{{\delta }_{t_1}}{\delta }f\left(\delta \right)=\frac{d_{t_1}}{d}{\tau }_{1c}z\left(d\right)\\ {\tau }_2=\frac{{\beta }_2^2}{{\alpha }_2}\frac{{\delta }_{t_2}}{\delta }f\left(\delta \right)=\frac{d_{t_2}}{d}{\tau }_{2c}z\left(d\right)\end{array}$$

(21)

where $z\left(d\right)$ is the normalized form of $f\left(\delta \right)$

$$z\left(d\right)=\frac{1}{{\sigma }_c}f\left(\delta \right)$$

(22)

In this paper, the heuristic softening function $z\left(d\right)$ introduced in [37] was utilized

$$z\left(d\right)=\left[1-\frac{a+b-1}{a+b}{e}^{d\left(a+cb/\left(\left(a+b\right)\left(1-a-b\right)\right)\right)}\right]\times \left[a\left(1-d\right)+b{\left(1-d\right)}^c\right]$$

(23)

where $a$, $b$, $c$ are the parameters chosen to fit a particular experimental curve.

To generalize the traction-separation model proposed in this work, Equation (21) is redefined as

$$\begin{array}{l}\sigma ={\sigma }_c{g}_n\left(d_n,d\right)z\left(d\right)\\ {\tau }_1={\tau }_{1c}{g}_{t_1}\left(d_{t_1},d\right)z\left(d\right)\\ {\tau }_2={\tau }_{2c}{g}_{t_2}\left(d_{t_2},d\right)z\left(d\right)\end{array}$$

(24)

where $g_n$ is a function of $d_n$ and $d$; $g_{t_1}$ is a function of $d_{t_1}$ and $d$; $g_{t_2}$ is a function of $d_{t_2}$ and $d$. In a simplified form, one can further assume that $g_n=g_{t_1}=g_{t_2}=1$, which yields

$$\begin{array}{l}\sigma ={\sigma }_{c}z\left(d\right)\\ {\tau }_1={\tau }_{1c}z\left(d\right)\\ {\tau }_2={\tau }_{2c}z\left(d\right)\end{array}$$

(25)

It is worth noting that the traction-separation law implemented in the 3D combined single smeared crack model [14,39] can be retrieved from Equation (25) when ${\delta }_{t_{1}c}={\delta }_{t_{2}c}$ and ${\tau }_{1c}={\tau }_{2c}$.

Since the users have more freedom to choose different forms for the functions $g_n$, $g_{t_1}$, $g_{t_2}$, and $z$, the traction-separation model defined through the Equations (19), (20), and (24) is referred to as a “generalized traction-separation law”. It is noted that, beside functions $g_n$, $g_{t_1}$, $g_{t_2}$, and $z$, the generalized traction-separation law needs six parameters, i.e., three critical separations ${\delta }_{nc}$, ${\delta }_{t_{1}c}$, ${\delta }_{t_{2}c}$, and three strengths ${\sigma }_c$, ${\tau }_{1c}$, ${\tau }_{2c}$. In addition, according to the derivation, one can conclude that the traction-separation law implemented in the combined single smeared crack model [37,38,39] and the traction-separation law defined through Equations (6) and (11) are two special cases of the proposed generalized traction-separation model.

4. Impact Fracture Process of Laminated Glass

Zang et al. experimentally studied the fracture generation and propagation of automotive glass under impact conditions [44]. In the experiments, a customized glass specimen, which had a thickness four times bigger than that of the general automotive glass (Figure 2), was used to capture the evolution of the fracture processes. The glass specimen was hit by an impactor at the mid-side. The fracture initialization and propagation processes near by the impact point were recorded by a photo-elastic device [44].

Photographic evidence of the fracture processes’ time sequences for two of the tests is shown in Figure 3. The fracture propagation and the stress distribution inside of the glass specimen (indicated by the fringes in the image) are shown in Figure 3a, where the time interval between images is 20 μs. Due to the field of view limit of the high speed camera, only the fracture process on the impact side of glass is captured [44]. Figure 3b, where the time interval between images is 100 μs, shows the complete fracture propagation process through both layers of glass [44].

The images presented in Figure 3 indicate that the failure phenomena are the same for both cases although there are some differences regarding the time at which the fracture processes occur. At the early stages of the impact process only the upper glass layer withstands the load and bends until fractures occur, while the stress levels on the lower glass layer are small. As time progresses, the lower glass layer starts to bend while the PVB layer is being compressed. This causes the impact load to reach the lower glass, which eventually fractures when the PVB layer is fully compressed. It is worth noting that this phenomenon will be used to qualitatively verify the numerical results.

The Y-code was used to simulate the impact fracture experiment described above, and the corresponding numerical model is shown in Figure 4. In this model, the bar of laminated glass is placed between four supports. The boundary conditions are such that the top surface of upper support and the bottom surface of lower support are fixed in space during the simulation. The weight and the initial velocity of impactor are 1 kg and 3.13 m/s respectively, which are the same as that in the experiment. The FDEM model includes 28,809 linear tetrahedron elements, which were generated using unstructured algorithms implemented in Gmsh [45].

Due to the relatively low impact velocity, the material models used to describe the deformation of the finite elements are linear elastic in nature. Fracture propagation is only allowed on the glass, while the rest of the components, i.e., impactor, support, and PVB layer, are considered to be continuum media. The material properties used in the simulations are shown in

Table 1. Material properties.

Parameter	Glass	PVB	Support	Impactor
Young’s module/GPa	75.0	0.5	5.0	200.0
Poisson’s ratio	0.2	0.4	0.4	0.27
Density/(kg/m3)	2456.0	100.0	2400.0	-
Tensile strength/MPa	100.0	-	-	-
Energy release rate/(N/m)	10.0	-	-	-

.

The numerical results for the impact fracture process on the laminated glass model are shown in Figure 5. In the figure all the pictures are front views, in accordance with the observation angle in the experimental work. A visualization perspective technology was used in the post processor to present the transparency of the glass parts, which enables the observation of the fracture states in the interior of glass layers. The black lines in the pictures indicate the cracks seen directly from the surface of the glass through the observation angle, while the shaded areas show the fracture surfaces in the interior of the glass layers.

The simulation results show that the first cracks occurs at 23 μs starting at the interface between the upper glass layer and the PVB layer right below the impact point (as shown in Figure 5a). As time progresses these initial cracks propagate toward the top surface of the upper glass layer in a very short period of time (about 12 μs). At 35 μs, the cracks have already penetrated the whole upper glass layer (see Figure 5b), however, the lower glass layer has no obvious cracks at this time and it takes a relatively long time for that to occur, i.e., 200 μs, as shown in Figure 5c. This occurs because the speed of propagation of the stress wave in the PVB layer is relatively slow, which makes that the load takes a longer time to reach the lower glass layer. The failure of the bottom glass layer starts at its free surface when the tensile stress exceeds the tensile strength of the glass material. Once this occurs, the fractures penetrate the whole glass layer at a very fast speed, as is the case for the upper glass layer (as shown in Figure 5d,e).

After comparing the numerical results to the experimental observations shown in Figure 3, it can be determined that the positions and the sequencing of cracks obtained from the model are in agreement with the test results. This strong evidence further confirms that the combined finite-discrete element method can be used to study fracture processes on glass assemblies.

It should be noted that the main purpose of this example is to qualitatively verify the numerical method. Therefore, the proper calibration of the material properties (and the correspondent sensitivity analysis) and quantitative analysis for specific applications are outside of the scope of this work. The numerical tests shown in this section demonstrate that, although the material parameters have an influence on the timing of the different cracks processes, the fracture propagation sequence and their relevant positions in the majority of the simulation results are in good agreement with the experimental observations, which proves the validity of the numerical method.

5. Fracture Patterns of Glass Plate under Impact

Researchers from various disciplines have conducted experimental investigations for glass plate impact fracture mechanisms. According to the experimental observations, it can be concluded that there are two main types of cracks observed when a glass plate is impacted by an object, namely radial cracks and circular cracks [46,47]. The radial cracks are those spreading radially outwards from the impact point, while the circular cracks are cyclic cracks centered around the impact point [46], as shown in Figure 6. Research has shown that both radial and circular cracks are generated by tensile stresses, i.e., they are tensile cracks. However, while circular cracks occur mainly on the surface of the glass that is facing the impactor, radial cracks occur mainly on the opposite surface of the glass (i.e., opposite to the impact side) [47].

In this work, the impact fracture processes on a glass plate is simulated by using Y-code. As shown in Figure 7, a glass plate with both ends constrained by two sets of supports is impacted by a cylindrical impactor with a mass of 50 g at the centre of its top surface. The size of the glass plate is 300 mm × 50 mm × 4.76 mm; and the radius of the impactor is 5 mm. The glass material properties used in this model are the same as those described in the last section. The supports and the impactor are considered to behave elastically and no fracture propagation is allowed in those parts. The Young’s modulus and the Poisson’s ratio for the supports and the impactor are 750 GPa and 0.2 respectively. In the FDEM model, the domain was discretized into 231,472 linear tetrahedron elements using unstructured algorithms implemented in Gmsh [45].

The sequence of the glass plate fracture propagation is shown in Figure 8. Here, the perspective technology is used again to present the transparency of glass materials, which enables to observe the fracture states in the interior of glass. As can be seen in Figure 8, both the radial and circular cracks are captured by the numerical model. Furthermore, radial cracks occur mainly on the opposite surface (i.e., the backside of the impactor), and circular cracks are observed on the impact side (Figure 8d). This phenomenological result is in very good agreement with the experimental observations reported by Bennett [46] and Bertino [47] et al.

6. Conclusions

In this work, we used 3D FDEM to simulate impact fracture processes on a laminated glass beam and on a monolithic glass plate. The numerical results show enriched fracture details, especially the circular and radial cracks in the monolithic glass plate that agree well with the corresponding experimental results. The main contribution of this work is the presentation of a generalized traction-separation model that further advances the fracture and fragmentation methodologies for FDEM. The simulation work presented utilized the open source Y-code which is based on the first generation FDEM algorithms. Given the recent code and method advancements in the field, it is expected these results, given due consideration, will offer automobile industry safety analysts a new opportunity for safety exploration.

This work demonstrated that 3D FDEM is an effective approach for modeling fracture and fragmentation of glass under dynamic loading. However, the 2.5D FDEM solution for fracture and fragmentation of thin shells and plates, as proposed by the authors [14,15], is a much more effective approach when shells made of glass (e.g., glass windows) need analysis.