0 引　言

结构与材料的大变形问题模拟与仿真具有实际意义，在许多技术和工程领域具有重要意义。如在造船和海洋工程中，典型的大变形问题包括船体过载和船体与浮冰之间的相互作用^[1]，以及材料成型和制造工程。

为了解决大变形问题，在过去几十年中提出了各种数值方法。有网格方法例如任意拉格朗日-欧拉理论^[2]，无网格方法例如颗粒有限元法^[3]、物质点法^[4]、光滑粒子流体动力学^[5]以及近场动力学方法^[6]。

近场动力学（Peridynamics，PD）作为经典连续介质力学（Classical Continuum Mechanics，CCM）的重构，引入了独特的非局部框架，与分子动力学类似，考虑了成对粒子之间的非局部相互作用力，并建立了近场动力学的拉格朗日框架。在近场动力学中，任何一对粒子都存在一个力密度函数，并且假设粒子可在一个有限的范围内相互作用，该范围称为近场（horizon）。材料点处的应变能密度由其相邻材料点的变形确定。与CCM不同，PD框架中不需要连续场。因此，PD在处理不连续问题（如裂纹和损伤）方面具有优势^[7-17]。

近场动力学根据键力的不同分为键基近场动力学和态基近场动力学，态基近场动力学又可分为常规态基近场动力学和非常规态基近场动力学。键基近场动力学模型普遍存在泊松比的限制，态基近场动力学通过构造非局部变形梯度建立与经典连续介质力学的联系，但变形梯度会存在奇异解。为此，本文从分子动力学角度从发，构造物质点间新型的微势函数，推导全新的本构模型，对研究近场动力学的适用性具有重要意义。

1 微势近场动力学模型

为了避免现有键基近场动力学和非常规态基近场动力学存在的问题，并丰富现有常规态基近场动力学键力模型，在Fan等^[18]工作基础上，进一步探究微势近场动力学模型（Micro-Potential Based Peridynamics，MPPD）的适用性。

1.1 本构模型推导 1.1.1 微观势能函数

对于冰、岩石、陶瓷和混凝土等岩土材料，破坏过程受粘性断裂定律控制。为了模拟粘性变形，采用指数-指数微势函数：

$w\left( {{\Delta _{{n}}},{\Delta _{{t}}}} \right) = \phi - \phi \left( {1 + \frac{{{\Delta _{{n}}}}}{{{\delta _{{n}}}}}} \right)\exp \left( { - \frac{{{\Delta _{{n}}}}}{{{\delta _{{n}}}}} - \frac{{\Delta _{{t}}^2}}{{\delta _{{t}}^2}}} \right)。$

(1)

式中： ${\Delta _{{n}}}$ 和 ${\Delta _{{t}}}$ 分别为键长变化法向变形分量和切向变化分量； $\phi$ ， ${\delta _{{n}}}$ 和 ${\delta _{{t}}}$ 为材料参数。其中， $\phi$ 的量纲为MT⁻²L⁻⁴， ${\delta _{{n}}}$ 和 ${\delta _{{t}}}$ 的量纲为L⁻²。注意，键基近场动力学中微势能仅取决于键的初始长度及其变形大小。在MPPD模型中，粒子的微势能不仅与键的初始长度有关，还与初始构型中的取向有关。此外，微势能是形变状态 ${{\boldsymbol{Y}}}$ 的泛函，它由其近场内所有相互关联的形变决定。

1.1.2 变形度量

在MPPD模型中引入键长的二阶变形：

$\Delta {\text{ = }}{\left| {{{\boldsymbol{Y}}}\left\langle {\xi } \right\rangle } \right|^2} - {\left| {\xi } \right|^2}。$

(2)

式中： ${\xi }$ 为参考构型下键的长度； ${{\boldsymbol{Y}}}\langle {\xi }\rangle$ 是当前构型下键的变形状态。在经典连续介质力学理论中，变形状态与未变形的键长有关，遵循Cauchy-Born法则：

${{\boldsymbol{Y}}}\left\langle {\xi } \right\rangle = {{\boldsymbol{F}}\xi } 。$

(3)

式中， ${{\boldsymbol{F}}}$ 为变形梯度。

用式（3）代替式（2）中的变形状态，得到：

