物质点法基础（二）：标准 MPM 的控制方程离散与计算流程

上一篇介绍了连续体的运动描述、变形梯度、变形率、柯西应力与焦曼应力率，但连续控制方程还不能直接交给计算机求解。物质点法的关键，是用物质点保存材料状态，用背景网格离散并求解动量方程，再把节点解传回物质点。这个过程同时包含空间离散、数值积分和粒子—网格双向映射，任一环节的符号或更新顺序不一致，都会造成非物理加速度、应力振荡或质量丢失。

本文只讨论最基础的单相、单套物质点标准物质点法（standard MPM）：物质点视为无尺寸积分点，背景网格采用有限元形函数，时间方向采用显式积分。为把离散主线讲清楚，暂不涉及多相与多套物质点、接触、GIMP、CPDI、APIC、隐式积分和复杂弹塑性算法，也不安排完整算例验证。

全文采用拉应力为正的连续介质符号约定。若程序采用岩土工程中常见的压应力为正约定，内力、本构和边界载荷的符号必须成套转换，不能只修改应力数组的正负号。

一、更新拉格朗日描述与控制方程

设单相连续体在时刻$t^n$占据现时构形$\Omega^n$，边界由给定位移（或速度）的$\Gamma_u^n$和给定面力的$\Gamma_t^n$组成。标准 MPM 在每个时间步内以$\Omega^n$为参考进行计算，空间导数均对当前坐标$\boldsymbol{x}^n$求取；物质点更新到$t^{n+1}$后，再以新的现时构形开始下一步。这就是更新拉格朗日格式。

1. 质量守恒

局部质量守恒方程为

$$\frac{D\rho}{Dt}+\rho\nabla\cdot\boldsymbol{v}=0, \tag{1}\label{eq:mass-local}$$

其中，$\rho$为当前密度，$\boldsymbol{v}$为速度。对随体区域$\Omega(t)$积分，可得

$$\frac{D}{Dt}\int_{\Omega(t)}\rho\,\mathrm{d}V=0. \tag{2}\label{eq:mass-integral}$$

MPM 用随材料运动的物质点表示连续体，因此每个物质点携带的质量$m_p$在计算中保持不变：

$$m_p^{n+1}=m_p^n=m_p. \tag{3}\label{eq:particle-mass}$$

质量方程通常不再像动量方程一样单独组装求解，而是通过物质点质量恒定以及体积、密度的更新自动满足：

$$\rho_p^{n+1}=\frac{m_p}{V_p^{n+1}}. \tag{4}\label{eq:density-update}$$

2. 动量守恒

忽略质量源和动量源时，现时构形中的局部线动量方程为

$$\rho\boldsymbol{a}=\nabla\cdot\boldsymbol{\sigma}+\rho\boldsymbol{b} \quad \text{in }\Omega^n, \tag{5}\label{eq:momentum-strong}$$

其中，$\boldsymbol{a}=D\boldsymbol{v}/Dt$为物质加速度，$\boldsymbol{\sigma}$为柯西应力，$\boldsymbol{b}$为单位质量体力。边界条件写成

$$\boldsymbol{v}=\bar{\boldsymbol{v}} \quad \text{on }\Gamma_u^n, \qquad \boldsymbol{\sigma}\cdot\boldsymbol{n}=\bar{\boldsymbol{t}} \quad \text{on }\Gamma_t^n, \tag{6}\label{eq:boundary-conditions}$$

其中，$\boldsymbol{n}$是边界外法向，$\bar{\boldsymbol{t}}$是给定面力。

二、动量方程的弱形式

局部方程包含应力的空间导数，直接离散要求应力场具有较高连续性。有限元和 MPM 都先将其转化为弱形式。取在$\Gamma_u^n$上为零的任意虚速度$\boldsymbol{w}$，用它点乘式 \eqref{eq:momentum-strong} 并在当前区域积分：

$$\int_{\Omega^n}\rho\boldsymbol{w}\cdot\boldsymbol{a}\,\mathrm{d}V = \int_{\Omega^n}\boldsymbol{w}\cdot(\nabla\cdot\boldsymbol{\sigma})\,\mathrm{d}V + \int_{\Omega^n}\rho\boldsymbol{w}\cdot\boldsymbol{b}\,\mathrm{d}V. \tag{7}\label{eq:weighted-momentum}$$

对第一项使用乘积求导公式

