![](https://csdnimg.cn/release/download_crawler_static/34023915/bg4.jpg)
4 Z. Sun et al.
where Ω is the current configuration of the continuum
body; ρ(x), v(x) and σ(x) denote the initial mass density,
velocity and initial stress of the continuum body,
respectively; ∂Ω
u
and ∂Ω
τ
are the parts of the boundary
subject to the prescribed displacement and traction,
respectively; u(t) and τ
b
(t) represent the prescribed
displacement and traction, respectively. n
b
is the outer
norm vector of the traction boundary.
2.1 MPM formulation
In the MPM, a continuum body is represented by a finite
number of discrete Lagrangian particles (also called
material points). The deformations of continuum body are
described by tracking the movements of material points.
Since the mass of each material point remains unchanged
over time, the conservation equation of mass equation (1) is
automatically satisfied. The disretised MPM equation of
momentum is developed by the weak form of equation (2),
i.e.,
:
τ
s
b
d ρ ddS
ρ d
ΩΩ ∂Ω
Ω
⋅Ω=− ∇Ω+ ⋅
+⋅Ω
∫∫ ∫
∫
wa σ w τ w
wb
(8)
where ρσ
s
= σ, and w denotes the test function.
At time level k, each material point carries the material
properties and state variables, including mass M
p
, density
,
k
velocity ,
k
v stress ,
k
σ strain ,
k
ε etc. Thus, the
integrals in equation (8) are calculated by the sums of
quantities at material points,
[][]
[]
[][]
1
()
1
1
(), (),
(), : ( , )|
(), (),
p
p
p
τ
p
N
pp p
p
N
s
pp Xt b
p
N
pp p
p
Mtttt
tt t dS
Mtttt
=
=
∂Ω
=
⋅=
−∇+⋅
+⋅
∑
∑
∫
∑
wX aX
σ Xwx wτ
wX bX
(9)
where X
p
(t) is the position of material point p at time t.
Meanwhile, an Eulerian background mesh as shown in
Figure 1 is defined to calculate the gradient terms and solve
the momentum equation.
By using the conventional finite element shape function,
the filed variables and test function at the material points are
interpolated by the corresponding values at the background
mesh nodes. Then, with the use of lumped mass, the
discretised MPM equation of momentum becomes
int, ,kextk
kk
ii i i
m =+af f (10)
int, ,
1
p
N
ksk
k
pp
iip
p
M
=
=− ⋅
∑
fGσ (11)
,
ˆ
ext k
kk
iii
+fbτ (12)
()
kk
ip
ip
N=∇GX
(13)
where
()
1
p
N
kk
pi p
i
p
mMN
=
=
∑
denotes the lumped mass the
i
th
grid node at time level k;
int,k
i
f and
,ext k
i
f represent the
i
th
grid node internal and external force tensor at time level
k, respectively;
()
k
ip
N
is the shape function of node i
associated with spatial material points X
p
at time level k.
In the present study, an explicit Euler algorithm is used
to solve equation (10) for nodal accelerations. The grid node
velocities at time level k + 1, are obtained by
1kkk
iii
t
+
+Δvva (14)
Then, the positions and velocities of the material point p at
time level k + 1,
1t
p
and
1
,
k
p
+
v can be calculated from the
interpolations of grid node velocities and accelerations,
respectively
()
1
1
n
N
kk kk
pp ip
i
i
Nt
+
=
+Δ
∑
XX vX (15)
()
1
1
n
N
kk k k
pp ip
i
i
Nt
+
=
+Δ
∑
vv aX (16)
where Δt is the time step for the explicit integration. The
strain increments are updated by mapping the updated
particle moment back to the grid nodes
()
()
()
1
1
1
p
N
k
k
k
ii p p i p
p
mMN
+
+
=
=
∑
vvX
(17)
()
111
2
p
T
kkkkk
ip i ip i
t
+++
Δ
⎤
Δ= +
⎥
⎦
ε Gv Gv
(18)
The stress increment can be derived from the given strain
increment based on the constitutive model, and then used to
update the particle stress tensor.
)
11kk k
pp p
f
++
=+Δσσ ε
(19)
where
f describes the constitutive law of the material.