Stochastic Material Point Method for Analysis in Non-Linear Dynamics of Metals

Abstract

A stochastic material point method is proposed for stochastic analysis in non-linear dynamics of metals with varying random material properties. The basic random variables are parameters of equation of state and those of constitutive equation. In conjunction with the material point method, the Taylor series expansion is employed to predict first- and second-moment characteristics of structural response. Unlike the traditional grid methods, the stochastic material point method does not require structured mesh; instead, only a scattered cluster of nodes is required in the computational domain. In addition, there is no need for fixed connectivity between nodes. Hence, the stochastic material point method is more suitable than the stochastic method based on grids, when solving dynamics problems of metals involving large deformations and strong nonlinearity. Numerical examples show good agreement between the results of the stochastic material point method and Monte Carlo simulation. This study examines the accuracy and convergence of the stochastic material point method. The stochastic material point method offers a new option when solving stochastic dynamics problems of metals involving large deformation and strong nonlinearity, since the method is convenient and efficient.

1. Introduction

The traditional grid methods have been successfully employed for several decades in engineering practice but has limitations when the considered non-linear dynamics problem of metals involves large deformation, such as hyper-velocity impact and explosion, due to distortions in the mesh when implementing a Lagrangian approach. The deformation of the grids can lead to numerical inaccuracies or even render the calculation impossible [1,2]. To overcome these difficulties, a great number of meshless methods—such as the smoothed particle hydrodynamics method [3,4], the diffuse element method [5], the element-free Galerkin method [6,7], partition of unity [8], the reproducing kernel particle method [9,10], the cracking particles [11], the dual-horizon peridynamics [12,13] and the material point method (MPM) [14,15,16]—have been proposed and developed in recent years to solve non-linear dynamics problems of metals [17,18].

Among these methods, MPM shows many advantages in tension stability and efficiency [19,20,21,22]. Compared with the traditional grid methods, MPM is specifically intriguing because of its simplicity, which have proven useful for solving solid mechanics problems involving large deformation and strong nonlinearity, such as hyper-velocity impact and explosion. [23,24,25,26,27]. MPM can be classified as a meshless method, inheriting the advantages of both Eulerian and Lagrangian methods. Hence, MPM avoids the mesh distortion in the Lagrangian methods and the convection problems in the Eulerian methods [28,29]. To date, MPM and its extensions have been applied successfully to explosion problems [30], impact and penetration problems [1,31,32], fluid-solid interaction problems [33,34] and so on. However, all of the developments in MPM have so far focused on deterministic problems. Probabilistic models using MPM have not received much attention. Furthermore, the stochastic analysis of non-linear dynamics problem involves large deformation and strong nonlinearity, such as hyper-velocity impact and explosion; depends mainly on Monte Carlo simulation, which has disadvantages of low calculating efficiency; and is not naturally suitable for complex engineering practice [35]. Hence, high-efficiency stochastic analysis involving MPM provides a rich, relatively unexplored field for future research [36].

This paper proposes a stochastic material point method (SMPM) for stochastic analysis in non-linear dynamics metal structure with varying random material properties. The samples of basic random variables are parameters of equation of state and those of constitutive equation. In conjunction with the MPM, the Taylor series expansion is employed to predict first- and second-moment characteristics (mean and variance) of structural response. Two numerical examples based on non-linear problems of steel are carried out to examine the accuracy and convergence of the proposed method.

2. Stochastic Material Point Method

The main idea of the SMPM is to introduce the stochastic theory into MPM to make stochastic analysis of non-linear dynamics. Thus, the SMPM has similar features to MPM. In the SMPM, the continuum body is divided into a set of $N_p$ infinitesimal mass elements and these infinitesimal mass elements are represented by a finite collection of material particles, as shown in Figure 1.

The governing equations of continuous materials are standard conservation equations for mass and momentum

$$\frac{d\rho }{dt}+\rho \mathsf{\nabla }\cdot v=0$$

(1)

$$\rho \frac{dv}{dt}=\mathsf{\nabla }\cdot \mathsf{\sigma }+\rho \mathit{b}$$

(2)

where ρ is the mass density, $v$ is the velocity, σ is the Cauchy stress tensor and b is the body force per unit mass.

The key steps for SMPM, of which includes the coupling between MPM module and stochastic module and the return of the needed response quantities are illuminated in Figure 2. The following detailed steps are performed at each time step increment Δt:

The first step is to input initial physical variables and at the same time, random variables are passed to the MPM module from the stochastic module.

The second step is mapping velocity and mass of the particles to the background grid nodes with the shape functions. The concrete forms are

$$v_i^t={\displaystyle \sum _{p=1}^{N_p}{v}_p^t}{N}_i\left({\mathit{x}}_p^t\right).$$

(3)

$$m_i^t={\displaystyle \sum _{p=1}^{N_p}{m}_p}{N}_i\left({\mathit{x}}_p^t\right)$$

(4)

where ${\mathit{x}}_p^t$ is the position vector of the particle at time $t$, $v_i^t$ is the nodal velocity at time $t$, $v_p^t$ is the particle velocity at time $t$, $N_i$ is the shape function for node $i$, $m_i^t$ is the nodal mass at time $t$ and $m_p$ is the particle mass.

The third step is calculating the nodal forces, the internal force vector ${\mathit{f}}_i^{\mathrm{int},\mathrm{t}}$ and the external force vector $f_i^{ext,t}$. The concrete forms are

$$f_i^{\mathrm{int},\mathrm{t}}=-{\displaystyle \sum _{p=1}^{N_p}\frac{m_p}{{\rho }_p^t}{\mathsf{\sigma }}_p^t\mathsf{\nabla }{N}_i\left({\mathit{x}}_p^t\right).}$$

(5)

$${\mathit{f}}_i^{ext,t}=-{\displaystyle \sum _{p=1}^{N_p}{m}_p\mathsf{\nabla }{N}_i\left({\mathit{x}}_p^t\right)}{\mathit{b}}_p^t+{\displaystyle {\int }_{{\mathsf{\Gamma }}_{\tau }}\overline{\mathsf{\tau }}}\cdot N_i\left({\mathit{x}}_p^t\right)\mathrm{ds}$$

(6)

where $\overline{\mathsf{\tau }}$ is the vectors of prescribed surface traction, ${\mathsf{\Gamma }}_{\tau }$ is the stress boundary condition, ${\rho }_p^t$.is the density of the material particle at time $t$, ${\mathsf{\sigma }}_p^t$ is the stress tensor of the material particle at time $t$ and ${\mathit{b}}_p^t$ is the body force of the material particle at time $t$.

