近场动力学数值模拟的程序实现（1）——引言

引 言

近场动力学已为在计算模拟代码中对材料破坏和破碎建模提供了最新技术。为了有效地对真实世界的系统进行建模，基础理论必须伴随有效的软件实现。本章的目的是提供一个在分析代码中实现近场动力学的路线图。这部分是由于存在于诸多有关近场动力学的文献之间的鸿沟而引起的，这通常侧重于理论考虑，而在实现中可能出现困难。

本章专注于介绍Silling and Askari ^[1] 的无网格方法。迄今为止，绝大多数近场动力学模拟都采用了这种方法。专注于此特定计算策略的决定，并不意味着它优于其它可替代方法。例如，在某些情况下，将有限元技术应用于近场动力学可能会提供更好的结果。尽管如此，Silling and Askari的无网格方法已被证明是一种解决固体力学中普遍存在的材料破坏和断裂问题的可靠而有效的策略。其优势包括易于实现、计算效率高以及存在用于材料分离的自然机制。Silling and Askari的方法是线性动量近场动力学平衡的强形式的直接离散化 ^[2] ：

\rho(\mathbf{x}) \ddot{\mathbf{u}}(\mathbf{x}, t)=\int_{\mathcal{H}_{\mathbf{x}}}(\underline{\mathbf{T}}[\mathbf{x}, t]\langle\mathbf{q}-\mathbf{x}\rangle-\underline{\mathbf{T}}[\mathbf{q}, t]\langle\mathbf{x}-\mathbf{q}\rangle) d V_{\mathbf{q}}+\mathbf{b}(\mathbf{x}, t)\quad\quad\quad(1)

其中 \rho 是材料密度， \ddot{\mathbf{u}} 是加速度， t 是时间， \mathbf{x} 和 \mathbf{q} 是材料点， \mathcal{H}_{\mathbf{x}} 是具有以 \mathbf{x} 为中心的半径 \delta 的球形邻域， \underline{\mathbf{T}} 表示近场动力学力态（force state）， d V_{\mathbf{q}} 表示与 \mathbf{q} 相关的无穷小体积，而 \mathbf{b} 是体力密度。对于三维空间中的问题，Silling and Askari的方法要求估计三维积分。相比之下，例如使用有限元离散化来求解线性动量的近场动力学平衡的弱形式需要估计六维积分。因此，相对于强形式策略，基于弱形式的近场动力学模型的求解策略需要相当大的几何复杂度和计算开销。

近场动力学模拟的目标是在指定的条件下、在指定的时间段内，对一个物体或一组物体中的材料响应进行建模。这是通过离散化时空中的控制方程来以数值方式实现的。根据Silling and Askari ^[1] ，物理系统被离散化为无限数量的节点，每个节点都具有与其关联的（标量）体积。采用了共置位（co-locational）方法，这意味着节点用于跟踪场变量（例如位移和速度），还用作本构模型评估的材料点。该策略可与 单点高斯积分格式 相比较，在该格式中，假定所有场在积分域上都是恒定的，而高斯点位于该域的形心处（类似于中点积分）。

对于离散化中的每个材料点，通过求解以下随时间变化的半离散运动方程，可以确定整个近场动力学系统的响应：

\rho(\mathbf{x}) \ddot{\mathbf{u}}(\mathbf{x}, t)=\sum_{\mathcal{H}_{\mathbf{x}}}(\underline{\mathbf{T}}[\mathbf{x}, t]\langle\mathbf{q}-\mathbf{x}\rangle-\underline{\mathbf{T}}[\mathbf{q}, t]\langle\mathbf{x}-\mathbf{q}\rangle) \Delta V_{\mathbf{q}}+\mathbf{b}(\mathbf{x}, t)\quad\quad\quad(2)

其中式（1）中的积分已由材料点 \mathbf{x} 的邻域的集合求和代替。力态 \underline{\mathbf{T}}[\mathbf{x}, t] 和 \underline{\mathbf{T}}[\mathbf{q}, t] 分别通过在点 \mathbf{x} 和 \mathbf{q} 评估本构模型来决定。接触相互作用（如果有的话）提供了附加的节点力密度，必须将其包含在式（2）中。通常，模拟可能涉及多个体，多种材料以及任意数量的指定初始条件和边界条件。式（2）对于键基本构模型和态基本构模型均成立，因为键基模型可以被视为更一般的态基模型族的子集。

模拟的进度由时间积分程序驱动。模拟的持续时间被分为有限数量的时间步长或负载步长，在该时间步长或负载步长上计算式（2）。在显式时间积分格式的情况下，式（2）用于直接求解节点加速度。然后将加速度应用于每个节点，以使仿真从一个时间步前进到下一个时间步。相比之下，用于准静态模拟的隐式时间积分格式将节点位移视为未知量。应用式（2）来求解那些不受运动学边界条件约束的节点自由度产生零加速度的位移值。无论采用何种时间积分策略，仿真都将通过在给定的时间步长规定一组已知值（例如，在域子集上的位置或速度）来进行，然后通过应用控制近场动力学方程，来确定剩余的场和材料状态数据。最终结果是对模拟中每个时间步的离散化中所有材料点（节点）的运动和材料状态数据的完整描述，从中可以确定任何数量的相关次级量。

