# Prediction of Pre- and Post-Breakage Behavior of Laminated Glass Using a Phase-Field Damage Model
## 3 Phase-Field Model
### Mindlin Plate Model  
A full 3D model is computationally demanding which suggests that application of a more effective spatially reduced plate model. We start from a single layer model which will be subsequently extended to laminated plates. The Mindlin-Reissner theory, which assumes that the cross-section remains planar after deformation but not necessarily perpendicular to the midsurface of the plate, is adopted. Commonly, the model is obtained by assuming a linear distribution of displacements and rotations and subsequent through-thickness integration. Owing to damage, this integration must be carried out numerically. The elastic energy then receives the form  
\[\Psi^{+}=\frac{1}{2}\int_{\Omega_{\mathrm{s}}}g(d)\left(\sum_{i} \Delta z_{i}\mathbf{\varepsilon}_{i}:\mathbf{\sigma}_{i}^{+}+\xi h\mathbf{\gamma}\cdot\mathbf{ \tau}\right)\mathrm{d}\mathbf{x}_{\mathrm{s}}, \tag{14}\] \[\Psi^{-}=\frac{1}{2}\int_{\Omega_{\mathrm{s}}}\sum_{i}\Delta z_{ i}\mathbf{\varepsilon}_{i}:\mathbf{\sigma}_{i}^{-}+(1-\xi)h\mathbf{\gamma}\cdot\mathbf{\tau} \mathrm{d}\mathbf{x}_{\mathrm{s}}, \tag{15}\]  
where the overall thickness \(h\) is numerically divided into individual layers \(i\) of thickness \(\Delta z_{i}\) with in-plane stress \(\mathbf{\sigma}_{i}\) and strain \(\mathbf{\varepsilon}_{i}\) vectors, see Figure 5. The midsurface \(\Omega_{\mathrm{s}}\) with the \(\mathbf{x}_{\mathrm{s}}=\{x,y\}\) coordinate system is the new spatially reduced domain of interest. Point out that the parameter \(\xi\) is an auxiliary material constant which indicates that the amount of shear subjected to damage, i.e., whether the shear stress is still transmitted across the crack. Notice that only the part of out-of-plane shear contribution that enters the active part of strain energy density \(\psi^{+}\) is degraded by function \(g(d)\).  
Figure 4: Flowchart of the central difference phase-field algorithm.  
The constitutive equations for in-plane stress \(\sigma\) and strain \(\epsilon\) fields and out-of-plane shear stresses \(\tau\) and strains \(\gamma\) are provided by  
\[\sigma(\epsilon) =\mathbf{D}^{\text{PS}}:\epsilon, \tag{16}\] \[\tau(\gamma) =2\mu\mathbf{I}:\gamma, \tag{17}\]  
where \(\mathbf{D}^{\text{PS}}\) represents the plane stress stiffness tensor expressed in terms of [PERSON]'s coefficients \(\lambda\) and \(\mu\) and \(\mathbf{I}\) is the identity tensor.  
The strains \(\epsilon\) and \(\gamma\) follow from standard strain-displacement relations and are written in terms of in-plane displacements \(u_{\text{s}}\), out of plate deflection \(w\), and rotations of the cross section \(\varphi\) as  
\[\epsilon(\mathbf{x}_{\text{s}},z) =\
abla_{\text{s}}\mathbf{u}(\mathbf{x}_{\text{s}})+\
abla_{\text{s}} \mathbf{S}\,\mathbf{\varphi}(\mathbf{x}_{\text{s}})z, \tag{18}\] \[\gamma(\mathbf{x}_{\text{s}}) =\frac{1}{2}(\
abla w(\mathbf{x}_{\text{s}})+\mathbf{S}\,\mathbf{ \varphi}(\mathbf{x}_{\text{s}})), \tag{19}\]  
where the matrix \(\mathbf{S}\) is given by  
\[\mathbf{S}=\begin{bmatrix}0&1\\ -1&0\end{bmatrix}. \tag{20}\]  
Equation (12) is also subjected to spatial reduction via through-thickness integration to obtain  
\[\int_{\Omega_{\text{s}}}\left(\frac{hG_{\text{f}}}{c_{\text{a}}}\left(\frac{ \alpha^{\prime}(d)}{l}+2l\
abla d\
abla\delta d\right)+g^{\prime}(d)\int_{-h/ 2}^{h/2}\psi^{+}\,\mathrm{d}z\right)\mathrm{d}\mathbf{x}=0,\quad\forall\,\delta d, \tag{21}\]  
where the integrated active part of the strain energy density  
\[\int_{-h/2}^{h/2}\psi^{+}\mathrm{d}z, \tag{22}\]  
drives the damage evolution. In this regard, it appears useful to replace the strain energy density \(\psi^{+}\) by a constant \(\psi^{+}_{\text{sur}}\) representing the strain energy density calculated on the basis of the maximum tensile stress. This rewrites Equation (21) as  
\[\int_{\Omega_{\text{s}}}\left(\frac{hG_{\text{f}}}{c_{\text{a}}}\left(\frac{ \alpha^{\prime}(d)}{l}+2l\
abla d\
abla\delta d\right)+g^{\prime}(d)\psi^{+}_ {\text{sur}}\cdot h\right)\mathrm{d}\mathbf{x}=0,\quad\forall\,\delta d. \tag{23}\]