The forth step is solving the momentum equation in the background grid nodes. The concrete form is

$${\mathit{P}}_i^{t+1}={\mathit{P}}_i^t+\left({\mathit{f}}_i^{ext,t}+{\mathit{f}}_i^{\mathrm{int},\mathrm{t}}\right)\mathsf{\Delta }t$$

(7)

where ${\mathit{P}}_i^t$ is the momentum of the node at time t and ${\mathit{P}}_i^{t+1}$.is the momentum of the node at time $t+1$.

The fifth step is to map the information from the background grid nodes back to the material particles. The velocity and displacement vectors of the material particles are mapped from the background grid nodes.

For history-dependent materials, it is convenient to carry strain and stress as well as history variables along with the material particles. Hence, the sixth step is to update the state of the material particles. The information of strain and stress can be updated with the constitutive equation and equation of state and the variables of strain, stress and pressure will be discussed in Section 3. Meanwhile, the first-order and second-order derivative of the response quantities are calculated and passed to the stochastic module from the MPM module. Then, the mean and variance of the response quantities are calculated in the stochastic module and the content of random analysis will be presented in Section 4.

The seventh step is to update the feature size and density of the material particles. The concrete forms are

$$l_{pj}^{t+1}=l_{pj}^t\left(1+\mathsf{\Delta }{\epsilon }_{jj}\right), (j=x,y,z)$$

(8)

$${\rho }_p^{k+1}=\frac{1}{\left(1+\mathsf{\Delta }{\mathsf{\epsilon }}_{xx}+\mathsf{\Delta }{\mathsf{\epsilon }}_{xx}+\mathsf{\Delta }{\mathsf{\epsilon }}_{yy}+\mathsf{\Delta }{\mathsf{\epsilon }}_{zz}\right)}$$

(9)

where $l_{pj}^t$ is feature size of the material particle at time $t$ in the direction $j$, $l_{pj}^{t+1}$ is feature size of the material particle at time $t+1$ in the direction $j$, $\mathsf{\Delta }{\mathsf{\epsilon }}_{jj}$ is the strain increment in the direction $j$ and ${\rho }_p^{k+1}$ is the density of the material particle at time $t+1$.

The eighth step is resetting the background grid elements to the original undeformed state.

This completes the computational cycle of a time step. Repeat steps 2–8 until the time has advanced to the desired value.

3. Material Property Equations

3.1. Equation of State

One century ago, Mie [37] and Grüneisen [38] developed a theory for metals which concluded that the pressure could be considered as a linear function of internal energy. The name Mie-Grüneisen is now associated to an equation of state which verifies this assumption. The Mie-Grüneisen equation of state is conveniently couched in an enthalpy framework, which contains two states of the materials defined by the index μ. And μ > 0 represents the compressed state, μ < 0 represents the expanded state. The concrete form is

