Simulating progressive failure in laminated glass beams with a layer-wise randomized phase-field solver

Abstract

Laminated glass achieves improved post-critical response through the composite effect of stiff glass layers and more compliant polymer films, manifested in progressive layer failure by multiple localized cracks. As a result, laminated glass exhibits greater ductility than non-laminated glass, making structures made with it suitable for safety-critical applications while maintaining their aesthetic qualities. However, such post-critical response is challenging to reproduce using deterministic failure models, which mostly predict failure through a single through-thickness crack localized simultaneously in all layers. This numerical-experimental study explores the extent to which progressive failure can be predicted by a simple randomized model, where layer-wise tensile strength is modeled by independent, identically distributed Weibull variables. On the numerical side, we employ a computationally efficient, dimensionally-reduced phase field formulation – with each layer considered to be a Timoshenko beam – to study progressive failure through combinatorial analysis and detailed Monte Carlo simulations. The reference experimental data were obtained from displacement-controlled four-point bending tests performed on multi-layer laminated glass beams. For certain combinations of the glass layer strengths, results show that the randomized model can reproduce progressive structural failure and the formation of multiple localized cracks in the glass layers. However, the predicted response was less ductile than that observed in experiments, and the model could not reproduce the most frequent glass layer failure sequence. These findings highlight the need to consider strength variability along the length of a beam and to include it in phase-field formulations.

1 Introduction

Laminated glass combines stiff glass layers with compliant, transparent polymer interlayers. The interlayers decrease structural stiffness in the pre-peak range but greatly enhance ductility because of the gradual failure of the brittle glass layers. This composite effect enables the deployment of laminated glass in architectural, structural, or automotive engineering for safety-critical applications (see, e.g., Haldimann et al. (2008); Bedon et al. (2018); Jozwik (2022)) while simultaneously making the mechanics of laminated glass structures intriguing.

The prevailing design approach to modeling the pre-fracture behavior of laminated glass is to replace a heterogeneous structure with a homogeneous one, with the overall thickness adjusted to achieve an equivalent response. This ”effective thickness” concept has been successfully implemented, e.g., for deflection Galuppi and Royer-Carfagni (2012a, b), buckling López-Aenlle et al. (2016); D’Ambrosio and Galuppi (2020), and free vibration Aenlle and Pelayo (2013); Zemanová et al. (2018) analyses of laminated glass beams and plates. Several analytical studies have also addressed the estimation of failure loads of laminated glass beams Foraboschi (2007), optimal design under mechanical and aesthetic criteria Foraboschi (2014), or the post-breakage stiffness of laminated glass structures under in-plane Galuppi and Royer-Carfagni (2016) and out-of-plane Galuppi and Royer-Carfagni (2018) loads. Although analytical approaches provide invaluable insight into the mechanics of laminated glass structures, the detailed post-fracture response has mostly been investigated with numerical means, including discrete element Zang et al. (2007); Baraldi et al. (2016), cohesive zone Chen et al. (2016); Vocialta et al. (2018), and peridynamics Wu et al. (2020); Naumenko et al. (2022) formulations (see also selected reviews Wang et al. (2017); Teotia and Soni (2018); Kuntsche et al. (2019) for a detailed comparison).

In this contribution, we focus on the computational modeling of laminated glass structures using dimensionally-reduced phase-field fracture models. The phase-field approach was first introduced by Bourdin et al. Bourdin et al. (2000) in 2000 as a continuous approximation to the free discontinuity problem arising from Francfort-Marigo’s variational reformulation Francfort and Marigo (1998) of Griffith’s fracture evolution theory Griffith (1921). Since then, this strategy has gained significant appeal in the computational mechanics community because (i) no apriori assumptions are made on the crack path, (ii) crack propagation criteria follow directly from the optimality conditions that replace crack tracking algorithms, and (iii) the formulation is objective under mesh refinement. Interested readers are referred to Bourdin et al. (2008); Ambati et al. (2015); Wu et al. (2020) for detailed overviews.

The original phase-field formulations were first developed for two- and three-dimensional bodies, and their adaptation to dimensionally-reduced models applicable to beams, plates, or shells appeared considerably later. Amiri et al.’s Amiri et al. (2014) pioneering formulation addressed Kirchhoff-Love thin shells with through-thickness cracks and employed maximum-entropy mesh-free approximation to discretize the governing equations. Kiendl et al. Kiendl et al. (2016) extended this work to account for the tension-compression asymmetry in the thickness direction and based discretization on isogeometric analysis. Further refinements focused on developing formulations for Reissner-Mindlin shells Kikis et al. (2021). A more accurate representation of partial fracture along thickness was achieved by using a specific ansatz for the phase field in Euler-Bernoulli beams Lai et al. (2020) or by employing a separate, dimensionally-reduced representation for kinematic variables and a full representation for the phase field that are coupled at quadrature points Ambati et al. (2022). For an illustration of the predictive power of dimensionally-reduced fracture models, see, e.g., a recent study Bijaya et al. (2023) on fracture propagation in architectured beam lattices.

Several earlier studies utilized a phase-field approach to simulate the post-fracture response of laminated glass structures. Alessi and Freddi obtained the first results for hybrid glass laminates, considering them as elastic-brittle glass layers connected with cohesive interfaces, employing one-dimensional Alessi and Freddi (2017) and two-dimensional Alessi and Freddi (2019) models. Freddi and Mingazzi extended this approach to tall laminated glass beams Freddi and Mingazzi (2020). In our two previous works, we used layer-wise beam and plate models to simulate the response of multi-layered glass samples under quasi-static Schmidt et al. (2020) and low-velocity impact Schmidt et al. (2023) loads.

This article extends our quasi-static study Schmidt et al. (2020) that employed several variants of phase-field formulations combined with 2D plane stress, a layer-wise Reissner-Mindlin plate, and Timoshenko beam structural models to simulate failure in laminated glass structures. The study confirmed that the layer-wise beam theories represent the most suitable options in terms of balancing accuracy and computational requirements. However, all the formulations considered failed to replicate the distributed cracking and the gradual failure sequence observed in the experimental campaign that complemented the numerical investigations.

Novelty and scope

This work investigates the gradual and distributed cracking in laminated glass structures with randomized phase field models. A well-established Weibull approach Weibull (1951) is adopted to model the intrinsic variability of the glass strength (see, e.g., a recent review for details) and combined with the Monte Carlo simulations (see, e.g., related studies Bonati et al. (2019, 2020) on probabilistic modeling based on the Weibull statistics and Biolzi et al. (2017); Casolo and Diana (2018) on numerical approaches based on the discrete element method or peridynamics). The dimensionally-reduced quasi-static layer-wise beam formulation Schmidt et al. (2020) is used as the deterministic solver to keep the simulation cost manageable. Note that our previous results Schmidt et al. (2023) show that a dynamic phase field formulation can reproduce crack branching in laminated glass under impact loads. However, its computational demands are prohibitive for stochastic simulations and a similar argument holds for the experimental campaign. Besides, because the first steps to the mesh-independent stochastic phase-field fracture approaches have appeared only recently Hai et al. (2023); Wu et al. (2023), the simplest stochastic model is adopted, which accounts for the variability of the tensile strength in the thickness direction by assuming that layer-wise strengths can be described by independent identically-distributed Weibull random variables.

