# Thermoelastic Investigation of Carbon-Fiber-Reinforced Composites Using a Drop-Weight Impact Test
## 3 Finite Difference Method (MATLAB(r))
### Mathematical Model  
The mathematical model of the experiment is given by governing Equations (7) and (8) and the initial-boundary values as given in Equations (9) and (10):  
\[\partial_{t}^{2}\epsilon-c^{2}\Delta\epsilon=0,\ (t,x)\in[0,\infty]\times \mathbb{D}_{R}, \tag{7}\]  
\[\partial_{t}T-\kappa\Delta T=-\alpha T\partial_{t}\epsilon,\ (t,x)\in[0,\infty] \times\mathbb{D}_{R}, \tag{8}\]  
\[\epsilon\big{|}_{t=0}=\epsilon_{0},\ \epsilon\big{|}_{\partial\mathbb{D}_{R}}=0, \ \epsilon_{t}\big{|}_{0}=\eta, \tag{9}\]  
\[T\big{|}_{t=0}=T_{i}\quad T\big{|}_{\partial\mathbb{D}_{R}}=T_{0} \tag{10}\]  
where \(T\) is the temperature and \(\epsilon\) is the trace of the stress tensor. The parameters \(c,\kappa,\alpha\) are the wave speed, thermal conductivity, and thermal expansion factor, respectively. To reach an analytical solution, we introduce \(\tau=T-T_{0}\) and linearize \(T\partial_{t}\epsilon\approx T_{0}\partial_{t}\epsilon\). Then we solve Equation (11) with the boundary conditions given in Equations (12) and (13):  
\[\partial_{t}\tau-\kappa\Delta\tau=-\alpha T_{0}\partial_{t}\epsilon,\ (t,x)\in[0,\infty]\times\mathbb{D}_{R}, \tag{11}\]  
\[\epsilon\big{|}_{t=0}=\epsilon_{0},\ \epsilon\big{|}_{\partial\mathbb{D}_{R}}=0,\ \epsilon_{t}\big{|}_{t=0}=\eta_{0}, \tag{12}\]  
\[\tau\big{|}_{t=0}=\tau_{0}\quad T\big{|}_{\partial\mathbb{D}_{R}}=0 \tag{13}\]  
To solve the linearized version of the problem with axisymmetric initial and boundary data, we introduce the eigenvalues and eigenfunctions of the axisymmetric Laplace operator on the circular CFRP sheet with radius \(R\), \(\mathbb{D}_{R}\). The eigenfunctions for \(m=1,\cdots,\infty\) are given in Equation (14):  
\[u_{m}(r)=J_{0}\Big{(}\frac{\rho_{m}}{R}r\Big{)}, \tag{14}\]  
where \(J_{0}\) is the 0 th Bessel function of the first kind, and \(\rho_{m}\) is the \(m\)th zero of \(J_{0}\). To simplify the notation, we introduce \(\bar{\rho}_{m}=\frac{\rho_{m}}{R}\), the eigenvalue corresponding to the eigenfunctions \(u_{m}\) is \(-\beta_{m}^{2}\). The solution is given in Equations (15) and (16):\[\epsilon(t,r)=\ \sum_{n=0}^{\infty}\sum_{m=1}^{\infty}\cos(c\bar{\rho}_{m}t)c_{m}u _{m}(r)+\sin(c\bar{\rho}_{m}t)d_{m}u_{m}(r), \tag{15}\]  
\[\begin{array}{c}\tau(t,r)=\ \sum_{n=0}^{\infty}\sum_{m=1}^{\infty}e^{-\kappa \bar{\rho}_{m}^{2}t}\tau_{m}u_{m}(r)+\alpha T_{0}c\ \bar{\rho}_{m}\frac{\kappa\bar{\rho}_{m}(\cos c\bar{\rho}_{m}t-1)+c\sin c\bar{ \rho}_{m}t}{\kappa^{2}\bar{\rho}_{m}^{2}+c^{2}}\ c_{m}u_{m}(r)\\ \qquad\qquad-\alpha T_{0}c\bar{\rho}_{m}\frac{c\left(\cos c\bar{\rho}_{m}t-1 \right)-\kappa\bar{\rho}_{m}\sin c\bar{\rho}_{m}t}{\bar{\rho}_{m}(\kappa^{2} \bar{\rho}_{m}^{2}+c^{2}\ )}d_{m}u_{m}(r)\end{array} \tag{16}\]  
The coefficients are then computed as follows in Equation (17) through (19):  
\[c_{m}=\frac{1}{\frac{1}{2}R^{2}I_{1}(\rho_{m})^{2}}\int_{0}^{\mathbb{R}} \epsilon_{0}(r)u_{m}(r)r\ dr, \tag{17}\]  
\[d_{m}=\frac{1}{\frac{1}{2}R^{2}I_{1}(\rho_{m})^{2}}\int_{0}^{\mathbb{R}}\eta_ {0}(r)u_{m}(r)r\ dr, \tag{18}\]  
\[\tau_{m}=\frac{1}{\frac{1}{2}R^{2}I_{1}(\rho_{m})^{2}}\int_{0}^{\mathbb{R}} \tau_{0}(r)u_{m}(r)r\ dr, \tag{19}\]  
where \(I_{1}\) is the first Bessel function of the first kind.