$$\{\begin{array}{cc}{p}^{*}=\frac{{\rho }_0{C}^2\mathsf{\mu }\left[1+\left(1-\frac{{\mathsf{\gamma }}_0}{2}\right)\mathsf{\mu }-\frac{a}{2}{\mathsf{\mu }}^2\right]}{{\left[1-\left(S_1-1\right)\mathsf{\mu }-S_2\frac{{\mathsf{\mu }}^2}{\mathsf{\mu }+1}-S_3\frac{{\mathsf{\mu }}^3}{{\left(\mathsf{\mu }+1\right)}^2}\right]}^2}+\left({\gamma }_0+a\mathsf{\mu }\right)e\hfill & \hfill \left(\mathsf{\mu }>0\right)\\ p^{*}={\rho }_0{C}^2\mathsf{\mu }+\left({\mathsf{\gamma }}_0+a\mathsf{\mu }\right)e\hfill & \hfill \left(\mathsf{\mu }<0\right)\end{array}$$

(10)

where $p^{*}$ is the pressure of a metal, S₁, S₂, S₃, γ₀ and a are constants, ρ₀ is the original value of the density and e is the internal energy of a metal. The index μ is defined as $\mathsf{\mu }=\frac{\rho }{{\rho }_0}-1$.

3.2. Constitutive Equation

In general, the response of metals under high-speed impact conditions involves consideration of effects of strain, strain rate and temperature. The Johnson-Cook plasticity model [39,40,41] is employed to model the flow stress behavior of ductile materials. The Johnson-Cook model represents the Von Mises flow stress σ_y as a function of some parameters as follows

$${\mathsf{\sigma }}_y=\left(A+B{\mathsf{\epsilon }}_P^n\right)\left(1+C\ln {\dot{\mathsf{\epsilon }}}_P^{*}\right)\left(1-T^{*m}\right)$$

(11)

where ε_P is equivalent plastic strain, ${\dot{\mathsf{\epsilon }}}_{\mathrm{P}}^{*}$ is equivalent plastic strain rate, A, B, C and m are constants, n is strain hardening exponent and $T^{*}$ is homologous temperature defined as

$$T^{*}=\frac{T-T_{\mathrm{room}}}{T_{\mathrm{melt}}-T_{\mathrm{room}}}$$

(12)

where T is the material temperature, T_melt is the melting temperature and T_room is the room temperature.

3.3. Von Mises Elastic-Plastic Material Model

Von Mises yield model can be used to describe the plastic behavior of metal materials [42]. The Von Mises elastic-plastic model can be used to describe deviatoric stress update procedure and the radial return mapping can be used at the same time. Under the assumption that the metal materials are in elastic state, the trial deviatoric stress tensor ${}^{*}{S}_{ij}^{n+1}$ can be calculated by

$${}^{*}{S}_{ij}^{n+1}=S_{ij}^{R^n}+2G\mathsf{\Delta }{{\mathsf{\epsilon }}^{\prime }}_{ij}$$

(13)

where $\mathsf{\Delta }{{\mathsf{\epsilon }}^{\prime }}_{ij}$ is increment of the deviatoric strain, $S_{ij}^{R^n}$ is the deviatoric stress tensor after the rotation, the superscript n is the previous iteration step at the same time step and the superscript n + 1 is the current iteration step at the same time step.

At this moment, the trial Von Mises flow stress ${}^{*}{\mathsf{\sigma }}_{\mathrm{eq}}^{n+1}$ can be calculated by

$${}^{*}{\mathsf{\sigma }}_{\mathrm{eq}}^{n+1}={\left(\frac{3}{2}{}^{*}{S}_{ij}^{n+1}{}^{*}{S}_{ij}^{n+1}\right)}^{\frac{1}{2}}$$

(14)

If the following condition holds ${}^{*}{\mathsf{\sigma }}_{\mathrm{eq}}^{n+1}<{\mathsf{\sigma }}_y^n$, the deviatoric stress is in the yield surface. And the trial deviatoric stress tensor acts as the deviatoric stress tensor $S_{ij}^{n+1}={}^{*}{S}_{ij}^{n+1}$.

If the following condition holds ${}^{*}{\mathsf{\sigma }}_{\mathrm{eq}}^{n+1}>{\mathsf{\sigma }}_y^n$, the deviatoric stress is beyond the yield surface. Therefore, ${}^{*}{S}_{ij}^{n+1}$ should be reduced scale representation with the purpose of making it in the yield surface, the concrete form as follows

$$S_{ij}^{n+1}=m{}^{*}{S}_{ij}^{n+1}$$

(15)

where the coefficient m can be calculated as

$$m=\frac{{\mathsf{\sigma }}_y^{n+1}}{{}^{*}{\mathsf{\sigma }}_{\mathrm{eq}}^{n+1}}$$

(16)

Considering Equations (13) and (15), the plastic strain increment $\mathsf{\Delta }{\mathsf{\epsilon }}_{ij}^p$ can be calculated as

$$\mathsf{\Delta }{\mathsf{\epsilon }}_{ij}^p=\frac{{}^{*}{S}_{ij}^{n+1}-S_{ij}^{n+1}}{2G}=\frac{1-m}{2G}$$

(17)

Considering Equations (14) and (16), the equivalent plastic strain increment is

$$\mathsf{\Delta }{\mathsf{\epsilon }}_p={\left(\frac{2}{3}\mathsf{\Delta }{\mathsf{\epsilon }}_{ij}^p\mathsf{\Delta }{\mathsf{\epsilon }}_{ij}^p\right)}^{\frac{1}{2}}=\frac{{}^{*}{\mathsf{\sigma }}_{\mathrm{eq}}^{n+1}-{\mathsf{\sigma }}_y^{n+1}}{3G}$$

(18)

The flow stress of Von Mises isotropic hardening materials can be updated by

$${\mathsf{\sigma }}_y^{n+1}={\mathsf{\sigma }}_y^n+H\mathsf{\Delta }{\mathsf{\epsilon }}_p$$

(19)

where H is the plastic modulus.

By substituting of Equations (17) and (20) into Equation (18), the equivalent plastic strain increment can be calculated by

$$\mathsf{\Delta }{\mathsf{\epsilon }}_p=\frac{{}^{*}{\mathsf{\sigma }}_{\mathrm{eq}}^{n+1}-{\mathsf{\sigma }}_y^n}{3G+H}$$

(20)

In a word, the calculation procedures of Von Mises elastic-plastic material model are summarized as follows

• Calculate the trial deviatoric stress tensor ${}^{*}{S}_{ij}^{n+1}$ and the trial Von Mises flow stress ${}^{*}{\mathsf{\sigma }}_{\mathrm{eq}}^{n+1}$.
• If the condition ${}^{*}{\mathsf{\sigma }}_{\mathrm{eq}}^{n+1}>{\mathsf{\sigma }}_y^n$ holds, the equivalent plastic strain increment can be calculated by Equation (20). Otherwise, the materials have no plastic deformation.
• Update equivalent plastic strain by ${\mathsf{\epsilon }}_p^{n+1}={\mathsf{\epsilon }}_p^n+\mathsf{\Delta }{\mathsf{\epsilon }}_p$.
• Update flow stress by Equation (19), get ${\mathsf{\sigma }}_y^{n+1}$.
• Calculate the coefficient m with the new flow stress ${\mathsf{\sigma }}_y^{n+1}$ by Equation (16) and update the deviatoric stress which satisfies the condition in the yield surface with the radial return mapping.

4. Random Method

4.1. The Basic Random Variables and Random Response

To introduce random parameters to describe the random characteristics of materials is one of the most important factors in stochastic analysis. The random parameters include parameters of equation of state, constitutive equation and those introduced to establish stochastic modeling by defining basic random variables $X={\left[x_1\dots ,x_k,\dots x_n\right]}^T$. Therefore, information for random variables must be provided, such as the mean $E\left(x_k\right)={\mathsf{\mu }}_k\left(k=1,2\dots n\right)$, the variance $D\left(x_k\right)={\mathsf{\sigma }}_k\left(k=1,2\dots n\right)$ and the correlation coefficient between two random variables ${\rho }_{kj}$.

The direct solution of structural response is challenging, as it comprises multidimensional integral and the structural response is denoted by $Y=f\left(x_1\dots ,x_k,\dots x_n\right)$. A possible solution to overcome this issue is approximating the structural response via an expression that is an explicit function of the basic random variables. Consequently, introducing the powerful technique of expansion is necessary to calculate first- and second-moment characteristics (mean and variance) of structural response [43], such as, Neumann series expansion [44], polynomial chaos expansion [45] and Taylor series expansion [46,47,48]. It is the first attempt to apply MPM in the field of stochastic analysis and this study is focused on applying the new method in stochastic mechanics and thus, Taylor series expansion has been employed because of its simple and convenient.

In order to calculate the approximation for first-moment of the structural response, assume $Y$ is expressed in terms of its corresponding second order Taylor series expansion about the point $\overline{X}={\left[{\mathsf{\mu }}_1\dots ,{\mathsf{\mu }}_k,\dots {\mathsf{\mu }}_n\right]}^T$. The concrete form is calculated as follows

$$Y=f\left({\mathsf{\mu }}_1\dots ,{\mathsf{\mu }}_k,\dots {\mathsf{\mu }}_n\right)+{\displaystyle \sum _{k=1}^n{\frac{\partial Y}{\partial x_k}|}_{\overline{X}}\left(x_k-{\mathsf{\mu }}_k\right)}+\frac{1}{2}{\displaystyle \sum _{k=1}^n{\displaystyle \sum _{j=1}^n{\frac{{\partial }^{2}Y}{\partial x_k\partial x_j}|}_{\overline{X}}\left(x_k-{\mathsf{\mu }}_k\right)}}\left(x_j-{\mathsf{\mu }}_j\right)+R_2\left(x_1\dots ,x_k,\dots x_n\right)$$

(21)

where $R_2\left(x_1\dots ,x_k,\dots x_n\right)$ is the remainder of second order Taylor series expansion.

Magnitude analysis proved that the remainder has a relatively weak intensity compared to the first-three terms of Taylor series expansion and hence, the influence of the remainder on deformation of the structure can be neglected. Calculating the mean of both sides of Equation (21) gives the mean of the structural response as follows

$$E\left[Y\right]=f\left({\mathsf{\mu }}_1\dots ,{\mathsf{\mu }}_k,\dots {\mathsf{\mu }}_n\right)+\frac{1}{2}{\displaystyle \sum _{k=1}^n{\displaystyle \sum _{j=1}^n{\frac{{\partial }^{2}Y}{\partial x_k\partial x_j}|}_{\overline{X}}{\mathsf{\sigma }}_k{\mathsf{\sigma }}_j}}{\rho }_{kj}$$

(22)

In a similar way, in order to calculated the approximation for the second-moment characteristics of structural response, assume $Y$ is expressed in terms of its corresponding first order of the Taylor series expansion about the point $\overline{X}={\left[{\mathsf{\mu }}_1\dots ,{\mathsf{\mu }}_k,\dots {\mathsf{\mu }}_n\right]}^T$. The concrete form is calculated as follows

$$Y=f\left({\mathsf{\mu }}_1\dots ,{\mathsf{\mu }}_k,\dots {\mathsf{\mu }}_n\right)+{\displaystyle \sum _{k=1}^n{\frac{\partial Y}{\partial x_k}|}_{\overline{X}}\left(x_k-{\mathsf{\mu }}_k\right)}+R_1\left(x_1\dots ,x_k,\dots x_n\right)$$

(23)

where $R_1\left(x_1\dots ,x_k,\dots x_n\right)$ is the remainder of first order Taylor series expansion.

Similarly, the variance of the structural response can be calculated by

$$D\left[Y\right]={\displaystyle \sum _{k=1}^n{{\displaystyle \sum _{j=1}^n\frac{\partial Y}{\partial x_k}}|}_{\overline{X}}{\frac{\partial Y}{\partial x_j}|}_{\overline{X}}}{\mathsf{\sigma }}_k{\mathsf{\sigma }}_j{\rho }_{kj}$$

(24)

The Taylor series expansion is applied to stochastic analysis while meeting certain conditions. If the variation coefficient is less than about 0.2 or 0.3, a better approximation can be obtained by using the first-two terms of the Taylor series expansion. Further, the condition is widened in approximation by using the first-three terms of the Taylor series expansion [49]. The first-moment characteristics of structural response can be obtained through the first-three term of the Taylor series expansion; correspondingly, the second-moment characteristics can be obtained from the first-two terms of the Taylor series expansion. Additionally, probabilistic convergence studies of second order stochastic perturbation technique show that the second order method cannot be used for the coefficient of variation larger than 0.15 to determine efficiently higher order statistics [43,50,51,52,53]. To avoid the obvious errors, a variation coefficient of no larger than 0.15 is selected for this paper. Meanwhile, through comparison with the Monte Carlo method, the reasonableness and correctness of the Taylor series expansion described in this paper will be proved through two concrete instances in Section 5.

Substituting concrete forms into Equation (22) and (24), then we can reach the stochastic structural response after rearranging.

Considering Equation (10) in Section 3.1, the mean and variance of the pressure can be calculated by

$$E\left[p^{*}\right]=f\left({\mathsf{\mu }}_1\dots ,{\mathsf{\mu }}_k,\dots {\mathsf{\mu }}_n\right)+\frac{1}{2}{\displaystyle \sum _{k=1}^n{\displaystyle \sum _{j=1}^n{\frac{{\partial }^2{p}^{*}}{\partial x_k\partial x_j}|}_{\overline{X}}{\mathsf{\sigma }}_k{\mathsf{\sigma }}_j}}{\rho }_{kj}$$

(25)

$$D\left[p^{*}\right]={\displaystyle \sum _{k=1}^n{{\displaystyle \sum _{j=1}^n\frac{\partial p^{*}}{\partial x_k}}|}_{\overline{X}}{\frac{\partial p^{*}}{\partial x_j}|}_{\overline{X}}}{\mathsf{\sigma }}_k{\mathsf{\sigma }}_j{\rho }_{kj}$$

(26)

Considering Equation (11) in Section 3.2 and the relationship between σ_y and $S_{ij}$ in Section 3.3, the mean and variance of the equivalent plastic stress and the deviatoric stress can be calculated by

$$E\left[{\mathsf{\sigma }}_y^{}\right]=f\left({\mathsf{\mu }}_1\dots ,{\mathsf{\mu }}_k,\dots {\mathsf{\mu }}_n\right)+\frac{1}{2}{\displaystyle \sum _{k=1}^n{\displaystyle \sum _{j=1}^n{\frac{{\partial }^2{\mathsf{\sigma }}_y^{}}{\partial x_k\partial x_j}|}_{\overline{X}}{\mathsf{\sigma }}_k{\mathsf{\sigma }}_j}}{\rho }_{kj}$$

(27)

$$D\left[{\mathsf{\sigma }}_y\right]={\displaystyle \sum _{k=1}^n{{\displaystyle \sum _{j=1}^n\frac{\partial {\mathsf{\sigma }}_y^{}}{\partial x_k}}|}_{\overline{X}}{\frac{\partial {\mathsf{\sigma }}_y^{}}{\partial x_j}|}_{\overline{X}}}{\mathsf{\sigma }}_k{\mathsf{\sigma }}_j{\rho }_{kj}$$

(28)

$$E\left[S_{ij}\right]=\frac{{}^{*}{S}_{ij}^{n+1}}{\sqrt{\frac{3}{2}{}^{*}{S}_{ij}^{n+1}{}^{*}{S}_{ij}^{n+1}}}\frac{3G}{3G+H}E\left[{\mathsf{\sigma }}_y^{}\right]+\frac{H}{3G+H}{}^{*}{S}_{ij}^{n+1}$$

(29)

$$D\left[S_{ij}\right]=\frac{2}{3}{\left(\frac{3G}{3G+H}\right)}^{2}D\left[{\mathsf{\sigma }}_y^{}\right]$$

(30)

where $E[\cdot ]$ is the operation of the mean, $D[\cdot ]$ is the operation of the variance.

4.2. The First-Order and Second-Order Derivatives of the Response Quantities

The first- and second-order derivatives of the response quantities with respect to the random input parameters are involved in Equations (25–28). The estimation of these derivatives is a key problem in stochastic analysis, which is discussed as follows. Suppose the parameters of Mie-Grüneisen equation of state model and those of Johnson-Cook plasticity model are random variables, such as S₁, S₂, S₃, γ₀, a, A, B, C, m and n.

If μ > 0, the derivatives in Equations (25) and (26) can be calculated by

$$\left[\begin{array}{ccccc}\frac{\partial p^{*}}{\partial S_1}& \frac{\partial p^{*}}{\partial S_2}& \frac{\partial p^{*}}{\partial S_3}& \frac{\partial p^{*}}{\partial {\mathsf{\gamma }}_0}& \frac{\partial p^{*}}{\partial a}\\ \frac{{\partial }^2{p}^{*}}{\partial S_1^2}& \frac{{\partial }^2{p}^{*}}{\partial S_1\partial S_2}& \frac{{\partial }^2{p}^{*}}{\partial S_1\partial S_3}& \frac{{\partial }^2{p}^{*}}{\partial S_1\partial {\mathsf{\gamma }}_0}& \frac{{\partial }^2{p}^{*}}{\partial S_1\partial a}\\ \frac{{\partial }^2{p}^{*}}{\partial S_2\partial S_1}& \frac{{\partial }^2{p}^{*}}{\partial S_2^2}& \frac{{\partial }^2{p}^{*}}{\partial S_2\partial S_3}& \frac{{\partial }^2{p}^{*}}{\partial S_2\partial {\mathsf{\gamma }}_0}& \frac{{\partial }^2{p}^{*}}{\partial S_2\partial a}\\ \frac{{\partial }^2{p}^{*}}{\partial S_3\partial S_1}& \frac{{\partial }^2{p}^{*}}{\partial S_3\partial S_2}& \frac{{\partial }^2{p}^{*}}{\partial S_3^2}& \frac{{\partial }^2{p}^{*}}{\partial S_3\partial {\mathsf{\gamma }}_0}& \frac{{\partial }^2{p}^{*}}{\partial S_3\partial a}\\ \frac{{\partial }^2{p}^{*}}{\partial {\mathsf{\gamma }}_0\partial S_1}& \frac{{\partial }^2{p}^{*}}{\partial {\mathsf{\gamma }}_0\partial S_2}& \frac{{\partial }^2{p}^{*}}{\partial {\mathsf{\gamma }}_0\partial S_3}& \frac{{\partial }^2{p}^{*}}{\partial {\mathsf{\gamma }}_0^2}& \frac{{\partial }^2{p}^{*}}{\partial {\mathsf{\gamma }}_0\partial a}\\ \frac{{\partial }^2{p}^{*}}{\partial a\partial S_1}& \frac{{\partial }^2{p}^{*}}{\partial a\partial S_2}& \frac{{\partial }^2{p}^{*}}{\partial a\partial S_3}& \frac{{\partial }^2{p}^{*}}{\partial a\partial {\mathsf{\gamma }}_0}& \frac{{\partial }^2{p}^{*}}{\partial a^2}\end{array}\right]=\left[\begin{array}{ccccc}2B_1& \frac{2\mathsf{\mu }{B}_1}{\left(\mathsf{\mu }+1\right)}& \frac{2{\mathsf{\mu }}^2{B}_1}{{\left(\mathsf{\mu }+1\right)}^2}& -\frac{B_2+2e}{2}& -\frac{\mathsf{\mu }{B}_2+2\mathsf{\mu }e}{2}\\ \mathsf{\mu }{B}_3& \frac{6{\mathsf{\mu }}^2{B}_3}{\left(\mathsf{\mu }+1\right)}& \frac{6{\mathsf{\mu }}^3{B}_3}{{\left(\mathsf{\mu }+1\right)}^2}& -\mathsf{\mu }{B}_4& -{\mathsf{\mu }}^2{B}_4\\ \frac{6{\mathsf{\mu }}^2{B}_3}{\left(\mathsf{\mu }+1\right)}& \frac{6{\mathsf{\mu }}^3{B}_3}{{\left(\mathsf{\mu }+1\right)}^2}& \frac{6{\mathsf{\mu }}^4{B}_3}{{\left(\mathsf{\mu }+1\right)}^3}& -\frac{{\mathsf{\mu }}^2{B}_4}{\left(\mathsf{\mu }+1\right)}& -\frac{{\mathsf{\mu }}^3{B}_4}{\left(\mathsf{\mu }+1\right)}\\ \frac{6{\mathsf{\mu }}^3{B}_3}{{\left(\mathsf{\mu }+1\right)}^2}& \frac{6{\mathsf{\mu }}^4{B}_3}{{\left(\mathsf{\mu }+1\right)}^3}& \frac{6{\mathsf{\mu }}^5{B}_3}{{\left(\mathsf{\mu }+1\right)}^4}& -\frac{{\mathsf{\mu }}^3{B}_4}{{\left(\mathsf{\mu }+1\right)}^2}& -\frac{{\mathsf{\mu }}^4{B}_4}{{\left(\mathsf{\mu }+1\right)}^2}\\ -\mathsf{\mu }{B}_4& -\frac{{\mathsf{\mu }}^2{B}_4}{\left(\mu +1\right)}& -\frac{{\mathsf{\mu }}^3{B}_4}{{\left(\mathsf{\mu }+1\right)}^2}& 0& 0\\ -{\mathsf{\mu }}^2{B}_4& -\frac{{\mathsf{\mu }}^3{B}_4}{\left(\mathsf{\mu }+1\right)}& -\frac{{\mathsf{\mu }}^4{B}_4}{{\left(\mathsf{\mu }+1\right)}^2}& 0& 0\end{array}\right]$$

(31)

where B₁, B₂, B₃ and B₄ are defined as follows

$$\begin{array}{l}{B}_1=\frac{{\rho }_0{C}^2{\mathsf{\mu }}^2\left[1+\left(1-\frac{{\mathsf{\gamma }}_0}{2}\right)\mathsf{\mu }-\frac{a}{2}{\mathsf{\mu }}^2\right]}{{\left[1-\left(S_1-1\right)\mathsf{\mu }-S_2\frac{{\mathsf{\mu }}^2}{\mathsf{\mu }+1}-S_3\frac{{\mathsf{\mu }}^3}{{\left(\mathsf{\mu }+1\right)}^2}\right]}^3}\hspace{1em}{B}_2=\frac{{\rho }_0{C}^2{\mathsf{\mu }}^2}{{\left[1-\left(S_1-1\right)\mathsf{\mu }-S_2\frac{{\mathsf{\mu }}^2}{\mathsf{\mu }+1}-S_3\frac{{\mathsf{\mu }}^3}{{\left(\mathsf{\mu }+1\right)}^2}\right]}^2}\\ B_3=\frac{{\rho }_0{C}^2{\mathsf{\mu }}^2\left[1+\left(1-\frac{{\mathsf{\gamma }}_0}{2}\right)\mathsf{\mu }-\frac{a}{2}{\mathsf{\mu }}^2\right]}{{\left[1-\left(S_1-1\right)\mathsf{\mu }-S_2\frac{{\mathsf{\mu }}^2}{\mathsf{\mu }+1}-S_3\frac{{\mathsf{\mu }}^3}{{\left(\mathsf{\mu }+1\right)}^2}\right]}^4}\hspace{1em}{B}_4=\frac{{\rho }_0{C}^2{\mathsf{\mu }}^2}{{\left[1-\left(S_1-1\right)\mathsf{\mu }-S_2\frac{{\mathsf{\mu }}^2}{\mathsf{\mu }+1}-S_3\frac{{\mathsf{\mu }}^3}{{\left(\mathsf{\mu }+1\right)}^2}\right]}^3}\end{array}$$

(32)

If μ < 0, the derivatives in Equations (25) and (26) can be calculated by

$$\left[\begin{array}{ccccc}\frac{\partial p^{*}}{\partial S_1}& \frac{\partial p^{*}}{\partial S_2}& \frac{\partial p^{*}}{\partial S_3}& \frac{\partial p^{*}}{\partial {\mathsf{\gamma }}_0}& \frac{\partial p^{*}}{\partial a}\\ \frac{{\partial }^2{p}^{*}}{\partial S_1^2}& \frac{{\partial }^2{p}^{*}}{\partial S_1\partial S_2}& \frac{{\partial }^2{p}^{*}}{\partial S_1\partial S_3}& \frac{{\partial }^2{p}^{*}}{\partial S_1\partial {\mathsf{\gamma }}_0}& \frac{{\partial }^2{p}^{*}}{\partial S_1\partial a}\\ \frac{{\partial }^2{p}^{*}}{\partial S_2\partial S_1}& \frac{{\partial }^2{p}^{*}}{\partial S_2^2}& \frac{{\partial }^2{p}^{*}}{\partial S_2\partial S_3}& \frac{{\partial }^2{p}^{*}}{\partial S_2\partial {\mathsf{\gamma }}_0}& \frac{{\partial }^2{p}^{*}}{\partial S_2\partial a}\\ \frac{{\partial }^2{p}^{*}}{\partial S_3\partial S_1}& \frac{{\partial }^2{p}^{*}}{\partial S_3\partial S_2}& \frac{{\partial }^2{p}^{*}}{\partial S_3^2}& \frac{{\partial }^2{p}^{*}}{\partial S_3\partial {\mathsf{\gamma }}_0}& \frac{{\partial }^2{p}^{*}}{\partial S_3\partial a}\\ \frac{{\partial }^2{p}^{*}}{\partial {\mathsf{\gamma }}_0\partial S_1}& \frac{{\partial }^2{p}^{*}}{\partial {\mathsf{\gamma }}_0\partial S_2}& \frac{{\partial }^2{p}^{*}}{\partial {\mathsf{\gamma }}_0\partial S_3}& \frac{{\partial }^2{p}^{*}}{\partial {\mathsf{\gamma }}_0^2}& \frac{{\partial }^2{p}^{*}}{\partial {\mathsf{\gamma }}_0\partial a}\\ \frac{{\partial }^2{p}^{*}}{\partial a\partial S_1}& \frac{{\partial }^2{p}^{*}}{\partial a\partial S_2}& \frac{{\partial }^2{p}^{*}}{\partial a\partial S_3}& \frac{{\partial }^2{p}^{*}}{\partial a\partial {\mathsf{\gamma }}_0}& \frac{{\partial }^2{p}^{*}}{\partial a^2}\end{array}\right]=\left[\begin{array}{ccccc}0& 0& 0& e& \mathsf{\mu }e\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right]$$

(33)

The derivatives in Equations (27) and (28) can be calculated by

$$\begin{array}{l}\left[\begin{array}{ccccc}\frac{\partial {\mathsf{\sigma }}_y}{\partial A}& \frac{\partial {\mathsf{\sigma }}_y}{\partial B}& \frac{\partial {\mathsf{\sigma }}_y}{\partial C}& \frac{\partial {\mathsf{\sigma }}_y}{\partial n}& \frac{\partial {\mathsf{\sigma }}_y}{\partial m}\\ \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial A^2}& \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial A\partial B}& \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial A\partial C}& \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial A\partial n}& \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial A\partial m}\\ \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial B\partial A}& \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial B^2}& \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial B\partial C}& \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial B\partial n}& \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial B\partial m}\\ \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial C\partial A}& \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial C\partial B}& \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial C^2}& \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial C\partial n}& \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial C\partial m}\\ \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial n\partial A}& \frac{{\partial }^2{\mathrm{s}}_y}{\partial n\partial B}& \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial n\partial C}& \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial n^2}& \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial n\partial m}\\ \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial m\partial A}& \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial m\partial B}& \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial m\partial C}& \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial m\partial n}& \frac{{\partial }^2{\mathsf{\sigma }}_y}{\partial m^2}\end{array}\right]\\ =\left[\begin{array}{ccccc}{C}_2{C}_3& {\mathsf{\epsilon }}_p^n{C}_2{C}_3& C_1{C}_3\ln {\dot{\mathsf{\epsilon }}}_p^{*}& B{\mathsf{\epsilon }}_p^n{C}_2{C}_3\ln {\mathsf{\epsilon }}_p& -T^{*}{}^m{C}_1{C}_2\ln T^{*}\\ 0& 0& C_3\ln {\dot{\mathsf{\epsilon }}}_p^{*}& 0& -T^{*}{}^m{C}_2\ln T^{*}\\ 0& 0& {\mathsf{\epsilon }}_p^n{C}_3\ln {\dot{\mathsf{\epsilon }}}_p^{*}& {\mathsf{\epsilon }}_p^n{C}_2{C}_3\ln {\mathsf{\epsilon }}_p& -{\mathsf{\epsilon }}_p^n{T}^{*m}{C}_2\ln T^{*}\\ C_3\ln {\dot{\mathsf{\epsilon }}}_p^{*}& {\mathsf{\epsilon }}_p^n{C}_3\ln {\dot{\mathsf{\epsilon }}}_p^{*}& 0& B{\mathsf{\epsilon }}_p^n{C}_3{C}_5& -T^{*}{}^m{C}_1{C}_4\\ 0& {\mathsf{\epsilon }}_p^n{C}_2{C}_3\ln {\mathsf{\epsilon }}_p& B{\mathsf{\epsilon }}_p^n{C}_3{C}_5& 0& -B{\mathsf{\epsilon }}_p^n{T}^{*}{}^m{C}_2{C}_6\\ -T^{*}{}^m{C}_2\ln T^{*}& -{\mathsf{\epsilon }}_p^n{T}^{*m}{C}_2\ln T^{*}& -T^{*}{}^m{C}_1{C}_4& -B{\mathsf{\epsilon }}_p^n{T}^{*}{}^m{C}_2{C}_6& 0\end{array}\right]\end{array}$$