$$\nabla\cdot(\boldsymbol{\sigma}\cdot\boldsymbol{w}) = \boldsymbol{w}\cdot(\nabla\cdot\boldsymbol{\sigma}) + \nabla\boldsymbol{w}:\boldsymbol{\sigma}, \tag{8}\label{eq:product-rule}$$

再使用散度定理，可得

$$\int_{\Omega^n}\boldsymbol{w}\cdot(\nabla\cdot\boldsymbol{\sigma})\,\mathrm{d}V = \int_{\Gamma^n}\boldsymbol{w}\cdot(\boldsymbol{\sigma}\cdot\boldsymbol{n})\,\mathrm{d}\Gamma - \int_{\Omega^n}\nabla\boldsymbol{w}:\boldsymbol{\sigma}\,\mathrm{d}V. \tag{9}\label{eq:integration-by-parts}$$

由于$\boldsymbol{w}=\boldsymbol{0}$位于$\Gamma_u^n$，而$\boldsymbol{\sigma}\cdot\boldsymbol{n}=\bar{\boldsymbol{t}}$位于$\Gamma_t^n$，最终得到动量方程弱式：

$$\int_{\Omega^n}\rho\boldsymbol{w}\cdot\boldsymbol{a}\,\mathrm{d}V = - \int_{\Omega^n}\nabla\boldsymbol{w}:\boldsymbol{\sigma}\,\mathrm{d}V + \int_{\Omega^n}\rho\boldsymbol{w}\cdot\boldsymbol{b}\,\mathrm{d}V + \int_{\Gamma_t^n}\boldsymbol{w}\cdot\bar{\boldsymbol{t}}\,\mathrm{d}\Gamma. \tag{10}\label{eq:momentum-weak}$$

式 \eqref{eq:momentum-weak} 的三项分别对应惯性、内力和外力。分部积分把应力的一阶导数转移到了虚速度上，因此只需形函数的一阶导数。

三、背景网格上的有限元插值

设当前时刻共有$N_n$个相关背景网格节点。用节点形函数$N_i(\boldsymbol{x})$近似虚速度和加速度：

$$\boldsymbol{w}_h(\boldsymbol{x}) = \sum_{i=1}^{N_n}N_i(\boldsymbol{x})\boldsymbol{w}_i, \qquad \boldsymbol{a}_h(\boldsymbol{x}) = \sum_{j=1}^{N_n}N_j(\boldsymbol{x})\boldsymbol{a}_j. \tag{11}\label{eq:grid-interpolation}$$

虚速度梯度为

$$\nabla\boldsymbol{w}_h = \sum_{i=1}^{N_n}\boldsymbol{w}_i\otimes\nabla N_i. \tag{12}\label{eq:test-gradient}$$

将式 \eqref{eq:grid-interpolation} 和式 \eqref{eq:test-gradient} 代入弱式，并利用每个自由节点上的$\boldsymbol{w}_i$都是任意量，可得半离散节点方程

$$\sum_{j=1}^{N_n}M_{ij}\boldsymbol{a}_j = \boldsymbol{f}_i^{\mathrm{int}} + \boldsymbol{f}_i^{\mathrm{body}} + \boldsymbol{f}_i^{\mathrm{trac}}, \tag{13}\label{eq:semidiscrete-momentum}$$

其中，一致质量矩阵为

$$M_{ij} = \int_{\Omega^n}\rho N_iN_j\,\mathrm{d}V, \tag{14}\label{eq:consistent-mass}$$

节点内力、体力和边界面力分别为

$$\boldsymbol{f}_i^{\mathrm{int}} = - \int_{\Omega^n}\boldsymbol{\sigma}\cdot\nabla N_i\,\mathrm{d}V, \tag{15}\label{eq:internal-force-integral}$$

$$\boldsymbol{f}_i^{\mathrm{body}} = \int_{\Omega^n}\rho N_i\boldsymbol{b}\,\mathrm{d}V, \tag{16}\label{eq:body-force-integral}$$

$$\boldsymbol{f}_i^{\mathrm{trac}} = \int_{\Gamma_t^n}N_i\bar{\boldsymbol{t}}\,\mathrm{d}\Gamma. \tag{17}\label{eq:traction-force-integral}$$

以分量形式检查式 \eqref{eq:internal-force-integral}，有

$$f_{i\alpha}^{\mathrm{int}} = - \int_{\Omega^n} \sigma_{\alpha\beta} \frac{\partial N_i}{\partial x_\beta} \,\mathrm{d}V. \tag{18}\label{eq:internal-force-components}$$

