# Prediction of Pre- and Post-Breakage Behavior of Laminated Glass Using a Phase-Field Damage Model
## 3 Phase-Field Model
### Model Selection and Decomposition  
The phase field model is determined by the decomposition of strain energy density into active \(\psi^{+}\) and passive \(\psi^{-}\) parts and by the selection of function \(\alpha(d)\).  
The most physically sound response is simulated using spectral decomposition of functional (5) proposed in [19]. However, this decomposition contains singularities in the finite element application and resulted in instabilities in our simulations. On the hand, stable solutions were experienced when considering the decomposition based on a volumetric-deviatoric split as suggested in [20]. The active and passive contributions to the strain energy density \(\psi\) are then given by  
\[\psi^{+}(\mathbf{u}) = \frac{1}{2}K\langle\mathrm{tr}(\varepsilon)\rangle^{2}+\mu \varepsilon_{\mathrm{D}}:\varepsilon_{\mathrm{D}}, \tag{29}\] \[\psi^{-}(\mathbf{u}) = \frac{1}{2}K\langle-\mathrm{tr}(\varepsilon)\rangle^{2}, \tag{30}\]  
where \(\varepsilon_{\mathrm{D}}\) is a deviatoric part of the strain tensor \(\varepsilon\), and \(K\) is the bulk modulus.  
As for the dissipation function, the model is valid for an arbitrary choice of \(\alpha(d)\). The present formulation follows [PERSON] [21] and considers a simple relation  
\[\alpha(d)=d,\quad c_{\alpha}=8/3, \tag{31}\]  
with the degradation function defined as  
\[g(d)=(1-d)^{2}. \tag{32}\]  
Unlike common representation of a dissipation function \(\alpha(d)=d^{2}\) suggested by [PERSON] et al. [17], which considers damage evolution already at the onset of loading, the adopted linear version allows for postponing the damage initiation beyond this point, and so not affecting the initial transfer of energy from the impactor into the plate.  
For the present choice of \(\alpha(d)\), the fracture energy \(G_{\mathrm{f}}\) can be written in terms of the tensile strength \(f_{\mathrm{t}}\), Young's modulus \(E\), and the length scale parameter \(l\) as [22]  
\[G_{\mathrm{f}}=\frac{8}{3}\frac{lf_{\mathrm{t}}^{2}}{E}. \tag{33}\]  
While strictly valid for one-dimensional (1D) setting, Equation (33) is generally adopted also for simulations in higher dimensions. Owing to the presence of length scale parameter \(l\), typically equal to two times the characteristic length of the smallest element in the finite element mesh [19]; this relation is also mesh dependent.  
For a typical edge length of a triangular element of 2 mm employed for the analysis of laminates discussed in Section 2, the fracture energy stored in Table 3 may considerably exceed a typical value of \(G_{\mathrm{f}}=8\) Jm\({}^{-2}\) used for glass. While the value of \(f_{\mathrm{t}}=45\) MPa corresponds to the characteristic quasi-static strength in tensile bending defined by the European Standard EN 16612 [23], the larger values, e.g., \(f_{\mathrm{t}}^{\mathrm{in}}=145\) MPa, used particularly in the analysis of a single glass ply, represents the initial value of tensile strength we adopted to postpone the onset of damage in numerical simulations and consequently to arrive at the glass response comparable with experimental observations. Point out that even larger values of the initial tensile strength were considered in [10]. Since at the onset of damage this value almost immediately drops down to the actual strength of 45 MPa, the associated fracture energy enters the analysis only at a very beginning stage of fracture process; see Section 4 ahead for further details. Nevertheless, realistic values of \(G_{\mathrm{f}}\) would still call for much smaller elements to avoid snap-back at a material level. On the contrary, this would lead to computationally unfeasible simulations. Thus, only the large values around 300 Jm\({}^{-2}\) to simulate the crack propagation were tested in all present calculations.