(34)

where C₁, C₂, C₃, C₄, C₅ and C₆ are defined as follows

$$C_1=A+B{\mathsf{\epsilon }}_P^n,C_2=1+C\ln {\dot{\mathsf{\epsilon }}}_P^{*},C_3=1-T^{*}{}^m,C_4=\ln \left(T^{*}+{\dot{\mathsf{\epsilon }}}_P^{*}\right),C_5=\ln \left({\mathsf{\epsilon }}_P^{}+{\dot{\mathsf{\epsilon }}}_P^{*}\right),C_6=\ln \left(T^{*}+{\mathsf{\epsilon }}_P^{}\right)$$

(35)

5. Results and Discussion

Two numerical examples based on steel with different stochastic parameters under explosive force are presented, in which the determined values of the parameters for steel are employed according to ref. [54]. These examples involve not only large deformation but also multi-field coupling, in which the uncertain parameters are all subject to normal distribution. In the first example, the uncertain parameters are the Mie-Grüneisen equation of state with respect to parameters $S_1$, $S_2$, $S_3$, γ₀ and $a$. In the second example, the uncertain parameters are the Johnson–Cook parameters $A$, $B$, $C$, $m$ and $n$. The probability distribution and correlation structure of the random properties should be defined through experimental measurements. However, in most cases, it is difficult to obtain the random material properties via experimental measurements. Hence, due to the lack of relevant experimental data, assumptions are made regarding these probabilistic characteristics. The Gaussian assumption is often used due to its simplicity and the lack of relevant experimental data in lots of research [35,55,56,57,58] and thus, it is adopted in this study because the focus of this paper is stochastic analysis in non-linear dynamics problem of metals with the SMPM. In practical applications, dependent random variables are often generated by some transformation of independent random variables, so it is assumed without much loss of generality that the components of the random variables, which mentioned in this study, are independent [59]. Hence, the correlation coefficient between two random variables ${\rho }_{kj}$ is 0.