式 \eqref{eq:internal-force-components} 对程序实现尤其重要：$i$是节点编号，$\alpha$是力的方向，$\beta$是求和的空间方向。内力数组的维度是“节点数$\times$空间维数”，而不是应力分量数。

四、物质点积分

MPM 与普通有限元的主要区别不在弱式，而在积分方式。连续体被离散为$N_p$个随材料运动的物质点，每个物质点$p$携带

$$\left\{ \boldsymbol{x}_p, m_p, V_p, \boldsymbol{v}_p, \boldsymbol{\sigma}_p, \boldsymbol{F}_p, \text{其他状态变量} \right\}. \tag{19}\label{eq:particle-state}$$

密度可用移动的 Dirac 函数近似表示：

$$\rho(\boldsymbol{x},t) \approx \sum_{p=1}^{N_p}m_p\delta(\boldsymbol{x}-\boldsymbol{x}_p). \tag{20}\label{eq:particle-density}$$

因此，含密度的体积分可转化为粒子求和：

$$\int_{\Omega^n}\rho f(\boldsymbol{x})\,\mathrm{d}V \approx \sum_{p=1}^{N_p}m_pf(\boldsymbol{x}_p), \tag{21}\label{eq:mass-quadrature}$$

普通体积分则使用当前物质点体积作为积分权重：

$$\int_{\Omega^n}g(\boldsymbol{x})\,\mathrm{d}V \approx \sum_{p=1}^{N_p}V_p^ng(\boldsymbol{x}_p^n). \tag{22}\label{eq:volume-quadrature}$$

定义当前物质点位置处的形函数及其空间梯度

$$N_{ip}^n=N_i(\boldsymbol{x}_p^n), \qquad \boldsymbol{G}_{ip}^n = \nabla N_i(\boldsymbol{x}_p^n). \tag{23}\label{eq:shape-at-particle}$$

标准 MPM 将物质点视为无尺寸积分点，$N_{ip}^n$和$\boldsymbol{G}_{ip}^n$只由点位置及其所在背景单元决定。例如，一维线性单元$[x_a,x_b]$的长度为$h=x_b-x_a$，则

$$N_a(x)=\frac{x_b-x}{h}, \qquad N_b(x)=\frac{x-x_a}{h}, \qquad \frac{\mathrm{d}N_a}{\mathrm{d}x}=-\frac{1}{h}, \qquad \frac{\mathrm{d}N_b}{\mathrm{d}x}=\frac{1}{h}. \tag{24}\label{eq:linear-shape}$$

当物质点越过单元边界时，其相关节点和形函数梯度会切换。这种点式特征正是标准 MPM 简单高效的原因，也是网格穿越误差的来源。

五、从物质点映射到背景网格

1. 节点质量

把粒子积分式 \eqref{eq:mass-quadrature} 代入一致质量矩阵式 \eqref{eq:consistent-mass}，得到

$$M_{ij}^n = \sum_{p=1}^{N_p}m_pN_{ip}^nN_{jp}^n. \tag{25}\label{eq:particle-consistent-mass}$$

显式 MPM 通常不求解完整质量矩阵，而是利用形函数单位分解性质$\sum_jN_{jp}=1$进行行和集中：

$$\begin{aligned} m_i^n &=\sum_jM_{ij}^n\\ &=\sum_p m_pN_{ip}^n\sum_jN_{jp}^n\\ &=\sum_p m_pN_{ip}^n. \end{aligned} \tag{26}\label{eq:lumped-mass}$$

节点质量只在当前时间步有效。背景网格不永久保存材料，下一步开始时必须清零并重新映射。

2. 节点动量与速度

物质点动量按相同形函数映射到节点：

$$\boldsymbol{p}_i^n = \sum_p m_pN_{ip}^n\boldsymbol{v}_p^n. \tag{27}\label{eq:nodal-momentum}$$

对有效质量节点，节点速度为

$$\boldsymbol{v}_i^n = \frac{\boldsymbol{p}_i^n}{m_i^n}. \tag{28}\label{eq:nodal-velocity}$$

程序中应设置相对于总质量或材料密度尺度的质量容差。若$m_i$小于容差，该节点不参与除法和后续更新，避免产生数值极大的速度与加速度。

3. 节点内力与外力

