## 3 Verification of the Combined TCK-JH2 Model
### Finite Element Model  
The drop-weight test described above was replicated numerically to verify the glass material model. Figure 3a shows the parts of the finite element model. The finite element model consists of the impactor, two layers of glass, the PVB interlayer, and the rubber cushion. The dimensions of these parts are in accordance with tests. However, due to symmetry, only one quarter of the model was simulated to decrease the analysis time. In numerical simulation, the combined TCK-JH2 constitutive model was employed to describe the material model of glass. The constants of this model are shown in Table 2. Based on the combined TCK-JH2 model, a user-defined material subroutine (VUMAT) was compiled by FORTRAN and connected with commercial finite element code ABAQUS/Explicit to calculate the tensile and compressive damage. The volumetric strain \(\varepsilon_{\
u}\) was used as a criterion to determine to element status. When \(\varepsilon_{\
u}>0\), it is determined that the material is under tension. Otherwise, it is determined to be under compression. Two damage scalars were adopted: \(D_{\mathrm{t}}\) (under tension) and \(D_{\mathrm{c}}\) (under compression). From tension to compression, the stiffness of the element should be multiplied by \((1-D_{\mathrm{c}})\). However, from Compression to tension, the stiffness of the element should be multiplied by \((1-D_{\mathrm{t}})(1-D_{\mathrm{c}})\). Details of implementation process of the combined TCK-JH2 model is presented by a flow chart in the Appendix A. The interlayer material PVB was revealed to be a rate-dependent material in laboratory tests. It is common that the elastic model is employed to describe PVB material under static loading. However, when PVB is under high strain rate, hardening and strengthening occur. For this reason, the elastoplastic model [7; 13; 14] is used in this paper. Yield stress (\(\sigma_{\mathrm{yield}}\)), initial Young's modulus (\(E_{\mathrm{initial}}\)), and failure stress (\(\sigma_{\mathrm{failure}}\)) of PVB at different strain rates are described as:  
\[\sigma_{\mathrm{yield}} =2.167(\dot{\epsilon})^{0.399}\mathrm{MPa}, \tag{18}\] \[E_{\mathrm{initial}} =30.591(\dot{\epsilon})^{0.271}\mathrm{MPa},\] (19) \[\sigma_{\mathrm{failure}} =27.689(\dot{\epsilon})^{0.040}\mathrm{MPa}. \tag{20}\]  
The density, Poisson's ratio, and plastic modulus (\(E_{\mathrm{yield}}\)) of PVB used in this study is shown in the Table 4[7; 13; 14]. The elastic model is used to describe the rubber cushion in this study. The Young's modulus and Poisson's ratio are 2.3 MPa and 0.495 respectively. The impactor was regarded as a non-deformed rigid object. The Young's modulus and Poisson's ratio are 210 GPa and 0.3, respectively. The average mesh size was 1 mm. At the area near the impact regions, glass layer, PVB interlayer and impactor were meshed with denser elements. As a result, the mesh size around the contact point was 0.5 mm. The hourglass mode was carefully checked in the mesh model. In order to prevent the zero-energy mode of the element, the mesh model with different number of elements through-thickness were pre-simulated. When the proportion of artificial strain energy is 10% less than the internal energy, the accuracy of the mesh model is satisfied. Through some pre-simulations, the element size in the thickness direction was set to 0.3 mm. In the thickness direction of PVB interlayer, it was meshed into four elements. An eight-node linear brick element with reduced integration was adopted in this finite \begin{table}
\begin{tabular}{c c c} \hline \hline \(\rho\) (kg/m\({}^{3}\)) & \(\
u\) & \(E_{\mathrm{yield}}\) (MPa) \\ \hline
**1100** & 0.495 & 11 \\ \hline \hline \end{tabular}
\end{table} 
Table 4: Polyvinyl Butyral (PVB) material properties. element model. The Figure 3b shows the mesh model of the glass layers and the PVB interlayer shares the same meshing pattern. Tie constraint was defined between the contact surfaces of glass and PVB interlayer.