Both the stochastic material point method and Monte Carlo simulation are used to calculate the mean and variance of structural response. Monte Carlo simulation is used as a validation tool. Notably, Monte Carlo simulation requires generating N random samples to guarantee result precision. Estimating N is based on the relationship between the confidence interval and sample size [60,61]. Hence, with 95% confidence degree, the sample size is 10000, which can meet the precision requirement of the two numerical examples. All calculations are processed on a personal computer with a 3.30 GHz processer and 8 GB of RAM. The calculated results are described as follows.

5.1. Example 1: The Uncertain Parameters are the State Equation Parameters

The computational domain is full of water with 2 m × 2 m × 2 m dimensions. A thick steel plate of 2 m × 2 m × 0.06 m dimensions is mounted within the domain, with the clamped edge condition assumed. TNT is a cuboid of 0.08 m × 0.08 m × 2 m dimensions and mounted on one side of the plate. Because of symmetry, the analogy model can be further simplified as a 2D planar model, as shown in Figure 3. To better monitor the stochastic response, two points, E and F, of the plate are chosen.

Table 1. The locations of two viewpoints for example1.

Location	E	F
Coordinates	(4.0 cm, 0.0 cm)	(6.0 cm, 0.0 cm)

shows the locations of points E and F of the plate. The situations of the coefficient of variation (C.V) 0.01 and 0.05 are under consideration and the values of the parameters are given in