由式 \eqref{eq:volume-quadrature}，节点内力离散为

$$\boldsymbol{f}_i^{\mathrm{int},n} = - \sum_pV_p^n \boldsymbol{\sigma}_p^n\cdot\boldsymbol{G}_{ip}^n. \tag{29}\label{eq:particle-internal-force}$$

单位质量体力离散为

$$\boldsymbol{f}_i^{\mathrm{body},n} = \sum_pm_pN_{ip}^n\boldsymbol{b}_p^n. \tag{30}\label{eq:particle-body-force}$$

若重力是唯一体力，则$\boldsymbol{b}_p=\boldsymbol{g}$。给定面力仍来源于边界积分式 \eqref{eq:traction-force-integral}。在实际程序中，可对边界段或边界面作数值积分，也可把已知合力按边界形函数分配到相关节点；不能把面力直接当作体力施加到所有物质点，否则载荷会随离散厚度或物质点数量变化。

节点总力为

$$\boldsymbol{f}_i^n = \boldsymbol{f}_i^{\mathrm{int},n} + \boldsymbol{f}_i^{\mathrm{body},n} + \boldsymbol{f}_i^{\mathrm{trac},n}. \tag{31}\label{eq:total-force}$$

六、背景网格上的显式动量更新

采用集中质量后，式 \eqref{eq:semidiscrete-momentum} 在每个节点上解耦：

$$\boldsymbol{a}_i^n = \frac{\boldsymbol{f}_i^n}{m_i^n}. \tag{32}\label{eq:nodal-acceleration}$$

以前向欧拉形式更新节点动量：

$$\boldsymbol{p}_i^{n+1} = \boldsymbol{p}_i^n +\Delta t\,\boldsymbol{f}_i^n, \tag{33}\label{eq:nodal-momentum-update}$$

相应的更新后节点速度为

$$\boldsymbol{v}_i^{n+1} = \frac{\boldsymbol{p}_i^{n+1}}{m_i^n} = \boldsymbol{v}_i^n+\Delta t\,\boldsymbol{a}_i^n. \tag{34}\label{eq:nodal-velocity-update}$$

速度边界条件在节点上施加。以固定边界的法向方向为例，应将该方向上的$\boldsymbol{v}_i^n$、$\boldsymbol{a}_i^n$和$\boldsymbol{v}_i^{n+1}$一致地约束为零。若只在更新前把速度清零，节点仍可能在内力作用下重新获得法向速度；若只在更新后清零，则用于物质点速度增量的节点加速度仍可能违反边界条件。

显式积分的时间步受稳定条件限制。线弹性问题可用下式估计：

$$\Delta t \leq \alpha\min_e\left(\frac{h_e}{c_e}\right), \qquad c_e=\sqrt{\frac{K_e+4G_e/3}{\rho_e}}, \tag{35}\label{eq:cfl-condition}$$

其中，$h_e$为单元特征长度，$K_e$和$G_e$分别为体积模量与剪切模量，$\alpha<1$为安全系数。大变形、强非线性或高速运动问题还应采用更保守的时间步估计。

七、从背景网格映射回物质点

为给出完整而单一的更新顺序，本文采用显式应力最后更新（update stress last，USL）流程，不比较其他时间更新变体。

1. 物质点速度

物质点速度使用节点加速度的插值增量更新：

$$\boldsymbol{v}_p^{n+1} = \boldsymbol{v}_p^n + \Delta t\sum_iN_{ip}^n\boldsymbol{a}_i^n. \tag{36}\label{eq:particle-velocity-update}$$

式 \eqref{eq:particle-velocity-update} 更新的是材料点自身原有速度，而不是简单用节点速度覆盖物质点速度，有利于减少映射造成的数值耗散。

2. 物质点位置

使用更新后的节点速度推进物质点位置：

$$\boldsymbol{x}_p^{n+1} = \boldsymbol{x}_p^n + \Delta t\sum_iN_{ip}^n\boldsymbol{v}_i^{n+1}. \tag{37}\label{eq:particle-position-update}$$

物质点移动后必须重新搜索所在单元，并在下一时间步重新计算相关节点、$N_{ip}$和$\boldsymbol{G}_{ip}$。规则结构网格可以由坐标直接计算单元编号，非规则网格则需要可靠的单元搜索和越界处理。

3. 速度梯度与变形率

背景网格速度场近似为

