# The Mesoscopic Numerical Simulation of GAP/CL20/AP Composite Solid Propellant Based on MPM and FEM
## 2 The Basic Theory of MPM  
MPM uses particles to represent the range of materials, and the particles contain historical variables such as stress, strain, position, speed, and energy. The interaction between particles and the updating of historical variables of particles depends on the background grid [28]. With the continuous development of the algorithm in recent decades [29; 30; 31; 32; 33; 34; 35; 36; 37; 38; 39; 40; 41; 42], the MPM has been widely used in numerical calculations of impact penetration, hypervelocity collisions, fluid-solid coupling, explosions, crack propagations, material failures, geotechnical impact failures, etc. [43; 44; 45; 46; 47; 48; 49; 50; 51; 52]. Particle position and other information can be easily mapped to the regular background grid, so it is very convenient to find the overlapping information of different parts by using its preprocessing technology.  
The control equations followed by the MPM include the mass conservation equation, the momentum conservation equation, and the energy conservation equation, which must be followed in the calculation process of the MPM.  
\[\left\{\begin{array}{l}\textit{Massconservation equation}:\ \rho(\textbf{X},t)J(\textbf{X},t)=\rho_{0}( \textbf{X})\\ \textit{Momentum conservation equation}:\ \sigma_{ij,j}+\rho b_{i}=\rho\ddot{u}_{i}\\ \textit{energy conservation equation}:\ \rho\ddot{w}^{\text{int}}=D_{ij} \sigma_{ij}\end{array}\right. \tag{1}\]  
where \(\rho(\textbf{X},t)\), \(J(\textbf{X},t)\), and \(\rho_{0}(\textbf{X})\) represent the density of the material at the current moment, the determinant of the deformation gradient, and the density at the initial moment, respectively. \(D_{ij}\) and \(\sigma_{ij}\) are the deformation rate tensor and the Cauchy stress tensor. \(u_{i}\) and \(b_{i}\) denote the displacement vector and physical force vector, and \(w^{\text{int}}\) is the internal energy of the material per unit mass.  
In the process of material point method calculation, the integration of the momentum equation and the calculation of a spatial derivative depend on the background grid, and its virtual work equation is as follows:  
\[\int_{V}\delta u_{i}\rho\ddot{u}_{i}dV+\int_{V}\delta u_{ij}\sigma_{ij}dV-\int _{V}\delta u_{i}\rho b_{i}dV-\int_{\Gamma_{t}}\delta u_{i}\ddot{\Gamma}_{i}d \Gamma=0 \tag{2}\]  
As shown in Figure 1, after a substance is dispersed into a series of particles, its density can be expressed as follows:  
\[\rho=\sum_{p=1}^{N_{p}}m_{p}\delta(x-x_{p}) \tag{3}\]  
where \(x_{ip}\) and \(m_{p}\) represented the mass and coordinates of the particles, respectively. \(\delta\) is the Dirac delta function. Formula (3) is substituted into the virtual work of Equation (2) and converted into the form of summation:  
\[\sum_{p=1}^{n_{p}}m_{p}\ddot{u}_{ip}\delta u_{ip}+\sum_{p=1}^{n_{p}}m_{p} \sigma_{ip}^{*}\delta u_{ip,j}-\sum_{p=1}^{n_{p}}m_{p}b_{ip}\delta u_{ip}-\sum _{p=1}^{n_{p}}m_{p}\overline{\dot{t}}_{ip}^{*}h^{-1}\delta u_{ip}=0 \tag{4}\]  
where \(u_{ip}=u_{i}(x_{p})\), \(\delta u_{ip,j}=\delta u_{i,j}(x_{p})\), \(\sigma_{ip}^{*}=\sigma_{ij}^{*}(x_{p})\), and \(b_{ip}=b_{i}(x_{p})\), \(\overline{\dot{t}}_{ip}^{*}=\overline{\dot{t}}_{i}^{*}(x_{p})\). \(h\) is the thickness of the assumed boundary, the integral term in Formula (2) is converted into a summation term by Formula (4), and the motion equation of the background grid node is as follows:  
\[\dot{p}_{iI}=\dot{f}_{iI}^{\text{int}}+f_{iI}^{\text{ext}},x_{I}\
otin\Gamma_{u} \tag{5}\]  
\[\dot{p}_{iI}=M_{I}\dot{u}_{iI} \tag{6}\]  
\[M_{IJ}=\sum_{p=1}^{n_{p}}m_{p}N_{Ip}N_{Ip} \tag{7}\]  
where \(\dot{p}_{iI}\) represents the momentum of grid nodes \(I\) in the direction of \(i\) and \(M_{IJ}\) is the mass matrix of the background grid.  
Figure 1: Schematic diagram of material point method.