Table 2. The values of Mie-Grüneisen equation of state parameters for steel in example 1.

Parameters	Determined Value	Mean	Standard Deviation
			C.V = 0.01	C.V = 0.15
S1	1.49	1.49	0.0149	0.2235
S2	0.0	0.0	0.0	0.0
S3	0.0	0.0	0.0	0.0
γ0	2.17	2.17	0.0217	0.3255
a	0.46	0.46	4.6 × 10−3	0.069

.

Particles distribution at the end time is shown in Figure 4. The result shows the characteristics of large deformation and multi-field coupling. Using the proposed stochastic material point method, the mean and variance of the pressure $p^{*}$ at E and F with C.V = 0.01 and C.V = 0.05 are calculated. The results of Monte Carlo simulation using 10000 samples are calculated in order to validate the accuracy of the proposed stochastic material point method. The comparative results of the two methods are shown in Figure 5 and Figure 6. A good agreement is obtained between the results of the stochastic material point method and Monte Carlo simulation. The results show the new method can well solve the stochastic analysis of metal structure involving large deformation and multi-field coupling.

Table 3. Computing time of the methods (example 1, C.V = 0.01).

Method	Grid Numbers	Particle Numbers	Computing Time
Monte Carlo	160,000	160,000	1,297,762 s
SMPM	160,000	160,000	138 s