$\Delta = 2{\xi ^{\text{T}}}{\boldsymbol{E}}\xi 。$

(4)

式中， ${\boldsymbol{E}}$ 为应变张量。

在MPPD模型中，提出了一种广义的非局部应变张量 ${{\boldsymbol{E}}}(x,t)$ ，定义为：

${{\boldsymbol{E}}}(x,t) = \frac{1}{4}\int_\mathcal{H} \omega \langle {\xi }\rangle \left( {\frac{{{{\boldsymbol{Y}}}\langle {\xi }\rangle {{\boldsymbol{Y}}}\langle {\xi }\rangle }}{{{{\xi }^{\text{T}}}{\xi }}} - 1} \right)\left( {5\frac{{{\xi } \otimes {\xi }}}{{{{\xi }^{\text{T}}}{\xi }}} - 1} \right){\rm{d}}{V_{\mathbf{\xi }}}。$

(5)

式中，影响函数 $\omega \langle \xi \rangle = \dfrac{{315}}{{16\pi }}{\delta ^{ - 7}}(\delta + 3\xi ){(\delta - \xi )^3}$ ，用于表示数值计算中的某些物质点之间相互作用权重。

结合CCM理论，非局部偏应变可表述为：

${{\boldsymbol{E}}'} = {{\boldsymbol{E}}} - \frac{{{E_{kk}}}}{3}{{\boldsymbol{I}}} 。$

(6)

式中， ${E_{kk}}$ 为应变张量的迹，定义为

${E_{kk}} = \frac{1}{2}\int_H \omega \langle {\xi }\rangle \left( {\frac{{{{\boldsymbol{Y}}}\langle {\xi }\rangle {{\boldsymbol{Y}}}\langle {\xi }\rangle }}{{{{\xi }^{\text{T}}}{\xi }}} - 1} \right){\rm{d}}{V_{\mathbf{\xi }}} 。$

(7)

把式（6）和式（7）代入式（4）可得：

$\Delta = \frac{2}{3}{E_{kk}}{{\xi }^{\text{T}}}{\xi } + 2{{\xi }^{\text{T}}}{{\boldsymbol{E}}'\xi } 。$

(8)

式（8）右边第一项可看作由于体应变导致的键的变形，记为键长变化的法向分量：

${\Delta _n} = \frac{2}{3}{E_{kk}}{{\xi }^{\text{T}}}{\xi } 。$

(9)

第二项可看作由于偏应变导致的键的变形，记为键长变化的切向分量：

${\Delta _t} = 2{{\xi }^{\text{T}}}{{\boldsymbol{E}}'\xi }。$

(10)

1.1.3 本构关系

键基近场动力学中当键长发生变化时，微势能被储存在键中，从而导致单个粒子的宏观非局部弹性应变能密度（Nonlocal Elastic Strain Energy Density，NESED），定义为：

$W = \frac{1}{2}\int_\mathcal{H} w (\eta ,\xi ){\rm{d}}{V_\xi }。$

(11)

为了得到MPPD的本构模型，首先取式（11）的变分：

$\delta {W^{non}} = \int_\mathcal{H} \delta w{\rm{d}}{V_\xi } = \frac{1}{2}\int_\mathcal{H} {\left( {{F_n}\delta {\Delta _n} + {F_t}\delta {\Delta _t}} \right)} {\rm{d}}{V_\xi } 。$

(12)

式中： ${F_n}\left( {{\Delta _n},{\Delta _t}} \right) = \dfrac{{\partial w}}{{\partial {\Delta _n}}}$ ， ${F_t}\left( {{\Delta _n},{\Delta _t}} \right) = \dfrac{{\partial w}}{{\partial {\Delta _t}}}$ 分别表示法向力态分量和切向力态分量。

将式（9）和式（10）代入式（12），经过计算可得：

$\delta {W^{{\text{non}}}} = \frac{1}{2}\int_\mathcal{H} {\left[ {\frac{2}{3}\left( {{F_n} - {F_t}} \right)|\xi {|^2}\delta {E_{kk}} + 2{F_t}{\xi ^{\text{T}}}\delta {{\boldsymbol{E}}}\xi } \right]} {\rm{d}}{V_\xi }。$

(13)

可以发现：

${\xi ^{\text{T}}}\delta {{\boldsymbol{E}}}\xi = \frac{1}{2}\delta \Delta = \frac{1}{2}\delta \left( {|{Y}\langle \xi \rangle {|^2} - {\xi ^2}} \right) = {{\boldsymbol{Y}}}\langle \xi \rangle \cdot \delta {{\boldsymbol{Y}}}\langle \xi \rangle。$

(14)

其中：

$\delta {{\boldsymbol{E}}} = \frac{1}{2}\int_\mathcal{H} \omega \langle {\xi }\rangle \frac{{{{\boldsymbol{Y}}}\langle {\xi }\rangle \delta {{\boldsymbol{Y}}}\langle {\xi }\rangle }}{{{{\xi }^{\text{T}}}{\xi }}}\left( {5\frac{{{\xi } \otimes {\xi }}}{{{{\xi }^{\text{T}}}{\xi }}} - 1} \right){\rm{d}}{V_{\xi }} ，$

(15)

$\delta {E_{kk}} = \frac{1}{2}\int_H \omega \langle {\xi }\rangle \frac{{{{\boldsymbol{Y}}}\langle {\mathbf{\xi }}\rangle \delta {{\boldsymbol{Y}}}\langle {\xi }\rangle }}{{{{\xi }^{\text{T}}}{\xi }}}{\rm{d}}{V_{\xi }} 。$

(16)

把式（15）和式（16）代入式（13），可得：

$\begin{split} \delta {W^{{\text{non}}}} =& \frac{1}{3}\int_\mathcal{H} {\left( {{F_n} - {F_t}} \right)} {\xi ^2}{\rm{d}}{V_\xi }\int_\mathcal{H} \omega \langle \beta \rangle \frac{{{{\boldsymbol{Y}}}\langle \beta \rangle \cdot \delta {{\boldsymbol{Y}}}\langle \beta \rangle }}{{{\beta ^{\text{T}}}\beta }}{\rm{d}}{V_\beta }+ \\ & \int_\mathcal{H} {{F_t}} {{\boldsymbol{Y}}}\langle \xi \rangle \cdot \delta {{\boldsymbol{Y}}}\langle \xi \rangle {\rm{d}}{V_\xi }。\\[-10pt] \end{split}$

(17)

式中，使用 ${A}$ 和 ${B}$ 两个状态的点积，定义为

${A} \cdot {{\boldsymbol{B}}}: = \int_\mathcal{H} {({{\boldsymbol{AB}}})} \langle \xi \rangle {\rm{d}}{V_\xi }。$

(18)

通过非局部弹性应变能密度的导数可得到本构关系，第一项进行积分变换 $\beta \leftrightarrow \xi$ ，可得：

$\begin{split} \delta {W^{{\text{non }}}} = & \int_\mathcal{H} {\left[ {\frac{1}{3}\omega \langle \xi \rangle {\xi ^{ - 2}}\int_\mathcal{H} {\left( {{F_n} - {F_t}} \right)} {\beta ^2}d{V_\beta } + {F_t}} \right]}\;\times\\ & {{\boldsymbol{Y}}}\langle \xi \rangle \delta {{\boldsymbol{Y}}}\langle \xi \rangle {\rm{d}}{V_\xi }。\end{split}$

(19)

根据式（19），可得本构关系：

${{\boldsymbol{T}}}\langle \xi \rangle = \left[ {\frac{1}{3}\omega (\xi ){\xi ^{ - 2}}\int_\mathcal{H} {\left( {{F_n} - {F_t}} \right)} {\beta ^2}{\rm{d}}{V_\beta } + {F_t}(\xi )} \right]{{\boldsymbol{Y}}}\langle \xi \rangle 。$

(20)

式中，力态可分解为

${{\boldsymbol{T}}}\langle \xi \rangle = {{{\boldsymbol{T}}}_1}\langle \xi \rangle + {{{\boldsymbol{T}}}_2}\langle \xi \rangle。$

(21)

式中： ${{{\boldsymbol{T}}}_1}(\xi ) = \dfrac{1}{3}\omega (\xi ){\xi ^{ - 2}}\left[ {\displaystyle\int_\mathcal{H} {\left( {{F_n} - {F_t}} \right)} {\beta ^2}{\rm{d}}{V_\beta }} \right]{{\boldsymbol{Y}}}\langle \xi \rangle$ ， ${{{\boldsymbol{T}}}_2}\langle \xi \rangle = {F_t}(\xi ){{\boldsymbol{Y}}}\langle \xi \rangle$ 。T₁表示族成员对单个键力密度的共同贡献，T₂可看作是单个键变形产生的力态作用。此外，T₁和T₂与变形键平行。因此，可将该模型看作常规态基近场动力学本构模型。

1.2 键参数标定

幂指数型微势的应用很难得到描述键合力密度的具体表达式。为使问题具体化，在MPPD模型中，通过分析各向同性均匀变形的平面应力场来标定键参数。本构方程 $\phi$ 、 ${\delta _n}$ 和 ${\delta _t}$ 是需要确定的3个未知参数。采用等效应变能密度法建立参数与宏观材料常数（如杨氏模量和泊松比）之间的关系。

假设材料经历一个各向同性的弹性无穷小变形，这时非局部弹性应变能密度与在经典的连续介质力学框架下各向同性下的线弹性应变能密度相等，可得：

$\frac{\phi }{{\delta _n^2}} = \frac{{21}}{{16{\text{π}}{\delta ^7}}}(3\lambda + 2\mu ) = \frac{{63}}{{16{\text{π}} {\delta ^7}}}K 。$

(22)

式中，K为体积模量， $\lambda$ 和 $\mu$ 为拉梅系数。

同样，假设材料不发生任何体积变化。可得如下关系：

$\frac{{\delta _n^2}}{{\delta _t^2}} = \frac{{5(1 - 2\nu )}}{{2(1 + \nu )}} 。$

(23)

1.3 键破坏准则

物体 $\mathcal{B}$ 通过变形状态 ${{{\boldsymbol{Y}}}_{\text{n}}}$ 从初始构型 ${\varOmega _0}$ 变换至当前构型 ${\varOmega _n}$ 。微势能随变形而产生，并储存在连接粒子 $x$ 和 $x'$ 的变形键中，粒子 $x$ 的微势能用 $w\left( {x,\xi ,{{{\boldsymbol{Y}}}_n}} \right)$ 计算，粒子 $x'$ 的微势能用 $w\left( {x', - \xi ,{{{\boldsymbol{Y}}}_n}} \right)$ 计算。基于能量建立判断键的失效准则。

对于二维问题，断裂穿过单位裂纹面所有键所需的能量如下：

${G_{{\text{II}}0}} = \int_0^\delta {\int_0^t {\int_z^\delta {\int_0^{{{\cos }^{ - 1}}z/\xi } {{w_c}\xi \sin \theta } } } } {\text{d}}\theta {\text{d}}\xi {\text{d}}x{\text{d}}z。$

(24)

式中， ${w_c}$ 键失效的临界微势能。图1显示了该计算的图解，其中t为厚度。

计算式（24），得到能量释放率为：

${G_{{\text{II}}0}} = \frac{{2{\delta ^3}t{w_c}}}{3} 。$

(25)

由于能量释放率可通过脆性材料的试验测试获得，求解式（25）得到二维问题中键失效的临界能为：

${w_c} = \frac{{3{G_{{\text{II}}0}}}}{{2{\delta ^3}t}}。$

(26)

在当前的断裂准则中，当总微势能 $w\left( {x,x',{{\boldsymbol{Y}}_n}} \right)$ 大于临界微势 ${w_c}$ 时，连接 $x$ 和 $x'$ 的键将永远失效。为此，引入一个依赖历史的标量值函数，记为：

$\mu \left(x,{x}^{\prime },t\right)=\left\{\begin{array}{ll}1，& w\left(x,{x}^{\prime },{\text{Y}}_{n}\right) < {w}_{c}，\\ 0，\text{}\text{} & \text{ }{\rm{others}}\text{ }。\end{array} \right.$

(27)

2 数值算法

基于近场动力学公式，Silling等开发了近场动力学计算框架，可用于研究材料和结构的变形和断裂行为，如EMU、peridigm^[19-21]。简要介绍近场动力学的数值实现技术。由于近场动力学的运动方程是积分-微分形式，理论求解难度很大。因此，基于空间离散格式的数值技术用于计算积分方程。在现有的近场动力学模型中，多是考虑到无网格方法空间离散的简单性而被采用。在该框架下，首先将连续体离散成一系列具有一定体积形状的物质点，然后对离散的物质点进行空间积分得到其相应的物理量。

2.1 运动方程离散

作为无网格方法的优点之一，在预处理阶段，即空间离散时不需要网格化。作为一种非局部理论，其本构模型自然包含了损伤和断裂的描述，从而避开了传统有限元失效分析中的不连续性和网格重构。在计算中，近场动力学运动方程离散为每个物质点的一系列离散运动方程：

${\rho _i}{\ddot u}\left[ {{x_i},{t_n}} \right] = \sum\limits_{j = 1}^N {\left( {{{\boldsymbol{T}}}\left[ {{x_i},{t_n}} \right]\left\langle {{\xi _{ij}}} \right\rangle - {{\boldsymbol{T}}}\left[ {{x_j},{t_n}} \right]\left\langle {{\xi _{ji}}} \right\rangle } \right)} \Delta {V_{{x_j}}} + {{\boldsymbol{b}}}\left[ {{x_i},{t_n}} \right]。$

(28)

式中：下标i, j分别为粒子标号；N为物质点近场内粒子总数； ${\xi _{ij}}{\text{ = }}{x_j} - {x_i}$ 。

MPPD的非局部体积应变计算：

${{\boldsymbol{E}}}\left[ {{x_i},{t_n}} \right] = \frac{1}{4}\sum\limits_{j = 1}^N \omega \langle {\xi }\rangle \left( {\frac{{{{\boldsymbol{Y}}}\langle {\mathbf{\xi }}\rangle {{\boldsymbol{Y}}}\langle {\xi }\rangle }}{{{{\xi }^{\text{T}}}{\xi }}} - 1} \right)\left( {5\frac{{{\xi } \otimes {\xi }}}{{{{\xi }^{\text{T}}}{\xi }}} - 1} \right)\Delta {V_{{x_j}}}。$

(29)

2.2 时间积分

利用速度Verlet算法对粒子的速度和位移进行更新：

$\left\{ {\begin{array}{*{20}{l}} {{\dot u}[{x_i},{t_{n + \frac{1}{2}}}] = {u}[{x_i},{t_n}] + {\ddot u}[{x_i},{t_n}] \cdot \dfrac{{\Delta t}}{2}}，\\ {{u}[{x_i},{t_{n + 1}}] = {u}[{x_i},{t_n}] + {\dot u}[{x_i},{t_{n + \frac{1}{2}}}] \cdot \Delta t} ，\\ {{\dot u}[{x_i},{t_{n + 1}}] = {\dot u}[{x_i},{t_{n + \frac{1}{2}}}] + {\ddot u}[{x_i},{t_{n + 1}}] \cdot \dfrac{{\Delta t}}{2}} 。\end{array}} \right.$

(30)

式中，∆t为计算时间步长。

显式格式求解的一个主要局限性是条件稳定，这意味着最大容许时间增量必须小于波传播到最小单元长度所需的时间，所用时间称为稳定时间增量。根据Von Neumann稳定性分析的经典假设，式（28）中使用的时间增量的稳定性条件如下：

$\Delta t < \sqrt {\frac{{2\rho }}{{\sum\limits_{{x_j} \in {{H}}x} {{V_j}} c\left( {{x_j} - {x_i}} \right)}}} 。$

(31)

2.3 数值修正

近场动力学在数值实现中将连续体离散为一定数量具有一定体积形状的物质点，这将导致位于近场范围边界处的物质点体积计算不准确。

由近场动力学理论的非局部特性可知，物质点 ${x_k}$ 与其近场范围内所有的物质点之间存在相互作用，这时需要计算相互作用物质点的体积，容易发现对于处在近场范围边界处的物质点其体积不一定能够处于物质点 ${x_k}$ 的近场范围内。如图2所示，显然有的区域没有完全在近场范围之内，倘若不进行体积修正，在实际计算时将导致体积计算不准确。

为了解决处于近场范围边界处物质点体积计算不准确的问题，一个可行的方法是对处于近场范围边界处的物质点在计算体积时进行特殊处理，即对物质点在计算体积时乘以相应的体积修正系数，称之为体积修正。本文采用的体积修正方案为：

${\alpha }_{{x}_{j}}(\xi )=\left\{\begin{split}&1，\text{if}\; \xi ⩽(\delta -0.5\Delta x)，\\& \dfrac{\delta -\xi +0.5\Delta x}{\Delta x}，\text{if}\;(\delta -0.5\Delta x) < \xi ⩽\delta ，\\ &0，{\rm{others}}。\end{split}\right.$

(32)

考虑体积修正，运动方程（28）写为：

$\begin{split} {\rho _i}{\ddot u}\left[ {{x_i},{t_n}} \right] = & \sum\limits_{j = 1}^N {\left( {{{\boldsymbol{T}}}\left[ {{x_i},{t_n}} \right]\left\langle {{\xi _{ij}}} \right\rangle - {{\boldsymbol{T}}}\left[ {{x_j},{t_n}} \right]\left\langle {{\xi _{ji}}} \right\rangle } \right)} ，\\ & ({\alpha _{{x_j}}}\Delta {V_{{x_j}}}) + {{\boldsymbol{b}}}\left[ {{x_i},{t_n}} \right]。\end{split}$

(33)

3 模型验证 3.1 无限变形算例

采用具有中心孔的矩形板两端受均布载荷来验证当前近场动力学模型，并与Le等^[22]结果对比，如图3所示。矩形板的上下边界采用自由边界，模型的左右两边界分别承受1MPa的均布拉伸载荷作用。模型的尺寸为100 mm×50 mm×10 mm，孔的半径为10 mm。材料属性为：弹性模量为100 GPa，泊松比为0.4，材料密度为2000 kg/m³。

在计算中模型离散粒子尺寸dx=0.025 cm，近场半径 $\delta$ =4×dx。材料采用均匀离散，施加载荷的边界有3层粒子，数值实现采用动态松弛法的显式时间积分方案，时间步长增量取为1。

将提出的模型得到的位移场与现有的常规态基近场动力学模型得到的水平和垂直位移场进行比较，验证当前模型计算的准确性。从图4可看出，本文改进的近场动力学得到的数值结果与常规态基近场动力学模型得到的收敛结果一致。同时也可发现当前模型的收敛速度较Le等^[22]常规态基近场动力学模型慢。这是因为当前模型在计算键力时一方面需多次计算粒子近场的积分解；另一方面本文采用较复杂的复合幂型指数函数，这些都增加了计算量。通过比较最终收敛后的结果，可验证本文模型的正确性和准确性。

图5和图6给出了A和B 2个点位移和速度的定量对比结果，其中黑色方块代表本文模型计算数据，黑色三角形代表Le等^[22]的计算结果。

由图5和图6可以发现，垂直和水平位移在大约1 000次迭代时达到稳态值，表明了本文模型准静态动力学问题的模拟是有效的。

3.2 有限变形算例

MPPD模拟的有限变形采用悬臂梁一侧受集中力另一侧固定来研究。悬臂梁的几何尺寸如图7所示。MPPD模型采用自适应动态松弛方法，有效解决准静态问题模拟。选取3个沿悬臂梁上侧面长度方向均匀分布的粒子P₁、P₂、P₃，研究其位移。力学参数为：弹性模量100.0 MPa，泊松比0.3，材料密度7800 kg/m³。长度L=10.0 m，高度H=1.0 m，集中力P=1 MPa。

图8给出了采用当前MPPD模型模拟悬臂梁变形过程的垂直和水平方向云图。由图8可以发现，当前MPPD模型能够比较稳定地模拟有限变形过程。

图9给出了采用MPPD，有限元和Le等模型计算的5个标记点在竖直方向和水平方向的位移对比结果。由图9可以发现，当前MPPD模型在垂直和水平位移均与有限元结果较为吻合，从而说明当前MPPD模型能够有效模拟有限变形问题。

4 结　语

本文推导了采用微势能的近场动力学本构模型，标定了平面应力问题的键参数，给出了数值实现的离散格式，并考虑了采用粒子离散时的数值修正因素，最后分析了无限变形和有限变形问题。结论如下：

1）采用二阶变形作为键长变形的度量，结合指数-指数微势函数推导的微势近场动力学模型能够有效模拟材料变形问题。

2）MPPD模型可以同时模拟无限变形和有限变形,与现有常规态基近场动力学模型相比，MPPD模型能够获得较为准确的结果。