## 5 3D FE Model and Schematic of the Problem  
To simulate the dynamic response of the viscoelastic composite plate in 3D, a commercially available numerical analysis package Comsol Multiphysics [43] based on the Finite Element method is used. Preliminary tests on guided waves dispersion curves were performed using the abovementioned software and revealed that all simulations matched well with the theory and also agreed fairly well with experiments [44]. In the proposed simulations, the PDE mode in coefficient form is used to model the composite plate, and Neumann boundary conditions have been applied on all faces of the plate. At the interface between the conical damage zone and the composite plate, the continuity of stress and displacement is used to solve the numerical model. Besides, one external routine is developed for computing normal displacement fields at the midplane around the point impact zone in a circular region at 360 different angular positions. The equation of the dynamic equilibrium for orthotropic viscoelastic material, with nine complex moduli, is written and solved in the Fourier domain in the 3D FE model as: [13]  
\[\sum_{j,\,k,l=1}^{3}C_{kjl}\Bigg{[}\frac{\partial^{2}u_{j}}{\partial x_{k} \partial x_{l}}\Bigg{]}+\rho\omega^{2}u_{i}=0,\,\,\,i=1,\,2,\,3, \tag{1}\]  
where \(u_{i}\) is the Fourier transform of the components of the displacement vector \(u\), with \(i=1,\,2,\,3\), representing the direction of the coordinate axis as shown in Figure 5. Here \(C_{kjl}\) is the component of the complex stiffness tensor, \(\rho\) is the material density and \(\omega\) is the angular frequency. The above partial differential equations (PDE) have been written in the following form as: [13]  
\[\
abla(c\
abla u)-au=0 \tag{2}\]  
where \(c\) is a 3 \(\times\) 3 matrix. For an orthotropic material, the axis of symmetry coincides with the coordinate axis, and this matrix is composed of the following sub-matrices \(c_{kl}\), \(k,l=1,2\), 3 as shown is Equation (3): [13]\[c_{11}=\left( \begin{array}{ccc}C_{11}&0&0\\ 0&C_{66}&0\\ 0&0&C_{55}\end{array} \right),\,c_{12}=\left( \begin{array}{ccc}0&C_{12}&0\\ C_{66}&0&0\\ 0&0&0\end{array} \right),\,c_{13}=\left( \begin{array}{ccc}0&0&C_{13}\\ 0&0&0\\ C_{55}&0&0\end{array} \right),\] \[c_{21}=\left( \begin{array}{ccc}0&C_{66}&0\\ C_{21}&0&0\\ 0&0&0\end{array} \right),\,c_{22}=\left( \begin{array}{ccc}C_{66}&0&0\\ 0&C_{22}&0\\ 0&0&C_{44}\end{array} \right),\,c_{23}=\left( \begin{array}{ccc}0&0&0\\ 0&0&C_{23}\\ 0&C_{44}&0\end{array} \right), \tag{3}\] \[c_{31}=\left( \begin{array}{ccc}0&0&C_{55}\\ 0&0&0\\ C_{31}&0&0\end{array} \right)\!c_{32}=\left( \begin{array}{ccc}0&0&0\\ 0&0&C_{44}\\ 0&C_{32}&0\end{array} \right),\,c_{33}=\left( \begin{array}{ccc}C_{55}&0&0\\ 0&C_{44}&0\\ 0&0&C_{33}\end{array} \right),\]  
The expression of '\(a\)' is written as:  
\[a=\left( \begin{array}{ccc}-\rho\omega^{2}&0&0\\ 0&-\rho\omega^{2}&0\\ 0&0&-\rho\omega^{2}\end{array} \right) \tag{4}\]  
where \(C_{ij}\) are the elastic moduli. As the 3D-FE model consists of modeling the propagation in the frequency domain, absorbing regions are used all around the plate to avoid undesired reflections. The absorbing regions (AR) in the 3D model are defined by a gradual increase of damping properties as explained in [45], where its formula is written as,  
\[C^{AR}_{ij}=C^{\prime}_{ij}\left(1+i*\left(\frac{Ra}{La}\right)^{3}\right)+i*C^ {\prime\prime}_{ij},\,C_{ij}=C^{\prime}_{ij}+iC^{\prime\prime}_{ij} \tag{5}\]  
Here, '\(Ra\)' and '\(La\)' are the interval and length of absorbing region respectively. The size of the AR used in this FE model is 1.5 \(\lambda_{MAX}\), where \(\lambda_{MAX}\) is the maximum wavelength for all modes existing in the plate in the studied frequency range. The full schematic of the 3D-FE model is shown in Figures 5 and 6.