# Material Point Simulation Method for Concrete Medium Fracture and Fragmentation under Blast Loading
## 2 Basic Theory of the MPM  
The MPM discretizes the continuum into a series of material points; each material point represents a material region and carries all the material information of that region, such as mass, velocity, stress, and strain. The collection of all material points represents the entire material area. Therefore, the density of the continuum can be approximated as:  
\[\rho(x_{i})=\sum_{p=1}^{n_{p}}m_{p}\delta\big{(}x_{i}-x_{ip}\big{)} \tag{1}\]  
In the equation, \(n_{p}\) represents the total number of material points, \(m_{p}\) is the mass of the material point, \(\delta\) is the Dirac Delta function, and \(x_{ip}\) represents the coordinates of the material point.  
Similar to the FEM, the solution format of the MPM is also based on the weak form. The solution format of the MPM is as follows:  
\[\int\limits_{\Omega}\rho\ddot{u}_{i}\dot{g}u_{i}dV+\int\limits_{\Omega}\rho \sigma_{ij}^{s}\delta u_{i,j}dV-\int\limits_{\Omega}\rho b_{i}\dot{o}u_{i}dV- \int\limits_{\Gamma_{\Gamma_{\Gamma_{\Gamma_{\Gamma}}}}}\rho\overline{t}_{i} ^{s}\delta u_{i}dA=0 \tag{2}\]  
In the equation, \(\rho\) is the density, \(\ddot{u}_{i}\) is the acceleration, \(b_{i}\) is the force acting on the unit mass of the object, \(\Omega\) is the solving area, \(\Gamma_{t}\) is the given surface force boundary, \(\sigma_{ij}^{s}=\sigma_{ij}/\rho\) is the comparing stress, and \(\overline{t}_{i}^{s}=\overline{t}_{i}/\rho\) is the comparing boundary surface force.  
Substituting Equation (1) into Equation (2), it can be obtained:  
\[\sum_{p=1}^{n_{p}}m_{p}\ddot{u}_{ip}\delta u_{ip}+\sum_{p=1}^{n_{p}}m_{p} \sigma_{ij}^{s}\delta u_{ip,j}-\sum_{p=1}^{n_{p}}m_{p}b_{ip}\delta u_{ip}-\sum _{p=1}^{n_{p}}m_{p}\overline{t}_{ip}^{s}h^{-1}\delta u_{ip}=0 \tag{3}\]  
In the equation, \(u_{ip}=u_{i}\big{(}x_{p}\big{)}\), \(\delta u_{ip,j}=\delta u_{i,j}\big{(}x_{p}\big{)}\), \(\sigma_{ip}^{s}=\sigma_{ij}^{s}\big{(}x_{p}\big{)}\), \(b_{ip}=b_{i}\big{(}x_{p}\big{)}\), \(\overline{t}_{ip}^{s}=\overline{t}_{i}^{s}\big{(}x_{p}\big{)}\), \(h\) is used to introduce a hypothetical boundary layer thickness in order to convert the final term of the left-hand side boundary integral in Equation (2) into a volume integral. As shown in Equation (3), the material point method converts the integrals in Equation (2) into the product of the value of the integrand at each material point and the volume represented by the material point, which is a form of point integration.  
The MPM arranges a background grid for the material domain and approximately solves the momentum equation using the virtual work principle. Within one computational time step, the mass and momentum of the material points are first mapped to the nodes of the background grid using shape functions to calculate the mass and momentum of the background grid nodes. Furthermore, the momentum equation is integrated at the background grid nodes, and the change in velocity and position of the background grid nodes is mapped back to the material points using shape functions to update their velocity and position. Finally, the deformed background grid is discarded, and the regular background grid is used again in the next time step. The solution process of the MPM within one computational time step is shown in Figure 1, and the MPM calculation algorithm is shown in Figure 2.  
Figure 1: Solution process of the MPM within one computational time step.  
Figure 2: MPM calculation algorithm.