## 2 Development of Full-Process Numerical Analyzing Method
### Development of the Progressive Low-Velocity Impact Damage Prediction Method  
#### 2.1.1 Transient Stress Analysis  
Figure 1 shows the transient analysis model of laminate and punch. At the \(t\) moment, the equilibrium equation can be deduced as:  
\[\sigma_{ij,j}^{t}=\rho^{t}\ddot{u}_{i}+\mu^{t}\dot{u}_{i} \tag{1}\]  
where, \(\rho^{t}\) and \(\mu^{t}\) are the density and dynamic friction coefficient of laminate at \(t\) moment, respectively. \(\dot{u}_{i}\) and \(\ddot{u}_{i}\) are the first and second derivative of the displacement to the time respectively.  
At this moment, boundary conditions of the load are:  
\[(\sigma_{ij}^{t})n_{j}-\overline{T}_{i}^{t}=0\ \ (\text{on}\ S_{\sigma}) \tag{2}\]  
where, \(S_{\sigma}\) represents the stress boundary.  
The equivalent integration of the equilibrium equation and the load boundary condition can be expressed, as follows:  
\[\int\limits_{V}\delta u_{i}\Big{(}\sigma_{ij,j}^{t}-\rho^{t}\ddot{u}_{i}-\mu^{ t}\dot{u}_{i}\Big{)}dV-\int\limits_{S\sigma}\delta u_{i}\Big{(}(\sigma_{ij}^{t})n_ {j}-\overline{T}_{i}^{t}\Big{)}dS=0 \tag{3}\]  
As there will be geometric nonlinearity during the deformation of composites under the low-velocity impact load, strain tensor at \(t\) moment can be expressed, as follows:  
\[\varepsilon_{ij}^{t}=B_{ij}U_{j}^{t}=\frac{1}{2}\Big{(}u_{i,j}^{t}+u_{j,i}^{t} +u_{k,i}^{t}u_{k,j}^{t}\Big{)}\ \ \ \ \ (i,j,k=x,y,z) \tag{4}\]  
Decompose the above equation into linear and nonlinear terms  
\[\varepsilon=\varepsilon_{L}+\varepsilon_{NL} \tag{5}\]  
There is the following relationship of \(\sigma_{ij}^{t}\) and \(\varepsilon_{kl}^{t}\) at \(t\) moment.  
\[\sigma_{ij}^{t}=\overline{Q}_{ijkl}^{t-\Delta t}\varepsilon_{kl}^{t} \tag{6}\]  
Figure 1: Transient analysis model of laminate and punch.  
where \(\overline{Q}_{ijkl}^{t-\Delta t}\) represents material elasticity matrix at \(t-\Delta t\) moment, and it can be obtained by coordinate transformation.  
\[\left[\overline{Q}_{ijkl}^{t-\Delta t}\right]=\left[\mathrm{T}\right]\left[Q^{t -\Delta t}\right]\left[\mathrm{T}\right]^{T} \tag{7}\]  
Figure 2 shows the included angle \(\theta\) between the direction of fiber and the global coordinate system, the transition matrix \(\left[\mathrm{T}\right]\) can be expressed, as follows:  
\[\mathbf{T}=\left[\begin{array}{cccccc}\cos^{2}\theta&\sin^{2}\theta&0&0&0&- 2\sin\theta\cos\theta\\ \sin^{2}\theta&\cos^{2}\theta&0&0&0&2\sin\theta\cos\theta\\ 0&0&1&0&0&0\\ 0&0&0&\cos\theta&\sin\theta&0\\ 0&0&0&-\sin\theta&\cos\theta&0\\ \sin\theta\cos\theta&-\sin\theta\cos\theta&0&0&0&\cos^{2}\theta-\sin^{2} \theta\end{array}\right] \tag{8}\]  
According to Equations (3) and (6), the stress equilibrium equation at each moment can finally be deduced as:  
\[\int_{V}\left(\delta\varepsilon_{ij}\overline{Q}_{ijkl}^{t-\Delta t} \varepsilon_{kl}+\delta u_{i}\rho^{t}\ddot{u}_{i}+\delta u_{i}\mu^{t}\dot{u}_{ i}\right)\mathrm{d}V+\int_{V_{n-1}}\sigma_{ij}^{n-1}\delta(\Delta\eta_{ij}) \mathrm{d}V=\int_{S_{\sigma}}\overline{T}_{i}^{n}\delta u_{i}\mathrm{d}S \tag{9}\]  
where, \(\overline{T}_{i}^{t}\) and \(\overline{T}_{i}^{n}\) are the surface force at \(t\) moment and \(n\)-th step in numerical analysis respectively; \(\delta\varepsilon_{ij}\) represents the strain at \(t\) moment; and, \(\Delta\eta_{ij}\) is the nonlinear term of strain increment.  
Equation (9) can be solved by finite element method (FEM). It should be noted that when damage occurs, the stiffness matrix of the damaged area should be degraded and updated; the stress and strain fields will redistribute and need to be recalculated with an updated stiffness matrix.