To this end, the remainder of the paper is organized as follows. Section 2 introduces the deterministic solver used to predict the response of laminated glass beams, valorizing the main insights from Schmidt et al. (2020). Section 3 provides details about the experimental campaign, material parameters for interlayer and glass layers, the setting of the deterministic solver, and its verification using a benchmark problem. Section 4 collects simulation results, namely a combinatorial analysis to explore the range of responses attainable by the model. It also provides insights into the experimental data inferred from the Monte Carlo simulations. Finally, Section 5 summarizes the most important findings.

Throughout the text, standard notation is employed, using $a$, $\boldsymbol{a}$, and $\boldsymbol{A}$ to denote a scalar quantity, a vector or a second-order symmetric tensor, and a symmetric fourth-order tensor, respectively. The single and double dots indicate the corresponding contractions, e.g., $\boldsymbol{a}\cdot\boldsymbol{b}$ or $\boldsymbol{A}:\boldsymbol{b}$, and $\boldsymbol{\nabla}$ and ${}^{\prime}$ represent the gradient and derivative operators. Other symbols are introduced as needed.

2 Deterministic phase-field fracture solver

The current section outlines the development of the computational tool used in fracture simulations, adapting the results of our previous study Schmidt et al. (2020) and adopting the following modeling assumptions:

1. the sample is subjected to quasi-static loading,

2. fracture occurs only in the glass layers,

3. viscoelastic response of the polymer interlayer is approximated by the quasi-elastic model with time-dependent stiffness, e.g., Galuppi and Royer-Carfagni (2012, 2013); Zemanová et al. (2017),

4. perfect bonding is maintained between glass and polymer layers, e.g., (Chen et al., 2017, Section 2.3),

5. each layer is modeled as a Timoshenko beam Timoshenko (1921, 1922) with through-thickness cracks.

2.1 Model setup

We consider the domain of a multi-layer laminated glass sample $\Omega$ to be divided into $M$ layers $\Omega_{m}$. Assuming the common setup with an odd number of layers, the odd indices $\{1,3,\ldots,M\}$ correspond to the glass layers, whereas the even indices $\{2,4,\ldots,M-1\}$ to the polymer interlayers; see Figure 1 for an illustration.

Figure 1: Subdomains in a multi-layer laminated glass beam sample occupying domain $\Omega$. The structure consists of $M$ layers $\Omega_{m}$, with the odd indices ($m=1,3,\ldots,M$) corresponding to the glass layers and even indices ($m=2,4,\ldots,M-1$) to the polymer interlayers.

Under the assumptions (A1)–(A4), the mechanical response of the laminated structure is encoded in the energy functional

$$\mathcal{E}(t,\boldsymbol{u},\boldsymbol{d})=\sum_{m=2,4}^{M-1}\mathcal{E}_{m}(t,\boldsymbol{u}_{m})+\sum_{m=1,3}^{M}\mathcal{E}_{m}(t,\boldsymbol{u}_{m},d_{m}),$$ (1)

where $t$ denotes the instant time varying over the time interval $[0,T]$, $\boldsymbol{u}=(\boldsymbol{u}_{m})_{m=1,2,\ldots,M}$ the displacement fields in all layers, and $\boldsymbol{d}=(\boldsymbol{d}_{m})_{m=1,3,\ldots,M}$ the damage phase fields in the glass layers (the values of $d_{m}=0$ and $d_{m}=1$ correspond to the intact and fully damage states, respectively). The layer-level energy functionals $\mathcal{E}_{m}$ have a different structure for glass and polymer layers as specified next.

For polymer layers, the assumptions (A1)–(A4) imply that their energies follow from

$$\mathcal{E}_{m}(t,\boldsymbol{u}_{m})=\Psi^{\text{e}}_{m}(t,\boldsymbol{u}_{m})-\mathcal{P}_{m}(t,\boldsymbol{u}_{m}), m=2,4,\ldots,M-1,$$ (2)

where $\mathcal{P}_{m}$ stands for the work done by external forces (if present). The elastically stored energy $\Psi^{\text{e}}_{m}$ is given as

$$\Psi^{\text{e}}_{m}(t,\boldsymbol{u}_{m})=\int_{\Omega_{m}}\psi^{\text{e}}_{m}\left(t,\boldsymbol{\varepsilon}_{m}\right)\,{\mathrm{d}}V, m=2,4,\ldots,M-1,$$ (3)

where $\boldsymbol{\varepsilon}_{m}=(\boldsymbol{\nabla}\boldsymbol{u}_{m}+(\boldsymbol{\nabla}\boldsymbol{u}_{m})^{\mathsf{T}})/2$ designates the small strain tensor and the elastic strain energy density $\psi^{\text{e}}_{m}$ has the standard quadratic form

$$\psi^{\text{e}}_{m}(t,\boldsymbol{\varepsilon}_{m})=\textstyle{\frac{1}{2}}\,\boldsymbol{\varepsilon}_{m}:\boldsymbol{C}_{m}(G_{m}(t),\nu_{m}):\boldsymbol{\varepsilon}_{m}, G_{m}(t)=G_{\text{Maxwell}}\left(\frac{t}{2},\theta\right),$$ (4)

with $\boldsymbol{\varepsilon}=\boldsymbol{\nabla}_{\text{s}}\boldsymbol{u}$ denoting the strain tensor and $\boldsymbol{C}_{m}$ the isotropic stiffness tensor determined from the time-dependent quasi-elastic shear modulus $G_{m}$ and time-independent Poisson ratio $\nu_{m}$, see Duser et al. (1999); Zemanová et al. (2017) for further discussion and justification. We assume that the shear modulus can be accurately predicted using the temperature-dependent generalized Maxwell chain model with parameters identified in Hána et al. (2019) are evaluated in the half of the current simulation time according to Duser et al. (1999). Because the temperature $\theta$ is considered constant during the whole experiment, we omit the dependence of $G_{m}$ on this parameter for notational simplicity.

Assumptions (A1) and (A2) imply that energies of the glass layers attain the form

$$\mathcal{E}_{m}(t,\boldsymbol{u}_{m},d_{m})=\Psi^{\text{e}}_{m}(\boldsymbol{u}_{m},d_{m})+\Psi^{\text{d}}_{m}(d_{m})-\mathcal{P}_{m}(t,\boldsymbol{u}_{m}), m=1,3,\ldots,M,$$ (5)

that incorporates the damage variable $d_{m}$ to account for dissipative processes. In particular, the elastically stored energy in the following form Miehe et al. (2010); Pham et al. (2011) is considered:

$$\bar{\Psi}^{\text{e}}_{m}(\boldsymbol{u}_{m},d_{m})=\int\limits_{\Omega_{m}}\Bigl{(}(1-d_{m})^{2}\,\psi^{\textrm{e}}_{+,m}(\boldsymbol{\varepsilon}_{m})+\psi^{\textrm{e}}_{-,m}(\boldsymbol{\varepsilon}_{m})\Bigr{)}\,{\mathrm{d}}V, m=1,3,\ldots,M,$$ (6)

where the split of the energy density into the positive, $\psi^{\textrm{e}}_{+,m}$, and negative, $\psi^{\textrm{e}}_{-,m}$, parts is discussed later. Finally, the dissipated energy in the $m$-th layer $\Psi^{\text{d}}_{m}$ is considered in the regularized form proposed in Pham et al. (2011)