在较高的层次上，计算近场动力学代码包含几个关键组成部分：时间积分器、评估内力的方法、评估键破坏定律的流程、施加初始条件和边界条件的流程、评估和施加接触力的方法，以及用于存储和检索数据的输入/输出例程。对于瞬态动态仿真，必须估算关键时间步。采用隐式时间积分的仿真通常需要实现非线性求解器。牛顿法是求解非线性方程组的一种标准方法，该方法要求构造正切刚度矩阵和一套求解线性方程组的方法。

近场动力学已在独立软件中实现，并且还作为更大、更全面的仿真代码的组成部分。 由Stewart Silling编写的EMU是第一个计算近场动力学代码 ^[1] 。先前的实现近场动力学的软件包实例包括LAMMPS和Sierra/SolidMechanics ^{[3] [4]} 。人们还证明了可以在商业有限元代码Abaqus中实现近场动力学 ^[5] 。近场动力学的最近实现是开源代码Peridigm ^{[6] [7]} 。

关于应用于近场动力学离散化的命名法存在混淆的可能。在本章中，术语“节点（node）”，“节点体积（node volume）”和“材料点（material point）”用于描述离散的近场动力学体。术语“节点体积”直接指Silling and Askari的离散化方法中使用的逐点位置和体积数据。节点既用作追踪模型几何的手段，又用作材料模型评估的位置。术语“节点”和“材料点”用于指代这些不同的作用，并与标准计算力学术语保持一致。Silling and Askari的方法通常称为 无网格（meshless/meshfree） ，这是因为没有考虑定义物理体积单元的节点和节点连接性（如在标准有限元方法中一样）。关于此术语的混淆主要源于以下事实：实际上，在离散化过程中经常使用定义体积单元的几何数据作为计算节点体积的方法（比较Henke and Shanbhag ^[8] ）。此外，与每个节点相关联的区域的几何形状，已被用于通过考虑这些实体与定义近场动力学邻域的球体的部分相交，来提高近场动力学求积的精度 ^{[9] [10] [11]} 。

以下各节介绍计算近场动力学代码的主要组成部分。在这样做的说明中，以Peridigm代码为例。讨论假定模拟是三维的，尽管概念通常也适用于一维和二维模拟。大部分内容都集中在串行代码上。尽管并行代码的实现增加了复杂性，但核心近场动力学算法并不需要根本性的改变。首先介绍了近场动力学本构模型的软件实现，然后介绍了粘结破坏、正切刚度矩阵、接触建模、无网格离散化、用于识别近场动力学粘结的邻近搜索和时间积分等部分。最后，给出了示例仿真，展示了显式瞬态动力学和准静态特性。

封面： ^[12]

参考

1. ^ ^{a b c} S. A. Silling and E. Askari. A meshfree method based on the peridynamic model of solid mechanics. Computers and Structures, 83: 1526–1535, 2005.
2. ^ S. A. Silling, M. Epton, O. Weckner, J. Xu, and E. Askari. Peridynamic states and constitutive modeling. Journal of Elasticity, 88: 151–184, 2007.
3. ^ M. L. Parks, R. B. Lehoucq, S. J. Plimpton, and S. A. Silling. Implementing peridynamics within a molecular dynamics code. Computer Physics Commu- nications, 179(11): 777–783, 2008.
4. ^ SIERRA Solid Mechanics Team. Sierra/SolidMechanics 4.40 User’s Guide. Report SAND2016-2707, Sandia National Laboratories, Albuquerque, NM and Livermore, CA, 2016.
5. ^ R. W. Macek and S. A. Silling. Peridynamics via finite element analysis. Finite Elements in Analysis and Design, 43: 1169–1178, 2007.
6. ^ Peridigm peridynamics code. http://peridigm.sandia.gov
7. ^ M. L. Parks, D. J. Littlewood, J. A. Mitchell, and S. A. Silling. Peridigm Users’ Guide v1.0.0. Report SAND2012-7800, Sandia National Laboratories, Albuquerque, NM and Livermore, CA, 2012.
8. ^ S. F. Henke and S. Shanbhag. Mesh sensitivity in peridynamic simulations. Computer Physics Communications, 185: 181–193, 2014.
9. ^ F. Bobaru and D. Y . Ha. Adaptive refinement and multiscale modeling in 2D peridynamics. International Journal for Multiscale Computational Engineer- ing, 9(6): 635–660, 2011.
10. ^ P . D. Seleson. Improved one-point quadrature algorithms for two-dimensional peridynamic models based on analytical calculations. Computer Methods in Applied Mechanics and Engineering, 282: 184–217, 2014.
11. ^ P . D. Seleson and D. J. Littlewood. Convergence studies in meshfree peridy- namic simulations. Computers and Mathematics with Applications, 77(11): 2432–2448, 2016.
12. ^ https://timgsa.baidu.com/timg?image&quality=80&size=b9999_10000&sec=1597687214612&di=0acfad12dccacd6f87e8e83430250328&imgtype=0&src=http%3A%2F%2Fyouimg1.c-ctrip.com%2Ftarget%2F100a0m000000dz9uxC264.jpg