$$\boldsymbol{v}_h(\boldsymbol{x}) = \sum_iN_i(\boldsymbol{x})\boldsymbol{v}_i^{n+1}. \tag{38}\label{eq:grid-velocity-field}$$

对空间坐标求梯度，并在物质点处取值，得到

$$\boldsymbol{L}_p^{n+1} = \left.\nabla\boldsymbol{v}_h\right|_{\boldsymbol{x}_p^n} = \sum_i\boldsymbol{v}_i^{n+1}\otimes\boldsymbol{G}_{ip}^n. \tag{39}\label{eq:particle-velocity-gradient}$$

速度梯度可分解为对称变形率张量和反对称旋率张量：

$$\boldsymbol{D}_p = \frac{1}{2}\left(\boldsymbol{L}_p+\boldsymbol{L}_p^{\mathrm{T}}\right), \qquad \boldsymbol{W}_p = \frac{1}{2}\left(\boldsymbol{L}_p-\boldsymbol{L}_p^{\mathrm{T}}\right). \tag{40}\label{eq:velocity-gradient-split}$$

4. 变形梯度、体积与密度

由$\dot{\boldsymbol{F}}=\boldsymbol{L}\boldsymbol{F}$，变形梯度可作一阶更新：

$$\boldsymbol{F}_p^{n+1} = \left(\boldsymbol{I}+\Delta t\,\boldsymbol{L}_p^{n+1}\right) \boldsymbol{F}_p^n. \tag{41}\label{eq:deformation-gradient-update}$$

当前体积由变形梯度行列式控制。若$V_p^0$是初始体积，则

$$J_p^{n+1}=\det\boldsymbol{F}_p^{n+1}, \qquad V_p^{n+1}=J_p^{n+1}V_p^0. \tag{42}\label{eq:particle-volume-from-j}$$

也可逐步写成

$$V_p^{n+1} = \det\left(\boldsymbol{I}+\Delta t\,\boldsymbol{L}_p^{n+1}\right)V_p^n. \tag{43}\label{eq:incremental-volume-update}$$

小变形条件下，式 \eqref{eq:incremental-volume-update} 的一阶近似为

$$V_p^{n+1} \approx \left(1+\Delta t\,\operatorname{tr}\boldsymbol{D}_p^{n+1}\right)V_p^n. \tag{44}\label{eq:small-strain-volume}$$

更新体积后，再由式 \eqref{eq:density-update} 计算密度。质量$m_p$、体积$V_p$和密度$\rho_p$三者中，质量必须保持不变，不能同时独立更新三个量。

5. 应变与应力

在小变形线弹性范围内，应变增量为

$$\Delta\boldsymbol{\varepsilon}_p = \Delta t\,\boldsymbol{D}_p^{n+1}, \tag{45}\label{eq:strain-increment}$$

各向同性材料的应力增量为

$$\Delta\boldsymbol{\sigma}_p = \lambda\operatorname{tr}(\Delta\boldsymbol{\varepsilon}_p)\boldsymbol{I} +2G\Delta\boldsymbol{\varepsilon}_p, \tag{46}\label{eq:linear-elastic-increment}$$

其中，$\lambda$和$G$为 Lamé 常数。于是

$$\boldsymbol{\sigma}_p^{n+1} = \boldsymbol{\sigma}_p^n+\Delta\boldsymbol{\sigma}_p. \tag{47}\label{eq:small-strain-stress-update}$$

若考虑有限转动，即使材料仍为弹性，也不能忽略应力客观性。采用上一篇介绍的焦曼应力率

$$\overset{\nabla}{\boldsymbol{\sigma}} = \dot{\boldsymbol{\sigma}} -\boldsymbol{W}\boldsymbol{\sigma} +\boldsymbol{\sigma}\boldsymbol{W}, \tag{48}\label{eq:jaumann-rate}$$

并令$\overset{\nabla}{\boldsymbol{\sigma}}=\mathbb{C}:\boldsymbol{D}$，则一阶显式更新可写为

$$\boldsymbol{\sigma}_p^{n+1} = \boldsymbol{\sigma}_p^n +\Delta t \left[ \mathbb{C}:\boldsymbol{D}_p +\boldsymbol{W}_p\boldsymbol{\sigma}_p^n -\boldsymbol{\sigma}_p^n\boldsymbol{W}_p \right]. \tag{49}\label{eq:jaumann-stress-update}$$