$$\Psi^{\text{d}}_{m}(d_{m})={\textstyle\frac{3}{8}}\,G^{\text{c}}_{m}\int\limits_{\Omega_{m}}\Bigl{(}\frac{d_{m}}{\ell_{m}}+\ell_{m}\|\boldsymbol{\nabla}d_{m}\|^{2}\Bigr{)}\,{\mathrm{d}}V, m=1,3,\ldots,M,$$ (7)

where $G^{\text{c}}_{m}$ stands for the fracture energy of the $m$-th layer and the regularization parameter $\ell_{m}$ sets the lengthscale of the damage profile.

To track the evolution of the structure over a time interval $[0,T]$, the time interval is discretized into $N$ time instants $0=t_{0}<t_{1}<t_{2}<\cdots<t_{N}=T$. . Then, given the initial data $(\boldsymbol{u}_{0},\boldsymbol{d}_{0})=((\boldsymbol{u}(t_{0}),\boldsymbol{d}(t_{0}))$, the displacement and damage fields $(\boldsymbol{u}_{n},\boldsymbol{d}_{n})=((\boldsymbol{u}(t_{n}),\boldsymbol{d}(t_{n}))$ at time $t_{n}$ follow from the recursive energy minimization

$$\mathcal{E}(t_{n},\boldsymbol{u}_{n},\boldsymbol{d}_{n})\leq\mathcal{E}(t_{n},\boldsymbol{u},\boldsymbol{d})\text{ for all admissible pairs }(\boldsymbol{u},\boldsymbol{d}), n=1,2,\ldots,N.$$ (8)

By the assumption (A4), the admissible displacement field $\boldsymbol{u}$ must remain continuous in $\Omega$ (and satisfy the kinematic boundary conditions on the boundary $\partial\Omega$ that we did not discuss for the sake of brevity). The admissible damage fields are subjected to the box constraints

$$d_{n-1,m}\leq d_{n,m}\leq 1, n=1,2,\ldots,N;\quad m=1,3,\ldots,M,$$ (9)

where the first inequality accounts for the damage irreversibility, i.e., the damage variable cannot decrease in time.

2.2 Dimensional reduction

Assumption (A5) implies the following distribution of normal, $\varepsilon_{m}$, and shear, $\gamma_{m}$, strains in the $m$-th layer:

$$\varepsilon_{m}(x,z_{m}) =u^{\prime}_{m}(x)+z_{m}\varphi_{m}^{\prime}(x)\,=\varepsilon_{\text{c},m}(x)+z_{m}\kappa_{m}(x),$$ (10a)

$$\gamma_{m}(x) =\varphi_{m}(x)+w_{m}^{\prime}(x), m=1,\ldots,M,$$ (10b)

where $x\in(0,L)$ and $z_{m}\in(-h_{m}/2,h_{m}/2)$ denote the position in the layer local coordinate system, $L$ and $h_{m}$ the beam length and the layer thickness, and $u_{m}$, $w_{m}$, and $\varphi_{m}$ the cross-section horizontal displacement, vertical displacement, and rotation; see also Figure 2 for illustration. In addition, $\varepsilon_{\text{c},m}$ designates the centerline strain and $\kappa_{m}$ the centerline pseudo-curvature. For later reference, we denote the layer (generalized) displacement fields by $\boldsymbol{u}_{m}=(u_{m},w_{m},\varphi_{m})$.

Figure 2: Local coordinate systems, geometrical parameters, and primary unknowns of the layer-wise beam model (left) and through-thickness distribution of normal strain in the $m$-th layer and the adopted numerical quadrature (right).

The rest of the model development follows the same steps as in Section 2.1. For polymer interlayers, the elastically stored energy (3) transforms into

$$\Psi^{\text{e}}_{m}(t,\boldsymbol{u}_{m})=\int_{0}^{L}\psi^{\text{e}}_{m}(t,\varepsilon_{\text{c},m},\kappa_{m},\gamma_{m})\,{\mathrm{d}}x, m=2,4,\ldots,M-1,$$ (11)

with the cross-section elastic energy

$$\psi^{\text{e}}_{m}(t,\varepsilon_{\text{c},m},\kappa_{m},\gamma_{m})={\textstyle\frac{1}{2}}\,G_{m}(t)\bigl{(}2(1+\nu_{m})\left(A_{m}\varepsilon_{\text{c},m}^{2}+I_{m}\kappa_{m}^{2}\right)+A^{*}_{m}\gamma_{m}^{2}\bigr{)}$$ (12)

for $m=2,4,\ldots,M-1$. Equation (12) additionally involves the layer cross-section characteristics

$$A_{m}=bh_{m}, I_{m}=\frac{1}{12}bh^{3}_{m}, A^{*}_{m}=\frac{5}{6}bh_{m},$$ (13)

where $A_{m}$ stands for the cross-section area, $I_{m}$ for the second moment of area, and $A^{*}_{m}$ for the effective shear area.

For glass layers ($m=1,3,\ldots,M$), the elastically stored energy (6) takes the form

$$\Psi^{\text{e}}_{m}(\boldsymbol{u}_{m},d_{m})=\int_{0}^{L}\Bigl{(}(1-d_{m})^{2}\psi^{\text{e}}_{+,m}(\varepsilon_{\text{c},m},\kappa_{m},\gamma_{m})+\psi^{\text{e}}_{-,m}(\varepsilon_{\text{c},m},\kappa_{m},\gamma_{m})\Bigr{)}\,{\mathrm{d}}x,$$ (14)

with the cross-section tensile and compressive energies following from a numerical quadrature in the $z_{m}$ variable, as first proposed in Kiendl et al. (2016):

$$\psi^{\text{e}}_{+,m}(\varepsilon_{\text{c},m},\kappa_{m},\gamma_{m}) ={\textstyle\frac{1}{2}}\,E_{m}b\,\sum_{j=1}^{J}\left\langle\varepsilon_{\text{c},m}+\kappa_{m}\,z_{m,j}\right\rangle^{2}_{+}\Delta z_{m}+{\textstyle\frac{1}{2}}\,G_{m}A^{*}_{m}\gamma_{m}^{2},$$ (15a)

$$\psi^{\text{e}}_{-,m}(\varepsilon_{\text{c},m},\kappa_{m},\gamma_{m}) ={\textstyle\frac{1}{2}}\,E_{m}b\,\sum_{j=1}^{J}\left\langle\varepsilon_{\text{c},m}+\kappa_{m}\,z_{m,j}\right\rangle^{2}_{-}\Delta z_{m},\quad m=1,3,\ldots,M.$$ (15b)

Here, $E_{m}$ stands for the Young modulus of the $m$-th layer, $J$ for the number of integration points located at $z_{m,j}$, $\Delta z_{m}=h_{m}/J$ their integration weight (see Figure 2), and $\langle\bullet\rangle_{+}$ and $\langle\bullet\rangle_{-}$ for the positive and negative parts of $\bullet$. The regularized dissipated energy, on the other hand, follows directly from the original expression (7):

$$\Psi^{\text{d}}_{m}(d_{m})={\textstyle\frac{3}{8}}\,G^{\text{c}}_{m}A_{m}\int_{0}^{L}\Bigl{(}\frac{d_{m}}{\ell_{m}}+\ell_{m}\left(d_{m}^{\prime}\right)^{2}\Bigr{)}\,{\mathrm{d}}x, m=1,3,\ldots,M$$ (16)

because the damage variable is assumed to be constant over a cross-section; recall assumption (A5).

Having reformulated all the components of the model introduced in Section 2.1, the evolution of the system follows from the same incremental energy minimization problem (8) under damage irreversibility enforced by (9). The interlayer displacement continuity implies that

$$w_{m} =w_{1}, m=2,3,\ldots,M,$$ (17a)

$$u_{m-1}+{\textstyle\frac{1}{2}}\varphi_{m-1}h_{m-1} =u_{m}-{\textstyle\frac{1}{2}}\varphi_{m}h_{m}, m=2,3,\ldots,M.$$ (17b)

Therefore, the centerline horizontal displacements $u_{m}$ and cross-section rotations $\varphi_{m}$ in the polymer layers can be expressed using the respective quantities in the glass layers:

$$u_{m} ={\textstyle\frac{1}{4}}\,\left(h_{m-1}\varphi_{m-1}-h_{m+1}\varphi_{m+1}\right)+{\textstyle\frac{1}{2}}\,\left(u_{m-1}+u_{m+1}\right), m=2,4,\ldots,M-1,$$ (18a)

$$\varphi_{m} =\frac{1}{h_{m}}\left(u_{m+1}-u_{m-1}-{\textstyle\frac{1}{2}}\,\varphi_{m-1}h_{m-1}-{\textstyle\frac{1}{2}}\,\varphi_{m+1}h_{m+1}\right), m=2,4,\ldots,M-1.$$ (18b)

For instance, the kinematics of a $5$-layer laminated glass beams ($3$ glass layers and $2$ interlayers) is specified by $7$ independent fields: $w_{1},u_{1},\varphi_{1},u_{3},\varphi_{3},u_{5}$, and $\varphi_{5}$.

2.3 Weak form of governing equations

The equations governing the evolution of displacement and damage fields in the layer-wise beam model follow from the optimality conditions of the variational problem (8). In particular, the optimality conditions with respect to the kinematic quantities yield:

$$\boldsymbol{\nabla}_{\boldsymbol{u}_{m}}\Psi^{\text{e}}_{m}(\boldsymbol{u}_{n,m},d_{n,m})\cdot\delta\boldsymbol{u}_{m}-\boldsymbol{\nabla}_{\boldsymbol{u}_{m}}\mathcal{P}_{m}(t_{n},\boldsymbol{u}_{n,m}) =0, m=1,3,\ldots,M,$$ (19a)

$$\boldsymbol{\nabla}_{\boldsymbol{u}_{m}}\Psi^{\text{e}}_{m}(t_{n},\boldsymbol{u}_{n,m})\cdot\delta\boldsymbol{u}_{m}-\boldsymbol{\nabla}_{\boldsymbol{u}_{m}}\mathcal{P}_{m}(t_{n},\boldsymbol{u}_{n,m}) =0, m=2,4,\ldots,M-1,$$ (19b)

where both the solution $\boldsymbol{u}_{n,m}$ at time $t_{n}$ and an admissible variation $\delta\boldsymbol{u}_{m}$ satisfy the continuity conditions (2.2). For reader’s convenience, the full form of the derivatives $\boldsymbol{\nabla}_{\boldsymbol{u}_{m}}\Psi^{\text{e}}_{m}$ is collected in A.

Because of the irreversibility constraints on the damage variables (9), the respective optimality conditions attain the form of a variational inequality

$$2\int_{0}^{L}\left(d_{n,m}-1\right)Y_{n,m}\delta{d}_{m}\,{\mathrm{d}}x+{\textstyle\frac{3}{8}}\,G^{\text{c}}_{m}A_{m}\int_{0}^{L}\Bigl{(}\frac{\delta d_{m}}{\ell_{m}}+2\ell_{m}d^{\prime}_{n,m}\delta d^{\prime}_{m}\Bigr{)}\,{\mathrm{d}}x\geq 0,$$ (20)

where $m=1,3,\ldots,M$ and the admissible variations $\delta d_{m}$ satisfy $d_{n-1,m}\leq d_{n,m}+\delta d_{m}\leq 1$, cf. Jirásek and Zeman (2015). The quantity $Y_{n,m}$ abbreviates the damage driving force obtained as

$$Y_{n,m}(x)=\psi^{\text{e}}_{+,m}\left(\varepsilon_{\text{c},m}(x),\kappa_{m}(x),\gamma_{m}(x)\right), x\in(0,L);$$ (21)

recall (15a) for the definition of the tensile cross-section energy $\psi^{\text{e}}_{+,m}$.

Our previous study Schmidt et al. (2020) revealed that the most efficient and robust phase field solver is obtained when replacing the variationally consistent driving force in Eq. (20) with

$$Y_{n,m}(x)\approx\textstyle{\frac{1}{2}}\,E_{m}A_{m}\max\bigl{(}\langle\varepsilon_{m}(x,h_{m}/2)\rangle_{+}^{2},\langle\varepsilon_{m}(x,-h_{m}/2)\rangle_{+}^{2}\bigr{)}, x\in(0,L),$$ (22)

where the normal strain in the $m$-th layer $\varepsilon_{m}$ follows from the displacement field $\boldsymbol{u}_{n,m}$ using the parameterization (10a) . This choice effectively introduces a Rankine-type failure condition, as first proposed in Miehe et al. (2015); see also Bilgen and Weinberg (2019) and Section 3.3 for further discussion. Notice that the modified version of the driving force $Y_{n,m}$ enters only the damage evolution equations (20), whereas the original form of the tensile energy $\psi^{\text{e}}_{+,m}$ from Eq. (15a) is used in the equilibrium equations (2.3).

2.4 Implementation

The set of two governing equations (2.3) and (20) is resolved with the Finite Element Method (FEM), employing a Python interface to the FEniCS library Logg et al. (2012) that allows an automatic treatment of the corresponding weak forms. This section briefly outlines the implementation strategy; the full details are available in the accompanying GitLab repository Schmidt (2023).

In particular, we employ the staggered approach to resolve the two optimality conditions sequentially. Following Bourdin et al. (2000), the kinematic quantities are first updated from (2.3) with the damage fields fixed to the values from the previous iteration. Because of the non-linearities induced by the positive-negative split (30), the resulting system of non-linear equations is solved with the Newton-Raphson method with a relative termination tolerance $10^{-12}$. The second step involves the update of the damage fields using (2.3) while freezing the kinematic quantities to the ones from the first step. For this step, we employ the FEniCS-SNES solver based on a semi-smooth Newton method for variational inequalities. The iterations terminate when the staggered error in the $i$-iteration

$$\xi^{(i)}=\max\Bigl{\{}\frac{\lVert\boldsymbol{w}^{(i)}-\boldsymbol{w}^{(i-1)}\rVert_{2}}{\lVert\boldsymbol{w}^{(i)}\rVert_{2}},\frac{\lVert\boldsymbol{d}^{(i)}-\boldsymbol{d}^{(i-1)}\rVert_{2}}{\lVert\boldsymbol{d}^{(i)}\rVert_{2}}\Bigr{\}}$$ (23)

drops below the tolerance set as $\mathrm{tol}=10^{-6}$. Here, column matrices $\boldsymbol{w}^{(i)}$ and $\boldsymbol{d}^{(i)}$ collect the vertical deflections and damage variables in all finite element nodes in the $i$-th iteration.

Figure 3: Flowchart of the staggered static phase-field algorithm.

An interested reader is referred to Figure 3 for a flowchart of the algorithm and to, e.g., Ambati et al. (2015); Wu et al. (2020); Dujc and Brank (2022) for more details on the staggered approach.

3 Test and solver setup

In this section, the structural-scale testing data (in Section 3.1), the determination of material parameters (Section 3.2), the setting of the finite element model (Section 3.3), and its verification using a single-layer benchmark (Section 3.4) are collected to provide the basis for the validation studies presented later. Whenever relevant, the observed or simulated failure sequences are compared with analytical approaches based on event tree analyses Bonati et al. (2019).

3.1 Experimental campaign

The experiments were carried out by our colleagues from the Experimental Centre and Department of Steel and Timber Structures at the Faculty of Civil Engineering, Czech Technical University in Prague. Similarly to our previous study Schmidt et al. (2020), the four-point bending setup was considered with the arrangement indicated in Figure 4 for one 5-layer and two 7-layer laminate configurations specified in Table 1.

Figure 4: Scheme of the four-point bending test on multi-layer laminated glass beams. $\overline{w}(t)$ denotes the time evolution of the prescribed displacement and $R(t)$ the corresponding reaction.

	5LG	7LG-1	7LG-2
number of layers	5	7	7
layer thicknesses [mm]	5/2.28/6/0.76/5	5/1.52/8/0.76/8/0.76/5	6/1.52/6/0.76/6/1.52/6
overall thickness [mm]	20.04	29.04	27.8

Table 1: Composition of laminated glass samples and nominal dimensions.

Laminated glass samples were placed into loading device and central deflections were measured by two displacement sensors. Following the standardized procedure specified in the EN1288-3 European Committee for Standardization (2000), the glass specimens were protected by rubber pads at contact with the loading and supporting steel cylinders to prevent stress concentration in these areas. All tests were displacement-controlled with a constant loading rate $\dot{\overline{w}}$ of 1 mm/min until the specimen failure. The prescribed displacements $\overline{w}(t)$, for which we determine the time evolution of the reaction force $R(t)$, increased monotonously in time without any loading/unloading cycles after the progressive failure of individual layers. Additional details on the tests are available in Hána et al. (2018, 2020); Konrád et al. (2022).

Three specimens of each beam sample specified in Table 1 were tested. Table 2 additionally provides the ambient temperatures during the test and the failure sequence leading to the structural collapse. The data acquired for the third specimen of sample 7LG-1 were defective and hence are excluded from further discussion.

sample	specimen	temperature [${}^{\circ}$C]	failure sequence
5LG	1	23.0	5 $\rightarrow$ 3 $\rightarrow$ 1
	2	23.4	5 $\rightarrow$ 3 $\rightarrow$ 1
	3	23.1	5 $\rightarrow$ 3 $\rightarrow$ 1
7LG-1	1	23.3	5 + 7 $\rightarrow$ 3 $\rightarrow$ 1
	2	23.3	5 + 7 $\rightarrow$ 3 $\rightarrow$ 1
	3	defective measurement
7LG-2	1	28.4	7 $\rightarrow$ 5 $\rightarrow$ 1 + 3
	2	29.4	7 $\rightarrow$ 1 + 3 + 5
	3	27.7	7 $\rightarrow$ 5 $\rightarrow$ 3 $\rightarrow$ 1

Table 2: Overview of glass samples and specimens, ambient temperature during the tests, and the failure sequence. Glass layers are indexed according to Figure 1.

Figure 5: Force-displacement diagrams for all laminated glass beam samples and specimens. Solid lines are used to distinguish different samples.

The load-displacement diagrams for all samples appear in Figure 5. Note that because of the indentation of steel cylinders into the rubber pads, the measured data included an initial non-linear branch that is omitted in the results presented in the next section. The gradual fracture of glass layers manifests in the load-displacement curves by the jumps in the reaction force at a constant displacement of the loading head; see also Figure 6 for an illustration of the crack development along the thickness. This information is complemented by the fracture patterns at the sample failure, collected in Figure 7, that document the development of distributed cracks in individual layers.

(a)

(b)

Figure 6: Progressive failure of glass layers in (a) specimen 1 of sample 5LG (failure sequence: $5\rightarrow 3\rightarrow 1$) and (b) specimen 1 of sample 7LG-1 (failure sequence: $5+7\rightarrow 3\rightarrow 1$). The layers are indexed according to Figure 1.

Figure 7: Top view of final fracture patterns on tested laminated glass beams. Row labels correspond to samples, and column labels to specimens.

Notice that in all cases, the failure process initiated in the bottom layer(s) and progressed upward. From the point view of the fault tree considerations Bonati et al. (2019), this suggests that the degree of shear coupling induced by interlayers was sufficient to ensure that the structural response corresponded more to the monolithic limit (in which the bottom layers carry the largest tensile stresses and, therefore, exhibit higher failure probabilities) than to the layered limit with independent layer responses and thus same layer failure probabilities.

3.2 Material parameters

In all the considered laminated glass structures, PolyVinyl Butyral (PVB) is used as an interlayer material and equivalent elastic response is determined from the temperature-depended generalized Maxwell model with parameters obtained in Hána et al. (2019) by an inverse analysis procedure and provided in the repository Schmidt (2023). Uniform but different temperatures were assumed for each simulation, as specified in Table 1 on page 1.

All glass layers are made of annealed float glass, with the standard Young modulus $E=70$ GPa and Poisson ratio $\nu=0.22$, Haldimann et al. (2008). To characterize the variability of the glass tensile strength $f_{\text{t}}$, we combined the data collected on monolithic glass samples in Veer et al. (2009) (24 samples under displacement control tests at a loading rate of 1.0 mm/min) with our previous independent experiments Zemanová et al. (2018) on three-layer laminated samples (10 samples); see the repository Schmidt (2023) for the full data set. Each sample in the set had the same geometry, loading type and rate, and type of edge treatment. Four-point bending tests were conducted using a MTS testing machine. Note that for this type of test, the strength exhibited by the glass sample is mainly governed by the damage caused by the edge finishing process instead of flaws randomly distributed on the surface Castori and Speranzini (2019). For a detailed discussion of the effect of edge treatment on edge strength, see, e.g., the following studies Müller-Braun et al. (2020); Bukieda et al. (2020).

These data were fitted with the Weibull distribution Weibull (1951) with the shape parameter $k=4.64$ and the scale parameter $\lambda=48.47$:

$$f_{\text{t}}\sim\text{Weibull}(k=4.64,\lambda=48.47)~{}\text{MPa};$$ (24)

see Figure 8 for the quality of the fit.

Figure 8: Histogram of tensile strengths complemented with the Weibull fit (left) and corresponding probability plot (right).

The corresponding $5\%$ or $95\%$ quantiles follow from the strength cumulative distribution function $F_{f_{\text{t}}}$ as:

$$f_{\mathrm{t,lo}}=F_{f_{\text{t}}}^{-1}(0.05)\doteq 25.6~{}\text{MPa}, f_{\mathrm{t,hi}}=F_{f_{\text{t}}}^{-1}(0.95)\doteq 61.4~{}\text{MPa},$$ (25)

and these values are considered to be the extreme values used in the deterministic simulations later in Section 4. Note that all the input data were determined from glass samples of a similar length as in our setup to eliminate the stochastic size effect, e.g., (Bažant et al., 2021, Chapter 6) for a general exposition and Bonati et al. (2019, 2020) for dedicated discussion for laminated glass structures. The tensile strength values, therefore, should be understood as effective layer-wise quantities.

3.3 Phase-field model and spatial discretization

As shown in Wu and Nguyen (2018), the regularized form of the dissipated energy (16) introduces an implicit coupling between the fracture energy $G^{\text{c}}$, the tensile strength $f_{\mathrm{f}}$, and the characteristic length $\ell$ in the form:

$$G^{\text{c}}=\frac{8}{3}\frac{f_{\text{t}}^{2}\ell}{E}.$$ (26)

In our simulation, the value of the layer-wise fracture energy is thus determined according to this relation, assuming the distribution of the tensile strength (24) and the characteristic length $\ell=2h_{\mathrm{e}}$, where $h_{\mathrm{e}}$ stands for the uniform element length.

The value of $h_{\mathrm{e}}=0.5$ mm yields a reasonable compromise between the adequate resolution for sharp cracks and the efficiency needed when generating many stochastic realizations. Because of the latter reason, only a symmetric half of the beam was considered in simulations and employed $J=40$ integration points for the evolution of energies $\psi^{\text{e}}_{+,m}$ in Eq. (2.2) in the thickness direction. For the discretization along the beam length, we used the mesh of $1,100$ elements per layer and employed linear basis functions for each unknown field $u_{m}$, $v_{m}$, $\varphi_{m}$, and $d_{m}$ (with one-point reduced integration to eliminate the shear locking). Overall, this setting resulted in typical times in the order of minutes for simulating a single experiment described in Section 3.1 on a conventional laptop with the Intel(R) Core(TM) i7-8550U processor, with larger simulation times needed for seven-layer samples exhibiting distributed cracking.

3.4 Verification benchmark

We consider a monolithic single-layer problem to verify the implementation of the solver, in particular, the cross-section non-penetration condition after fracture. The pre-fracture behavior of the single-layer beam setup from Figure 4 is governed by the elementary theory of thin beams, e.g., (Gross et al., 2018, Chapter 4), which predicts the following relation between the prescribed displacement $\overline{w}$ and the maximum normal stress in the beam $\sigma_{\max}$:

$$\overline{w}=\frac{l^{2}}{3Eh}\left(3\left(\frac{a}{l}\right)-4\left(\frac{a}{l}\right)^{2}\right)\sigma_{\max},$$ (27)

where $E=70$ GPa is the Young modulus of glass Haldimann et al. (2008), $h=20$ mm is the beam height, $l=1,000$ mm denotes the beam span, and $a=400$ mm stands for the location of the loading cylinders relative to supports; recall Figure 4. Setting the maximum stress to the mean value of the tensile strength $\sigma_{\max}=f_{\mathrm{t},\mathrm{mean}}=45$ MPa according to Eq. (24) yields the prescribed displacement at failure $\overline{w}_{\mathrm{f}}=6$ mm. After the fracture, the structure behaves as the mechanism consisting of two rigid bodies connected by an eccentric hinge, see Figure 9a, leading to the jump of the horizontal displacement at the mid-span

$$\llbracket u\rrbracket=h\sin\alpha, \alpha=\arctan\frac{\overline{w}_{\mathrm{f}}}{a};$$ (28)

for the current data, we obtain $\alpha\doteq 1.500\cdot 10^{-2}$ rad and $\llbracket u\rrbracket\doteq 0.300$ mm.

(a)	(b)

Figure 9: Post-fracture kinematics of single-layer glass beam (a) without and (b) with cross-section overlapping.

The same problem was simulated using the finite element model, with the Young modulus of the central finite element decreased by $0.1\%$ to ensure that a single crack localized at the mid-span. We obtained the values of $\overline{w}_{\mathrm{f}}^{\mathrm{FE}}\doteq 6.006$ mm and $\llbracket u\rrbracket^{\mathrm{FE}}\doteq 0.3003$ mm. Among others, these results confirm the cross-section overlapping has been eliminated by the positive-negative split of the stored energy (15b), without which $\llbracket u\rrbracket^{\mathrm{FE}}=0$ mm; see Figure 9b.

4 Results

The present section collects the results of numerical simulations of the progressive failure of laminated glass samples introduced in Section 3.1. We depart from validating our implementation of a single-layer beam benchmark with a closed-form solution in Section 3.4. Section 4.1 briefly discusses the results with mean strength data, while Section 4.2 assesses the range of low- and high-strength layers and the corresponding failure behavior. Finally, in Section 4.3, we collect the results of Monte-Carlo simulations and compare them with the available experimental data.

4.1 Mean value analysis

We begin our discussion by briefly comparing experimentally determined load-displacement curves with simulation data for deterministic glass layer strengths with mean values $f_{\mathrm{t,mean}}\doteq 45$ MPa (recall Section 3.2). The results shown in Figure 10 confirms that the initial linear response is well-captured by the numerical model, while the maximum displacement $\overline{w}$ and reaction force $R$ are overestimated for samples 5LG and 7LG-2 and underestimated for sample 7LG-1. More importantly, all models predict perfect elastic-brittle failure through simultaneous fracture of all glass layers by a single crack. Such behavior is consistent with our previous study Schmidt et al. (2020).

Figure 10: Force-displacement diagrams for all laminated glass beam samples and specimens. Solid lines show experimental results, and dashed lines show simulated results with average tensile strength in glass layers.

4.2 Combinatorial analysis

To investigate whether the difference in failure modes among individual specimens and distributed cracking observed in Section 3.1 can be explained by variable strength distribution, this section collects the mechanical response of all $2^{3}$ and $2^{4}$ combinations of the strengths $f_{\mathrm{t,lo}}$ and $f_{\mathrm{t,hi}}$ from Eq. (25) for the samples 5LG in Section 4.2.1, and 7LG-1 in Section 4.2.2, respectively (the results for the sample 7LG-2 are similar thus omitted for the sake of brevity).

4.2.1 5-layer beam

All simulation results for the sample 5LG are collected in Figure 11 and are visualized using the force-displacement diagrams (Figure 10(a)) and bar charts showing the integrity of individual glass layers until failure (Figure 10(b)). The particular combinations are distinguished by the distribution of the low and high strengths in the top-to-bottom direction, according to Figure 1. For example, ’lo-hi-hi’ combination indicates a structure with the top glass layer of strength $f_{\mathrm{t,lo}}$ and the middle and bottom layers of strength $f_{\mathrm{t,hi}}$.

(a)

(b)

Figure 11: Evolution of (a) reaction $R$ and (b) layer integrity as a function of the prescribed displacement $\overline{w}$ for all combinations of minimum and maximum glass layer strengths in the 5-layer beam sample 5LG. The layers are ordered in the top-to-bottom direction and enumerated by an index $m$ according to Figure 1.

lo-lo-lo
	$\rightarrow$		$\rightarrow$
lo-hi-lo
	$\rightarrow$		$\rightarrow$
hi-hi-lo
	$\rightarrow$		$\rightarrow$

Figure 12: Crack development for three representative combinations of minimum and maximum glass layer strengths: lo-lo-lo, lo-hi-lo, and lo-hi-hi in 5-layer beam sample 5LG. The violet color corresponds to an undamaged cross-section ($d=0$) and yellow to a fully fractured cross-section ($d=1$). The layers are ordered in the top-to-bottom direction and enumerated by an index $m$ according to Figure 1.

As expected, the simulated responses remain confined between the uniform combinations of low (red curve in Figure 10(a)) and high (blue curve) strengths for which the crack propagates through all three layers simultaneously; the same behavior is observed for the lo-lo-hi combination. These results correspond to our previous findings in Schmidt et al. (2020). All the remaining combinations exhibit progressive collapse and gradual crack development in individual layers, although the difference from the sudden failure is only minor for the hi-lo-lo combination. We particularly highlight the lo-hi-hi, hi-hi-lo, and lo-hi-lo combinations that may represent the situation in which one or both outer glass layers are unintentionally scratched during transportation and manipulation. Nevertheless, their failure behavior differs. For the lo-hi-lo combination, scratched from both sides, the cracks first propagate through the bottom layer ($m=5$). The top layer ($m=1$) fractures under an additional load increase, followed by structural failure after the rupture of the middle layer ($m=3$). A different behavior is observed for configurations with a single low-strength layer adjacent to the top or bottom surface. Here, the first crack appears in the low-strength layer. Still, upon an increasing load, an array of equally spaced cracks starts to develop until the remaining two layers fail by the simultaneous formation of a single crack per layer. The lo-hi-lo combination exhibits similar behavior, except the fragmentation remains confined to the middle layer. The fracture evolution for the selected configurations is additionally shown in Figure 12.

4.2.2 7-layer beam

An analogous analysis for the 7-layer beam 7LG-1 reveals the same mechanisms of the failure processes, as shown for the selected layer combinations in Figure 12(a) and all layer combinations in Figure 12(b).

(a)

(b)

Figure 13: Evolution of (a) reaction $R$ and (b) layer integrity as a function of the prescribed displacement $\overline{w}$ for the (a) selected and (b) all combinations of minimum and maximum glass layer strengths in the 7-layer beam sample 7LG-1. The layers are ordered in the top-to-bottom direction and enumerated by an index $m$ according to Figure 1.

In particular, the uniform strength distribution through the beam thickness leads to a brittle failure by a single crack aligned in all glass layers. The same brittle response occurs for configurations lo-hi-hi-hi, hi-lo-lo-hi, lo-lo-lo-hi, and hi-lo-lo-lo with a uniform strength distribution in the adjacent glass layers; cf. Figure 12(b). The alternating distributions of minimum and maximum layer strengths lead to a more progressive collapse and even the formation of equally-spaced cracks in specific glass layers, as illustrated in Figure 14.

lo-hi-lo-hi
	$\rightarrow$		$\rightarrow$

Figure 14: Crack development for the lo-hi-lo-hi combination of glass layer strengths in the 5-layer laminated glass beam (sample 7LG-1). The violet color corresponds to an undamaged cross-section ($d=0$) and yellow to a fully fractured cross-section ($d=1$). The layers are ordered in the top-to-bottom direction and enumerated by an index $m$ according to Figure 1.

4.3 Stochastic analysis

The results of the combinatorial analysis confirmed the numerical model’s potential to predict the progressive cracking phenomena; this section, therefore, deals with the extension to randomized simulations. Note that the results for all samples are presented simultaneously because, according to the combinatorial analysis, 5- and 7-layer samples revealed similar failure mechanisms. Recall that the outcome of the experimental campaign in Section 3.1 involved the measured force-displacement diagrams (shown in Figure 5) with failure sequences (Table 2 and Figures 6 and 7). These data are compared here with the results of 1,000 Monte-Carlo simulations, in which the layer strengths were sampled from the Weibull distribution displayed in Figure 8.

Figure 15: Results and statistics reaction-displacement curves obtained from $1,000$ Monte-Carlo (MC) simulations of laminated glass beam samples 5LG, 7LG-1, and 7LG-2. The thick red curves correspond to experimental data. Thin gray curves denote individual realizations, blue curves the mean values, and the green curves the upper and lower $5~{}\%$ quantiles determined from the MC simulations.

The simulated force-displacement diagrams are shown in Figure 15 with thin gray lines, accompanied by the experimental results with thick red lines. The range of predicted responses is substantial and aligns well with the results from the combinatorial analysis conducted in Section 4.2 as seen by comparing Figure 15 with Figures 10(a) and 12(a). To analyze the simulated data, the reaction values $R$ were sampled for a set of prescribed displacements $\overline{w}$ at 1/30 mm intervals. For each $\overline{w}$ value, we constructed an empirical distribution function from the complete set of 1,000 Monte-Carlo simulations, from which the median values (shown in Figure 15 as green curves) and the $5\%$ and $95\%$ quantiles (blue solid and dashed curves) were extracted. Note that the $5\%$ quantile response exhibits gradual cracking, and the $95\%$ quantile transitions to the brittle response. For samples 5LG and 7LG-2, the experimental response falls between the two quantiles, and the median response approximates the experimental data well. For sample 7LG-1, the median response underestimates the experimental response, which is close to the $95\%$ quantile curve. However, the smoothing nature of the median response (especially near the peak) and the consideration of extreme value statistics over all realizations disregard the information about the failure sequences, and more refined interpretations of the Monte-Carlo data are needed.

Figure 16: Top: Histograms of displacements inducing the first crack $\overline{w}_{\mathrm{i}}$ (blue) and the final crack $\overline{w}_{\mathrm{f}}$ (red) displacements and their Weibull fits for the 5-layer sample 5LG and 7-layer samples 7LG-1 and 7LG-2. Bottom: Plots of the final $\overline{w}_{\mathrm{f}}$ against the first $\overline{w}_{\mathrm{i;}}$ crack displacements.

For this purpose, let us analyze displacements $\overline{w}_{\mathrm{i}}$ leading to the formation of the initial cracked layer (corresponding to the occurrence of a fully-cracked cross-section with $d=1$ in at least a single layer) and $\overline{w}_{\mathrm{f}}$ leading to the complete structural failure (the occurrence of a fully-cracked cross-section at all layers). The resulting data for all three samples are visualized in Figure 16 using two means. The first one is based on histograms of the initial crack displacements $\overline{w}_{\mathrm{i}}$ (in blue) and the complete failure $\overline{w}_{\mathrm{f}}$ (in red), along with their fits with the Weibull distributions (24). The fitted distributions display visible differences between the two quantities as the ratios between the most likely values (modes) $\overline{w}^{\mathrm{mod}}_{\mathrm{i}}/\overline{w}^{\mathrm{mod}}_{\mathrm{f}}$ reach 0.93 (for sample 5LG), 0.93 (7LG-1), and 0.91 (7LG-2). However, a different conclusion follows from the bottom graphs in Figure 16, in which we plotted the failure displacement against the displacements at the initial failure for all realizations. These graphs show that for most realizations $\overline{w}_{\mathrm{i}}=\overline{w}_{\mathrm{f}}$, indicating the brittle failure by crack(s) running through all glass layers.

Figure 17: Probability of occurrence of different failure sequences for the 5-layer sample 5LG and 7-layer samples 7LG-1 and 7LG-2. The layers are ordered in the top-to-bottom direction.

This claim is further substantiated by Figure 17, showing the probabilities of occurrence for the three most frequent failure sequences. The most common failure scenario for all samples involves simultaneous failure for all layers, with a higher likelihood for the 5-layer sample than for the 7-layer configurations. The second most likely failure sequence proceeds from the collapse of the bottom layer, followed by the simultaneous failure of all remaining layers. This scenario is the closest to the ones observed in experiments for all specimens, in which the failure initiated in the bottom layer gradually propagated, layer-by-layer, to the top surface under small displacement increments (exceptionally through the simultaneous failure of two layers); recall Section 3.1 for additional details. Hence, the response predicted by the proposed model is less ductile due to neglecting variability in material properties along the lengths of individual layers; see Biolzi et al. (2017); Casolo and Diana (2018) for related results in this direction.

Similarly to the experimental results, the numerical simulations confirm that the mechanical response of the laminated beam is much closer to the monolithic than to the laminated limit failure-wise Bonati et al. (2019). In particular, the scenarios corresponding to the failure initiation in the two bottommost layers account for more than $80\%$ of the simulated cases ($95.8\%$ for 5LG, $89.3\%$ for 7LG-1, and $81.3\%$ for 7LG-2), cf. Figure 17.

Figure 18: Layer-wise distribution of tensile strengths corresponding to the second most likely failure sequences $5\rightarrow 1+3$ (sample 5LG) and $7\rightarrow 1+3+5$ (samples 7LG-1 and 7LG-2). The layers are ordered in the top-to-bottom direction and enumerated by an index $m$ according to Figure 1 and the histograms are fitted with the Weibull distributions.

Nevertheless, we find it instructive to assess which strength combinations lead to the two-step cracking mechanisms by collecting the layer strength values for the realizations exhibiting the second most likely failure sequences in Figure 18. The data show that for all samples, the bottom layers exhibit the lowest strengths, followed by the top and interior layers. Such results support the fact that, in real-world applications, the outer layers have inferior mechanical properties to the interior ones because they are exposed to the surrounding environment and hence more prone to damage from transport and manipulation; see (Bonati et al., 2020, Section 2) and (Alter et al., 2017, Section 6.2) for further discussion.

5 Conclusions

This study used randomized phase field models to address simulations of distributed and gradual fracture patterns in laminated glass beams. As a building block, we used a layer-wise dimensionally-reduced phase-field model identified in our previous study Schmidt et al. (2020) as an efficient and accurate alternative to more refined two-dimensional plane strain or plate formulations. In particular, we adopted Timoshenko kinematics for each layer, ensured interlayer compatibility by means of an elimination technique, approximated the time- and the temperature-dependent response of an interlayer with a quasi-elastic model, and performed the numerical integration in the thickness direction to ensure that glass layer cross-sections were non-penetrable. Variability in the material properties was accounted for through a simple stochastic model, using the assumption that the strengths of individual layers are independent identical random variables with the Weibull distribution. Parameters of this distribution were determined from the data available in the literature and our own data for same-sized three-layer specimens.

From the performed combinatorial and Monte-Carlo simulations and experimental campaign for 5- and 7-layer glass beams, we concluded that:

1. With a suitable combination of low- and high-strength layers, the model exhibits distributed cracking via an array of localized cracks.

2. The second most likely failure patterns, which are in reasonable agreement with experimental results, correspond to lower strengths in the bottom and top layers, which are most exposed to the environment.

3. The single-layer benchmark proposed here can be used to test implementations of phase-field models for beam structures.

4. The model fails to predict gradual cracking in individual layers in the bottom-up direction, as observed in our experiments. We attribute this limitation to neglecting strength variability along the length of the layers, which promotes the crack alignment in the adjacent layers.

Extensions of this work could involve an extensive experimental campaign to obtain a representative set of structural responses, but this is beyond the scope of this paper. On the modeling side, a more refined description of strength variability would involve modeling strengths using random fields Gerasimov et al. (2020); Hai and Li (2022); Nagaraja et al. (2023) and ensuring objectivity for the probability of failure under mesh refinement in the spirit of recent studies on crack band Gorgogianni et al. (2022) and phase-field Hai et al. (2023); Wu et al. (2023) models. We plan to address this limitation in a future study, benefiting from the existing works on discrete element/peridynamics approaches to laminated glass fracture Biolzi et al. (2017); Casolo and Diana (2018) and insights from dedicated physically-based models for glass strength Symoens et al. (2023); Rudshaug et al. (2023) and stochastic size effect considerations Bonati et al. (2019, 2020).

Author contributions

J. Schmidt: Conceptualization, Methodology, Software, Formal analysis, Data Curation, Writing - Original Draft, Visualization; A. Zemanová: Methodology, Formal analysis, Data Curation, Writing - Original Draft, Writing - Review & Editing; J. Zeman: Conceptualization, Methodology, Formal analysis, Writing - Original Draft, Writing - Review & Editing

Appendix A Derivatives of elastic energies

Because the polymer interlayers are assumed to be quasi-elastic, the energy derivative (19b), for $m=2,4,\ldots,M$, attains the standard linear form:

$$\boldsymbol{\nabla}_{\boldsymbol{u}_{m}}\Psi^{\text{e}}_{m}(t_{n},\boldsymbol{u}_{m})\cdot\delta\boldsymbol{u}_{m}= \,E_{m}\int_{0}^{L}\Bigl{(}A_{m}u_{m}^{\prime}\delta u_{m}^{\prime}+I_{m}\varphi^{\prime}_{m}\delta\varphi^{\prime}_{m}\Bigr{)}\,{\mathrm{d}}x$$ (29)

$$+G_{m}(t_{n})A^{*}_{m}\int_{0}^{L}\left(\varphi_{m}+w^{\prime}_{m}\right)\left(\delta\varphi_{m}+\delta w^{\prime}_{m}\right)\,{\mathrm{d}}x,$$

For the glass layers, the explicit form of (19a) for $m=1,3,\ldots,M$ is more involved:

$$\boldsymbol{\nabla}_{\boldsymbol{u}_{m}}\Psi^{\text{e}}_{m}(\boldsymbol{u}_{m},d_{m})= E_{m}\,b\int_{0}^{L}\sum_{j=1}^{J}\left((1-d_{m})^{2}\langle u^{\prime}_{m}+\varphi^{\prime}_{m}z_{m,j}\rangle_{+}-\langle u^{\prime}_{m}+\varphi^{\prime}_{m}z_{m,j}\rangle_{-}\right)\times$$ (30)

$$\left(\delta u_{m}^{\prime}+\delta\varphi^{\prime}_{m}z_{m,j}\right)\Delta z_{m}\,{\mathrm{d}}x+G_{m}A_{m}^{*}\int_{0}^{L}(1-d_{m})^{2}\left(\varphi_{m}+w_{m}^{\prime}\right)\left(\delta\varphi_{m}+\delta w_{m}^{\prime}\right)\,{\mathrm{d}}x,$$

because of the numerical integration and the non-smooth non-linearity induced by the positive and negative parts $\langle\bullet\rangle_{+}$ and $\langle\bullet\rangle_{-}$.