counts the computing time of the stochastic material point method and Monte Carlo simulation with C.V = 0.01 and the comparative result of the two methods shows that the computing time for the stochastic material point method is 138 s and the Monte Carlo simulation is 1297762 s. Compared with Monte Carlo simulation, computing time of the SMPM save more than 99%; the computation efficiency of the new method is greatly increased.

5.2. Example 2: The Uncertain Parameters are the Constitutive Equation Parameters

The computational domain is full of water with 2 m × 2 m × 2 m dimensions. A thin steel plate of 2 m × 2 m × 0.02 m dimensions is mounted within the domain, with the clamped edge condition assumed. TNT is a cuboid of 0.08 m × 0.08 m × 2 m dimensions and the geometric center of the cube is the origin of coordinates. The cuboid of TNT is mounted on one side of the plate. Because of symmetry, only one section of the computing model is considered for the analysis, as shown in Figure 7. Four points-H, I, J and K are chosen to monitor the plate’s stochastic response.

Table 4. The locations of two viewpoints for example2.

Location	H	I	G	K
Coordinates	(4.5 cm, 0.0 cm)	(4.5 cm, 3.0 cm)	(5.0 cm, 0.0 cm)	(5.5 cm, 0.0 cm)

shows the locations of four points (H, I, J and K) of the plate. The situations of the coefficient of variation (C.V) 0.1 and 0.15 are under consideration. $A$, $B$, $C$, $m$ and $n$ are subject to normal distribution; mean and standard deviation of the parameters are given in

Table 5. The values of Johnson-Cook plasticity model parameters for steel in example 2.

Parameters	Determined Value	Mean	Standard Deviation
			C.V = 0.1	C.V = 0.15
A (MPa)	792	792	79.2	118.8
B (MPa)	510	510	51.0	76.5
$C$	0.014	0.014	0.0014	2.1 × 10−3
$n$	0.26	0.26	0.026	0.039
$m$	1.03	1.03	0.103	0.1545

