## 3 Verification of the Interface Model with a Benchmark Problem  
The interlaminar stress calculation strategy based on linear elastic traction-separation behavior was validated by the consideration of a benchmark problem. This benchmark problem was analyzed by [PERSON] [17], where the 3D exact elasticity solution of a square, simply supported, symmetric cross-ply [0/90/0] laminated plate under a bisinusoidally distributed transverse load \(q\), was provided. As shown in Figure 2, the span and thickness of the laminate is denoted by \"\(a\)\" and \"\(h\)\", and \(q\) can be described as  
\[q=q_{0}\sin(\pi x/a)\sin(\pi y/a) \tag{5}\]  
where \(q_{0}\) is a constant. A laminate of span-to-thickness ratio \(S=a/h=4\) are considered, and can be considered as a thick laminate. The thickness of the ply is 0.25 mm when \(a=3\) mm. The following lamina material properties are used: \(E_{1}=25\)\(\times\)\(10^{6}\) psi, \(E_{2}=E_{3}=10^{6}\) psi, \(G_{12}=G_{23}=0.2\)\(\times\)\(10^{6}\) psi, \(\
u_{12}=\
u_{13}=\
u_{23}=0.25\). The resulting stresses are normalized according to the following formulae:  
\[\left(\overline{\tau}_{xz},\overline{\tau}_{yz}\right)=\frac{1}{q_{0}S}\left( \tau_{xz},\tau_{yz}\right) \tag{6}\]  
The interlaminar stresses analyses were performed using the commercial FE code, Abaqus (Version 6.11, Dassault Systemes Simulia Corp., Providence, RI, USA, 2011). A 3D 8-node linear brick, with an incompatible modes element, was used to model the ply. The interface element based on the linear traction-separation description was introduced into two adjacent piles. This constitutive response is readily available in Abaqus and therefore does not require user-defined subroutines. This 3D FE model is generated using uniform meshes, where the in-plane mesh consisted of 30 \(\times\) 30 elements. For comparison purposes, two kinds of mesh density in the thickness directions are considered, namely coarse mesh and fine mesh, respectively. For the case, \(a=3\) mm, \(S=4\), the former of which is just one element for each ply in the thickness direction, and the latter of which indicates four elements for each ply. According to the recommended rage of the interface stiffness by Reference [38], \(K\) is taken as \(5\times 10^{6}\) N/mm\({}^{3}\). The distribution of various non-dimensionalized stresses \(\tau_{xz}\), \(\tau_{yz}\) respectively, through the thickness of the plate at the position where each stress attained a maximum is shown in  
Figure 2: Square laminated plate subjected sinusoidally distributed transverse load. Figure 3. In this procedure, the interlaminar stresses were calculated by the interface element, which used the traction-separation description. The transverse stresses in layers were analyzed based on 3D elasticity and extrapolated from the element integration point. As shown in Figure 3, for the coarse mesh, the percentage error of \(\tau_{xz}\) at the point of interface, compared with the exact solution is 4.3%, and \(\tau_{yz}\) is 2.9%. However, the percentage error at the points of maximum value of \(\tau_{xz}\) and \(\tau_{yz}\) compared with the exact solution is 21%, 20.9%, respectively. Refined mesh in thickness direction can improve the precision of the transverse stress calculation precision within the ply, for the maximum value of \(\tau_{xz}\) and \(\tau_{yz}\), the percentage error compared with the exact solution is 3.1% and 5.3%, respectively, and the stresses obtained at interface are very close to the exact solution. In summary, for the transverse stresses at interface (namely interlaminar stresses), both of these meshing strategies showed a reasonable accuracy (within 5%). It was noted that the stress distributions within the ply predicted by the coarse mesh model showed considerable error for thick plate. As the transverse stresses at the interface are only considered in this paper, the modeling strategy with coarse mesh in thickness direction is still used in this paper.