式 \eqref{eq:jaumann-stress-update} 只用于说明标准 MPM 中运动学量如何进入本构更新，并不代替针对具体材料设计和验证的应力积分算法。

八、标准 MPM 的单步计算流程

将前面的离散关系组合起来，一个 USL 显式时间步可整理为以下顺序：

已知 tⁿ 时刻的物质点状态：
位置、质量、体积、速度、应力、变形梯度和材料状态变量

1. 清空背景网格上的质量、动量、内力和外力
2. 搜索每个物质点所在单元，计算 Nᵢₚ 和 ∇Nᵢₚ
3. 物质点 → 网格：
   3.1 映射节点质量 mᵢ
   3.2 映射节点动量 pᵢ
   3.3 组装节点内力、体力和面力
4. 对有效质量节点计算速度、总力和加速度
5. 在节点上施加速度边界条件
6. 显式更新节点动量和节点速度
7. 网格 → 物质点：
   7.1 更新物质点速度
   7.2 更新物质点位置
   7.3 计算物质点速度梯度
   7.4 更新变形梯度、体积和密度
   7.5 调用本构关系更新应力与状态变量
8. 丢弃当前节点状态，进入下一时间步

背景网格在步骤 8 中被“丢弃”，是指清除节点上的临时材料量，而不是删除网格拓扑。下一步仍可使用原来的规则背景网格，但节点质量、动量和力必须由更新后的物质点重新生成。

九、离散关系的基本自检

本文不展开完整验证算例，但实现时仍可利用形函数性质检查最基本的守恒关系。若形函数满足单位分解

$$\sum_iN_{ip}=1, \qquad \sum_i\boldsymbol{G}_{ip}=\boldsymbol{0}, \tag{50}\label{eq:shape-properties}$$

则粒子到网格映射应满足总质量守恒：

$$\sum_im_i = \sum_i\sum_pm_pN_{ip} = \sum_pm_p. \tag{51}\label{eq:grid-mass-conservation}$$

节点总动量也应等于映射前的物质点总动量：

$$\sum_i\boldsymbol{p}_i = \sum_i\sum_pm_pN_{ip}\boldsymbol{v}_p = \sum_pm_p\boldsymbol{v}_p. \tag{52}\label{eq:grid-momentum-conservation}$$

在没有外力的完整物体内部，节点内力之和应为零：

$$\sum_i\boldsymbol{f}_i^{\mathrm{int}} = - \sum_pV_p\boldsymbol{\sigma}_p\cdot \left(\sum_i\boldsymbol{G}_{ip}\right) = \boldsymbol{0}. \tag{53}\label{eq:internal-force-balance}$$

这些等式不能替代解析解、网格收敛性和能量误差验证，但能快速发现单元搜索遗漏、形函数编号错误、梯度方向错误以及节点清零不完整等实现问题。

十、标准 MPM 的适用范围与局限

标准 MPM 保留了有限元弱式和形函数插值，同时利用移动物质点避免网格随材料发生严重畸变。对于大位移问题，背景网格能够反复使用，历史变量始终跟随物质点保存。

它的局限也与这种最简离散直接相关。物质点被视为无尺寸积分点，线性单元的形函数梯度在单元边界不连续；物质点跨越单元时，相关节点和梯度突然改变，可能引起网格穿越误差。物质点分布不均还会降低粒子积分精度，单元内物质点过少时可能出现积分不足。本文只说明这些问题的来源，不引入任何改进形函数或粒子域方法。

从控制方程到程序数组，标准 MPM 的主线可以概括为：先在现时构形建立动量方程弱式，再用背景网格形函数完成空间离散，用物质点完成积分，随后在节点上显式求解动量，最后把运动和变形信息传回物质点。只要质量、内力、边界条件和更新时间层保持一致，这套基础流程就构成后续接触、多相耦合和改进 MPM 方法的共同起点。

参考资料

1. Sulsky, D., Chen, Z., & Schreyer, H. L. A particle method for history-dependent materials. Computer Methods in Applied Mechanics and Engineering, 1994, 118(1–2): 179–196.
2. Sulsky, D., Zhou, S. J., & Schreyer, H. L. Application of a particle-in-cell method to solid mechanics. Computer Physics Communications, 1995, 87(1–2): 236–252.
3. Bardenhagen, S. G., & Kober, E. M. The generalized interpolation material point method. Computer Modeling in Engineering & Sciences, 2004, 5(6): 477–495.

本文含有 AI 生成内容