, respectively.

Particles distribution at the end time is shown in Figure 8. The result shows the characteristics of large deformation and multi-field coupling. Using the proposed stochastic material point method, the mean and variance of the equivalent plastic stress σ_y at locations H and I with C.V = 0.1 and C.V = 0.15 are calculated, as shown in Figure 9 and Figure 10. At the same time, the deviatoric stress $S_{11}$, $S_{22}$ at locations G with C.V = 0.1 and C.V = 0.15 are calculated and shown in Figure 11 and Figure 12. The results of Monte Carlo simulation using 10000 samples are also given in these figures. A good agreement is obtained between the results of the stochastic material point method and Monte Carlo simulation.

In order to better prove the accuracy of the proposed stochastic material point method, the mean and variance of the response are also counted. The mean and variance of the equivalent plastic stress σ_y and the deviatoric stress $S_{11}$, $S_{22}$ are shown in

Table 6. The mean and variance of the equivalent plastic stress σ_y.

Location	C.V	Mean			Variance
		Monte Carlo (×10−2)	SMPM (×10−2)	Relative Errors (%)	Monte Carlo (×10−6)	SMPM (×10−6)	Relative Errors (%)
H	0.1	0.91345878	0.91237754	0.12	0.61061724	0.61631472	0.93
	0.15	0.91512557	0.91353660	0.17	1.37753594	1.38670813	0.67
I	0.1	0.94050252	0.93940429	0.12	0.62727196	0.62988312	0.42
	0.15	0.94204592	0.94042639	0.17	1.41455632	1.41723704	0.19
G	0.1	0.95229046	0.95117980	0.12	0.63159823	0.63391311	0.37
	0.15	0.95376624	0.95212650	0.17	1.42320844	1.42630451	0.22
K	0.1	0.92692249	0.92582361	0.12	0.61809317	0.62174278	0.59
	0.15	0.92852135	0.92690313	0.17	1.39413861	1.39892126	0.34

,

Table 7. The mean and variance of the deviatoric stress $S_{11}$.

Location	C.V	Mean			Variance
		Monte Carlo (×10−2)	SMPM (×10−2)	Relative Errors (%)	Monte Carlo (×10−6)	SMPM (×10−6)	Relative Errors (%)
H	0.1	−0.6301295	−0.6293716	0.12	0.29225172	0.29327005	0.35
	0.15	−0.6312935	−0.6301712	0.18	0.65820611	0.65985761	0.25
I	0.1	−0.3286691	−0.3283024	0.11	0.07606357	0.07693130	1.14
	0.15	−0.3291897	−0.3286596	0.16	0.17150211	0.17309543	0.93
G	0.1	−0.6628665	−0.6620950	0.12	0.30607402	0.30714616	0.35
	0.15	−0.6638928	−0.6627540	0.17	0.68912090	0.69107885	0.28
K	0.1	−0.6529968	−0.6522325	0.12	0.30601019	0.30857383	0.84
	0.15	−0.6541126	−0.6529930	0.17	0.69014858	0.69429113	0.60

and

Table 8. The mean and variance of the deviatoric stress S₂₂.

Location	C.V	Mean			Variance
		Monte Carlo (×10−2)	SMPM (×10−2)	Relative Errors (%)	Monte Carlo (×10−6)	SMPM (×10−6)	Relative Errors (%)
H	0.1	0.39198667	0.39154413	0.11	0.11069985	0.11350495	2.53
	0.15	0.39267753	0.39204154	0.16	0.24985084	0.25538614	2.22
I	0.1	0.12396910	0.12390616	0.05	0.01075824	0.01095824	1.86
	0.15	0.12408433	0.12404097	0.03	0.02422766	0.02465603	1.77
G	0.1	0.40287346	0.40240417	0.12	0.11283936	0.11345652	0.55
	0.15	0.40299641	0.40367725	0.17	0.25450105	0.25527716	0.30
K	0.1	0.37985846	0.37939486	0.12	0.10439505	0.10440881	0.38
	0.15	0.38052868	0.37983724	0.18	0.23554921	0.23491981	0.27

, respectively. The maximum relative error of the mean is less than 0.3% and the maximum relative error of the variance is about 3%. The results prove the correctness and accuracy of the proposed stochastic material point method.

.

In addition, the accuracy of the mean is higher than the variance. The main reason for the accuracy difference is that different number terms of the Taylor series expansion have been taken to obtain the approximate computation. In this paper, the first-three terms of the Taylor series expansion are used to calculate the approximate value of the mean but the first-two terms of the Taylor series expansion are used to calculated the approximate value of the variance. Hence, the accuracy of the mean is higher than the variance.

Table 9. Computing time of the methods (example 2, C.V = 0.1).

Method	Grid Numbers	Particle Numbers	Computing Time
Monte Carlo	84,864	160,000	874,579 s
SMPM	84,864	160,000	93 s

records the computing time of the stochastic material point method and Monte Carlo simulation with C.V = 0.1. Monte Carlo simulation takes 874,579 s and the proposed method takes 93 s. The speed is hundreds of times as fast as that of Monte Carlo simulation and the computation efficiency of the new method is greatly increased. Hence, the SMPM makes it more convenient and efficient to solve stochastic dynamics problems of metals involving large-scale structure and large deformation.

6. Conclusions

A new stochastic material point method is presented for stochastic analysis of nonlinear dynamics problems of metals. This new method can solve stochastic problems of various random factors, such as parameters of constitutive relationship and equation of state. In this aspect, the SMPM is more advanced than other stochastic meshless method which are only applicable for linear problems involving nothing but elasticity modulus of material properties. Therefore, this new method is more realistic for demonstrating the stochastic response of metal structure in stochastic dynamic analysis. To date, MPM has not received enough attention with respect to probabilistic models. Therefore, the stochastic material point method enlarges the application field of MPM.

This paper presents numerical examples to examine the accuracy and convergence of the SMPM. Good agreement is observed between the results of the SMPM and Monte Carlo simulation. Additionally, the computing time of the SMPM is much less than that of Monte Carlo simulation and the computation efficiency of the SMPM is greatly increased. Furthermore, the proposed stochastic material point method provides a new method for solving stochastic dynamics problems of metals involving large deformation and strong nonlinearity, such as hyper-velocity impact and explosion.

It is the first attempt to apply MPM in the field of stochastic analysis and thus some limitations exist, such as without considering the parameters varying across the space and inapplicability to complex models due to the limitation of two dimensions. These valuable topics are currently under development and will be subject of future works.