周围动力学(PD)模型通常通过利用称为标准方案的基于粒子的方法来实现。与基于经典理论的数值方法(如有限元方法)相比，使用无网格标准方案的PD模型通常计算成本更高，主要有两个原因。首先，PD的非局域性需要先进的正交方案。其次，标准格式的非均匀离散化是不准确的，因此通常可以避免。因此，非常精细的均匀离散化应用于整个域，即使是在只有一小部分(例如，接近不连续面和界面)本身需要精细分辨率的情况下。在本研究中，设计了一个新的框架，大大提高了PD模型的计算性能。它仅将标准格式应用于出现不连续和界面的局部区域，并将要求较低的正交格式应用于该区域的其余部分。此外，该方法采用多网格方法，仅在关键区域采用精细的网格间距。由于这些区域是随时间动态识别的，因此我们的框架被称为多自适应框架。通过两个现实世界的问题，Kalthoff-Winkler实验和镁基骨植入螺钉的生物降解来检验所提出的方法的性能。结果表明，与标准方案的简单应用相比，我们的新框架可以大大降低计算成本(对于给定的精度要求)。

　　材料失效和断裂建模是计算力学研究人员和实践者面临的一个主要问题(Ravi-Chandar 2004;安德森2017年)。在过去的几年里，基于经典连续介质力学(CCM)的几种方法被用于这一目的。然而，它们的应用引入了一些困难，这些困难与控制方程中存在位移的空间导数有关，当位移场的连续性没有得到验证时，这些导数是未定义的(Silling 2000;Bobaru et al. 2016)。因此，许多科学家试图使这些方法具备模拟裂纹形成和扩展的能力(Kuhn and m

　　ller 2008;Amor et al. 2009;Mo?s et al. 1999;Belytschko和Black 1999;Ortiz et al. 1987;Pepe等人，2020;Reccia et al. 2018)，但所有提出的策略都存在一些缺点(Bobaru et al. 2016)。

　　近年来，基于周期动力学(PD)理论的创新方法(Silling 2000;Silling等人(2007)提出了基于积分-微分方程的经典连续介质力学的非局部重新表述，以解决基于ccm的方法在处理不连续方面的缺点。

　　PD理论处理的是积分方程而不是空间微分。因此，即使存在不连续的位移场，PD控制方程也是有效的，允许将裂纹的成核、扩展和分支等复杂现象视为自然的材料响应。

　　空间积分在局部局部模型的数值实现中起着至关重要的作用，因此积分的精确计算对局部局部模型的数值实现具有重要意义。

　　PD模型的大多数实现采用了Silling和Askari(2005)中引入的基于粒子的方法，这里称为标准离散化方案，它是无网格的，因为它不依赖于物体离散化所需的底层节点网格之间的任何元素或几何连接。

　　近年来，减少局部放电模型的计算负担一直是各种研究的主题。事实上，一系列的工作都试图通过高性能的计算方法来增强PD模型的计算机实现(Mossaiby et al. 2017;Diehl et al. 2020;男孩等人2021;Brothers et al. 2014;Fan and Li 2017)，而其他研究提出了自适应网格细化算法来优化计算资源的使用(Bobaru et al. 2009;Bobaru and Ha 2011;Dipasquale et al. 2014;Shojaei et al. 2018;Bazazzadeh et al. 2020;Gu et al. 2017)。另一个研究分支将注意力集中在将基于pd的模型与基于CCM的方法相结合的方法的发展上，最常见的是有限元法(FEM) (Oterkus et al. 2012;Han and Lubineau 2012;Lubineau et al. 2012;Yu et al. 2018;太阳和鱼2019;Galvanetto et al. 2016;怀尔德曼和Gazonas 2014;Shojaei et al. 2016;Zaccariotto et al. 2018;Ongaro等，2021,2022,2023;Seleson et al. 2013;D 'Elia et al. 2021;Wang et al. 2019b;Seleson et al. 2015)。尽管这些方法被证明是有效的，但它们都存在一些局限性(Mossaiby et al. 2022;ongaro et al. 2021;Zaccariotto et al. 2018)，这主要是由于耦合方法的色散关系不同(Shojaei et al. 2023;Hermann et al. 2023)。

　　在Shojaei et al.(2022)中，为了提高PD模型的数值性能，作者引入了一种有效的无网格配置方案来离散弹性和扩散问题中的PD控制方程。正如Shojaei等人(2022)所全面讨论的那样，尽管所提出的配置方案在计算成本和空间积分计算方面具有优势，但在关键区域(即以存在不连续或具有复杂几何细节为特征的区域)附近可能会受到精度损失的影响，而标准方案在这些区域的表现更好。因此，在Shojaei等人(2022)中，作者引入了一种混合无网格离散化，通过将标准方案和配置方案结合在一个自适应框架内，利用了它们的优点。在混合构型中，标准方案的应用实际上仅限于局部的关键区域和边界，而身体的其余区域则通过更有效的搭配方案进行离散化。

　　在本工作中，为了进一步提高PD模型的数值性能，将混合离散化的使用与Shojaei等人(2018)提出的多网格方法的实现相结合，其中通过耦合不同间距的网格来实现具有可变网格尺寸特征的PD模型。这样，精细网格的使用仅限于需要精细网格的区域的一小部分(即关键区域和边界，因为沿着边界采用更精细的网格实际上是减少众所周知的PD模型表面效应的有效方法)，而区域的其余部分则通过更粗的网格进行离散化。因此，在本研究中，关键区域和边界采用标准方案和细网格进行离散化，其余区域采用搭配方案和粗网格进行离散化。所提出的方法被称为多自适应，因为受Shojaei等人(2016)的启发，设计了一种自适应算法，可以动态切换配置方案所代表的区域的离散化方案和网格间距，并将粗网格离散化为标准方案和细网格间距。通过该算法，采用标准方案和精细网格离散的零件可以随时间演化，并跟踪关键区域(如动态断裂问题中出现的裂纹形态)的演化，同时相应地更新耦合配置。

　　这导致了一个高性能和高效的纯PD模型，该模型最佳地利用了计算资源。

　　通过一些数值例子，包括二维和三维的脆性断裂和腐蚀问题，对该方法的能力进行了评估。通过将所获得的解与由标准方案表示的PD模型生成的解进行比较，突出了所开发框架的高性能，标准方案分别采用均匀粗离散化、均匀精细离散化和基于上述多网格方法的非均匀离散化(Shojaei et al. 2018)。

　　目前的工作是将Shojaei等人(2022)的混合方案扩展到空间中的非均匀离散化。与Shojaei等人(2022)提出的方法不同，这项工作更侧重于证明将所提出的策略应用于研究大型复杂结构的可行性，这些结构需要大的水平尺寸，并且受到域的某些局部区域存在不连续和界面的影响，需要更精细的分辨率。另一个创新方面在于将该方法应用于具有物理参数的现实世界生物降解示例。

　　本文的内容组织如下:在第2节中，简要概述了用于脆性断裂和扩散问题的基于键的PD公式。第3节对Shojaei et al.(2022)中介绍并在本研究中使用的标准方案和搭配方案进行简要总结。在本节的最后一部分，概述了多自适应方案的实现，重点描述了Shojaei等人(2018)提出的多网格方法，以及Shojaei等人(2016)引入的开关策略，并在本工作中加以利用。在第4节中，对所提出的多自适应方案的性能进行了数值研究，其中包括第4.1节中二维模型中的动态脆性断裂和第4.2节中的三维腐蚀问题。第5节是一些结束语。

　　图1

　　

　　变形前后一般PD域的表示;两个材料点和之间的相对位置向量(初始和当前)和相对位移向量也被报道

　　在空间维度为p的域中，用PD模型描述，每个材料点与位于该材料点的有限邻域内的所有其他材料点相互作用(见图1)。Silling(2000)给出了任意材料点在某一时刻的基于键的PD运动方程:

　　（1）

　　式中为质量密度，为位移场的时间二阶导数，为成对力函数，表示质点对质点施加的单位体积力的平方(或微力密度)，为规定的体力密度场。该社区的定义如下:

　　（2）

　　地平线在哪里?对于物体体积中的质点，即远离物体边界的质点，邻域表示一维的线段，二维的圆盘，三维的圆心为的球。

　　参考构型中两个质点的相对位置矢量(或初始相对位置矢量)表示为(见图1):

　　（3）

　　它代表了键的标准周动力学符号。

　　在变形构型时，两个材料点和将分别被和移位。如图1所示，对应的相对位移矢量定义为。

　　力矢量，也称为粘结力，作用于连接两个质点的线的方向，即在变形构型中它们的相对位置矢量的方向，通常称为当前相对位置矢量(见图1)。

　　包含材料全部本构信息的成对力函数的一般形式为:

　　（4）

　　其中是基于材料类型定义的标量值偶函数。

　　以Silling和Askari(2005)中引入的原型微弹性脆性(PMB)材料模型为例，为:

　　（5）

　　式中，c为微模量函数，表示粘结弹性刚度;s为粘结拉伸，表示为:

　　（6）

　　在变形较小的情况下，考虑到(5)，假设(4)的线性化形式为:

　　（7）

　　可以将微模量函数c与可测量的宏观量联系起来，例如材料的拉伸模量E和泊松比。

　　有必要强调的是，由于基于键的PD公式，其中键仅基于成对相互作用来表征，泊松比被限制为固定值。对于三维和二维平面应变情况，泊松比固定为，而对于二维平面应力情况，泊松比被限制为(Silling 2000;Gerstle et al. 2005)。在Silling等人(2007)引入的基于状态的理论版本中，这一限制已被消除。Silling和Askari(2005)、Gerstle等人(2005)、Bobaru等人(2009)、Ha和Bobaru(2010)全面概述了用上述宏观量获得粘结弹性刚度的过程。在本工作中，采用恒定微模函数。

　　在PMB材料模型中，可以通过为键拉伸建立一个预定义的极限值(通常称为临界拉伸)来引入失效，并考虑当当前拉伸超过该极限值时键断裂(Silling和Askari 2005)。因此，所采用的破坏准则称为最大拉伸准则。

　　对于该材料模型，(4)可改写为(Silling 2000):

　　（8）

　　式中为依赖历史的标值函数，作为断键参数引入，因此可以取以下两个值之一:

　　（9）

　　临界拉伸可以与可测量的宏观量有关，如材料的临界能量释放率(见Silling和Askari(2005))。

　　感兴趣的读者可以参考Dipasquale等人(2022)确定基于状态的PD的临界拉伸。本构模型的历史依赖性是一个事实的结果，即一个键的断裂是一个不可逆的过程。

　　然后，引入局部损伤指数来表示时间t点处的损伤概念，该局部损伤指数由下式定义:

　　（10）

　　式中，0表示材料的完好状态，1表示该材料点与其周边所有材料点完全断开(Silling and Askari 2005)。

　　周期动力学的非局域性引入了一个关于边界条件应用的重要问题。考虑到PD平衡方程是基于积分算符而不是偏微分算符，它们的变分性质不会像经典连续介质力学模型那样产生自然(诺伊曼)边界条件(Macek and Silling 2007)。因此，在用PD模型描述的域中，规定的位移或载荷需要通过有限体积的边界层施加，而不是在表面上施加。根据Macek和Silling(2007)的数值研究，沿物质边界的边界层深度通常设为等于。

　　在下一节中，简要概述了扩散型问题的基于键的PD的非局部公式。在Bobaru和Duangpanya(2010年，2012年)中提出了该配方的全面推导。该模型是通过将浓度值关联到时间t的每个物质点来构建的。每个物质点通过通常称为扩散键的管道状导体与位于有限邻域内的所有物质点相连。在这个框架中，每个键可以被认为是一个管道，它在两个桶之间传递浓度，即两个材料点和两端。在基于键的PD扩散公式中，即本研究中使用的PD版本，每个键中浓度的输运独立于其他键中的输运(Bobaru and Duangpanya 2012)。

　　Bobaru和Duangpanya(2010)给出了任意时间点基于键的PD扩散控制方程:

　　（11）

　　式中，C表示浓度场，为C的一阶导数，S为给定源函数，J为积分算子或浓度流密度的核，其形式为(Zhao et al. 2018):

　　（12）

　　其中为视界(见2.1节)，q是影响核J形状的整数，通常取为等于0、1或2(详见Chen and Bobaru (2015b))，称为PD微扩散系数。提出了几种形式。在这项工作中，采用恒定的微扩散系数(Zhao et al. 2018)。Bobaru和Duangpanya(2010)将局部扩散通量等同于局部扩散通量，从而将局部扩散系数与经典局部扩散方程的扩散系数联系起来;Shojaei et al.(2020)。在二维扩散中，假设，可表示为:

　　（13）

　　与弹性情况2.1节中介绍的类似，在PD扩散中，边界条件需要通过有限体积的边界层来应用，其深度通常取为。(11)中源函数S的一部分包含neumann型边界条件(参见Chen and Bobaru (2015b)， Wang et al. (2019a))。

　　摘要。

　　1 介绍

　　2 bond-ba概述

　　sed peridynamics

　　3 标准和建议的离散化方案

　　4 数值例子

　　5 结论

　　参考文献。

　　致谢。

　　作者信息

　　道德声明

　　相关的内容

　　# # # # #

　　图2

　　

　　标准离散化方案中与节点相关的邻域的示意图

　　可以采用不同的数值方法来实现周动力模型的数值离散化(Emmrich and Weckner 2007)。用于离散PD控制方程的最常用方法是Silling和Askari(2005)引入的无网格方案，此处称为标准离散化方案(Shojaei et al. 2022)。在这个框架中，周期动力学方程的数值近似需要将域细分为称为节点的点网格:在用PD模型描述的域中，每个节点在参考构型中与某个有限体积相关联，使得节点集形成一个均匀的网格，其中最近的两个相邻节点之间的距离称为网格间距(见图2)。在均匀网格构型中，，其中，和分别为x方向、y方向和z方向的网格间距。没有元素或节点之间的其他几何连接被预见(Silling和Askari 2005)。与每个节点相关联的体积，即节点体积，对于三维模型计算为立方体体积，对于平面情况计算为正方形体积单元，对于一维系统计算为正方形体积单元，其中和a分别为二维和一维域的厚度和恒定横截面积。卷位于中心的相关节点称为源节点。并且，将时间离散为瞬间，即，其中n表示时间步长数。

　　在标准离散化方案中，任意源节点在时间瞬间的基于键的运动PD方程(参见(1))的离散化形式是将(1)中的积分替换为有限和，如下所示:

　　（14）

　　同样，将(11)中的积分代入有限和，得到(11)中PD扩散方程(cf.(11))的离散化形式，有:

　　（15）

　　式(14)和(15)中，i为源节点的索引，j为包含节点有限邻域内所有节点的族节点集的关联索引，为节点关联的离散化体积，为体积校正因子，用于逼近源节点邻域内部分体积(见图2)。对于二维情况，采用如下形式(Yu et al. 2011):

　　（16）

　　采用一点高斯正交规则进行空间积分。在本工作中，Yu et al.(2011)提出的算法在目前可用的提高空间整合精度的方法中被选择，因为它是PD社区中最常用的算法，因为它易于实现且结果相当准确。Seleson(2014)提出的正交算法和Scabbia et al.(2022)引入的创新算法当然可以交替使用，也许会导致稍微更准确的结果。然而，考虑到工作的重点更多地是效率而不是计算的准确性，所选择的算法为所考虑的问题类型提供了足够准确的解决方案。

　　离散化PD的一个基本参数是地平线与网格间距的比值，即;因此，和m是决定离散化PD模型中每个节点需要考虑的相互作用数量的两个参数。

　　在目前的工作中，弹性问题的时间积分是通过利用速度- verlet时间积分(Swope et al. 1982)来完成的。该显式算法具有鲁棒性和数值稳定性，在PD框架中得到广泛应用。

　　采用显式正演欧拉算法代替扩散问题的时间积分。给定时刻瞬间的节点值，即和，即可实现浓度场的更新过程:

　　(17)

　　出于数值稳定性的考虑，模拟过程中采用的时间步长必须小于预定义的临界时间步长。这一关键步长的计算一直是文献中许多研究的主题，例如(Silling和Askari 2005;Bobaru et al. 2016;2015年Littlewood)。在Silling和Askari(2005)中，作者对关键时间步的评估如下:

　　(18）

　　其中j遍历源节点的所有族节点。这种稳定性条件对于基于键的PD模型特别有用(Silling和Askari 2005)。

　　计算临界步长的另一种选择是考虑众所周知的Courant-Friedrichs-Lewy条件，例如(Bobaru et al. 2016;2015年Littlewood):

　　（19）

　　式中为离散域的最小网格尺寸，为纵波在材料中的传播速度，用关系式表示:

　　（20）

　　其中K为材料的体积模量。

　　然而，后一种情况对于PD模型来说是相当保守的(Bobaru et al. 2016);因此，在本工作中，考虑Silling和Askari(2005)中导出的稳定性条件来执行临界时间步长计算。

　　标准离散化方案易于实现，可以方便地计算非线性PD模型所描述的包含物质分离(如裂纹)的区域所需的积分项。

　　然而，一旦取一个大值，它可能会变得计算昂贵;在三维问题中，当需要大的水平尺寸和精细的网格间距时，例如在裂纹尖端区域或腐蚀界面附近的对应关系中，这个问题会加剧。

　　这是因为在标准方案中，所有相邻节点都充当积分的正交点。在Shojaei et al.(2022)中，为了提高PD解算器的数值性能，建议限制标准格式的应用，从而在不一定需要应用标准格式的区域采用有效的无网格配置格式。下一节将描述并置方案。

　　图3

　　

　　与节点相关的邻域示意图，其中(左)为标准方案中的族节点及其对应的积分单元，(中)为被笛卡尔辅助网格包围的族节点，(右)为搭配方案中选择的搭配节点

　　该搭配方案是真正的无网格的，因为它以半解析的方式计算积分，并且不像标准方案那样需要积分的背景单元。在该方案中，在积分之前，在每个邻域处，主要的场变量用一组节点的节点值来近似，这些节点被称为搭配节点，用表示，它们是族节点的子集(即)。图3(左)显示了一个用标准方案离散化的通用PD邻域。受Shojaei et al.(2022)的启发，我们考虑了一个以源节点为中心的辅助节点网格，其笛卡尔模式如图3(中)所示。辅助网格表示搭配节点的近似位置(即搭配方案中积分的正交点)。同样地，对每个辅助节点，关联最近的家族节点，从而将其加入到搭配节点集合中;见图3(右)。通过加权最小二乘(WLS)方法逼近，该方法使用单项式作为基函数。根据Shojaei et al.(2022)，搭配节点的数量取决于近似的阶数。为了使计算有效率，采用尽可能少的基数;例如，在三维空间中:

　　（21）

　　式中是一个包含单项式基函数的向量，在这种情况下，它表示一组完整的直至二阶的单项式。在对场变量进行近似后，在具有满邻域的节点上，控制方程可以用线性PD模型描述(参见(7))，相应的积分可以根据搭配节点的节点值解析计算。由此得到节点处的运动方程如下:

　　（22)

　　其中集合时刻沿、、和方向的位移:

　　(23)

　　式(22)中，该项产生正交权值;为矩阵，为(Shojaei et al. 2022):

　　(24)

　　地点:

　　（25）

　　是矩阵，为:

　　（26）

　　其中通过:

　　（27）

　　其中，和是获得WLS近似中的矩矩阵所需的两个矩阵，为(Shojaei et al. 2017,2016):

　　（28）

　　和:

　　（29）

　　在(28)和(29)中，是影响每个搭配节点相对于其与源节点的距离的贡献的权重函数，在这里取为(Shojaei et al. 2017;Mossaiby et al. 2020;Shojaei et al. 2016):

　　(30)

　　其中表示搭配节点到源节点的最大可能距离。对于(11)中的扩散方程，可以采用类似的策略进行离散化，得到内近似的下式:

　　(31)

　　式中为时刻节点的集中度值:

　　(32)

　　表示得到的向量如下:

　　(33)

　　因此，与标准方案相比，该方案不需要计算任何部分体积，并且在更少的节点上进行积分。读者可以参考Shojaei et al.(2022)来了解搭配方案的收敛性和准确性。

　　如第1节所述，Shojaei等人(2022)引入的混合离散化的使用在这里与Shojaei等人(2018)开发的多网格方法的实现相结合，目的是最大化PD模型的数值性能。因此，在下面的章节中，将概述由此产生的框架(称为多自适应方案)的实施，重点描述用于耦合不同网格大小的网格的多网格方法(Shojaei et al. 2018)以及用于自适应跟踪关键区域演变的切换策略(Shojaei et al. 2016)。

　　3.3.1 不同网格大小的网格耦合

　　图4

　　

　　关键区域附近虚拟物质点的示意图表示和对应活动域相应数据的计算

　　图5

　　

　　域的示意图;(左)细网格，(右)粗网格。为清晰起见，在域的精细区域采用，而在域的粗糙部分考虑，其中为地平线。注意，这里只作为一个例子

　　多自适应方案可以通过考虑包含一个或多个不连续点的一般PD域来提出，这些不连续点可能由于变形过程而传播、分支或合并。在提出的框架中，使用标准方案和精细网格对不连续点周围的边界和关键区域进行离散化，而通过搭配方案和较粗网格对区域的其余部分进行离散化(参见第1节)。

　　区域的粗、细两部分分别称为和，需要进行耦合，以确保关键区域和边界附近的材料点具有完整的邻域。

　　Shojaei等人(2018)引入的多网格方法背后的思想是考虑一系列虚构的材料点，包括在和部分中，其中靠近关键区域边界的材料点和边界期望完整的邻域。

　　需要强调的是，虚拟材料点的功能是完成其他材料点的邻域，从而使粗网格和细网格相互独立。

　　如图4所示，网格的独立性是通过考虑两个相似但不同的域来实现的，一个是通过标准方案和细网格离散化的PD模型，另一个是通过搭配方案表示的PD模型，用粗网格离散化。

　　参照图4，为属于域的物质点，即属于细网格间距离散的域，为细网格中虚构的物质点，对应于的缺失邻居。

　　和之间的相互作用是在与相关的量(如位移、速度、加速度)的计算之后建立的。鉴于这些数据仅在域中可用(见图4)，它们是通过在与in相邻的材料点上插值所需数量的值来计算的。与粗网格中虚拟材料点相关的数量同样是通过插值中相邻材料点的相应数量的值来计算的。

　　这就造成了区域之间的相互约束，因为图4中临界区域外部分的解被约束在其对应的区域中，而临界区域对应的解被约束在其对应的区域中。

　　如图4所示，区域分为活动区域和非活动区域，其中非活动区域表示解中不再需要的约束部分，因此在计算中不考虑。这种结构的优点体现在可以保留系统的总质量，并且可以独立地提前离散每个域，从而提高了实现自适应细化过程的效率，以跟踪关键区域的演变。多网格方案的离散形式可以很容易地引入，通过考虑图4所示的和两个域中靠近临界区域边界的区域的相同部分，如图5所示。在下面，和分别表示用于离散化和域的网格间距。在离散域内识别出四种不同类型的节点，即细网格中的主动和虚拟节点和粗网格中的主动和虚拟节点。如图5所示，细网格中活动节点“的邻居”由活动节点和虚拟节点组成，分别用“和”表示。类似地，粗网格中活动节点的邻居由活动节点和虚构节点组成，分别标记为“和”。

　　下面，在弹性问题中概述细网格中一般活动节点”的时间推进方案，如下:

　　(34a) (34b) (34c)

　　在扩散问题中，这样:

　　(35a) (35b)

　　在(34b)和(35b)中，右侧的圆点表示活动节点”附近其他活动节点和虚构节点的贡献，类似于已经提供的“和”。

　　同样，在弹性问题中定义粗网格中一般活动节点的时间推进方案，有:

　　(36a) (36b) (36c)

　　在扩散问题中，这样:

　　(37a) (37b)

　　在(36b)和(37b)中，右侧的圆点表示活动节点”附近其他活动节点和虚构节点的贡献，类似于已经提供的“和”。

　　根据时间推进格式的定义，在弹性和扩散问题中分别提出了用于评估粗网格和细网格中虚拟节点相关数量的插值格式。如图5所示，“”表示粗网格中的一个虚拟节点，其关联值(例如位移、速度、加速度)是通过对其周围细网格中活动节点的对应值进行逆距离加权的简单插值得到的，例如:

　　(38a) (38b)

　　在弹性问题上，如:

　　(39)

　　在扩散问题中。

　　图6

　　

　　关键区域对应切换策略的实现示意图;(上)时间瞬间和(下)时间瞬间的配置

　　相应的，如图5所示，“表示细网格中的一个虚拟节点，其所需值是通过对其周围粗网格中活动节点的对应值进行逆距离加权的简单插值得到的。”，例如:

　　(40 a) (40 b)

　　在弹性问题上，如:

　　(41)

　　在扩散问题中。

　　3.3.2 切换策略

　　在本研究中，通过实现一种自适应算法来优化计算资源的使用，该算法能够跟踪关键区域(例如，损伤的时间形态)的演变，通过动态替换由配置方案表示的区域，用粗网格离散与标准方案描述的区域并用细网格离散。回顾上一节介绍的域，该策略的主要目的是自适应地在不连续点可能合并或传播的活动区域切换配置，以确保不连续点的合并和传播发生在该域的活动区域内。因此，自适应切换是通过取消激活中的节点，从而相应地激活中的节点来执行的，同时在两个域中识别适当的虚拟节点(见图6)。识别需要切换配置的关键区域，需要定义一个合适的标准来触发自适应算法。受Shojaei等人(2016)、Shojaei等人(2022)、Mossaiby等人(2022)、Sheikhbahaei等人(2023)的启发，本研究中采用的标准是基于每个键的拉伸状态。回顾式(6)中键张力的定义，考虑两个一般相邻节点，在系统的初始配置或参考配置中，在某一时刻，如果两个节点之间的连接键的张力值在以下范围内，则认为该连接键是临界的:

　　(42)

　　式中是表示切换准则保守程度的因子，通常选择在0.9左右(Mossaiby et al. 2022)，而为临界拉伸，如第2.1节所定义。如果连接节点和的键被识别为临界，则将其邻域的联合视为临界区域，从而将用于描述该区域的离散化方案和网格间距分别切换为标准方案和精细网格间距。

　　切换策略的离散形式如图6所示，以便更清楚地了解算法的工作原理。图6的顶部描述了通用域在时间瞬间的关键区域及其附近的结构。在该框架中，一旦自适应算法识别出一对节点之间和粗网格中的连接键为关键(参见(42))，在继续进行下一个时间瞬间之前，需要更新两者的配置。如图6所示，这是通过使位于交换区域内的节点失活，即在节点的邻域合并形成的区域内，并相应地激活在同一区域内的一些节点来实现的。

　　节点的激活和停用之后是在两个域中识别适当的虚拟节点。图6底部显示了两者在时间瞬间的更新配置。Mossaiby等人(2022)提供了有关自适应算法实现的进一步细节。

　　重要的是要强调，即使切换策略限制了标准方案和精细网格间距仅适用于域的局部区域，自适应算法的实现对于涉及不连续的快速和不确定演化的问题，例如在动态裂纹扩展的情况下，比关键区域的演化较慢且更可预测的问题更方便。比如在生物降解过程中。这是由于算法的触发给模型增加了一定的开销。鉴于此，在第4节给出的数值示例中，切换策略仅在动态脆性断裂案例研究中实现，而在腐蚀问题中，和域是先验定义的，并且在模拟过程中保持不变。除此之外，在腐蚀问题中，引入与相边界对应的更细的网格，无论是固定的还是自适应的配置，对于准确模拟腐蚀层附近的体积损失是必要的，也是至关重要的。

　　在本节中，通过两个关于动态裂纹扩展和腐蚀的实际例子来检验所提出的多自适应方案的性能。在第一个例子Kalthoff-Winkler实验中，将关键区域内网格间距等于，其余区域网格间距等于的(adaptive, hybrid)方案的解与其他三种方法的解进行比较:(均匀，粗)，将标准方案应用于网格间距等于的均匀离散化;(均匀，细)，将标准方案应用于网格间距等于的均匀网格;(自适应，标准)，基于Shojaei等人(2018)引入的多网格方法，并将标准方案应用于类似于所提出方法的非均匀离散化。

　　在验证了(自适应，混合)方案与其他方法的能力之后，在第二个例子中，即三维腐蚀问题，所提出的方法的解决方案仅与(均匀，精细)模型的解决方案进行比较。

　　(均匀，精细)方法的解决方案实际上在这里被视为参考，因为由于用于离散化的网格的高度细化，它是最准确地再现问题物理的模型。

　　在这项工作中，选择了较高的m值(特别是在第一个例子中)来证明所提出的多自适应方法在处理所研究的材料或现象的特征长度尺度较大的问题时的能力和有效性，即视界值较大，所需的分辨率也很高。

　　所有的方法都是通过使用c++编程语言开发的内部代码实现的。受Mossaiby等人(202,2017)的研究和Mossaiby等人在Mossaiby(2017,2022)中的PD代码的启发，开发的代码被大规模并行化。所有的模拟都是在一个内部高性能计算(HPC)集群的计算节点上进行的，该集群包括两个带有24核2.1 GHz英特尔至强可扩展白金8160处理器的插槽。通过使用OpenMP指令，可以在多达48个内核上实现共享内存多处理并行性。通过将I/O时间从总时间中排除，可以比较不同解决方案的运行时间。为了研究并行性对所提出的方法实现的速度的影响，在第一个示例中，报告了关于12核和48核并行性的两种场景的运行时间。在第二个示例中，只报告48核并行情况下的运行时间。

　　第一个例子中考虑的问题是Kalthoff - winkler实验(Kalthoff 2000)，这是一个众所周知的问题，也是PD框架中关于脆性材料动态断裂建模的众多研究的主题(Shojaei et al. 2018;Zaccariotto et al. 2018;Dipasquale et al. 2014;Ren et al. 2017,2019;Islam and Shaw 2020;Ning et al. 2022)。图7(左)图解地表示了该问题的几何结构，该问题由一块矩形钢板组成，钢板上有两个平行的预裂缝，钢板受到来自刚性弹丸的冲击载荷，该弹丸在预裂缝之间的边界区域横向撞击钢板。正如Kalthoff(2000)所全面描述的那样，实验观察到断裂类型取决于冲击速度。如果板由18Ni1900钢制成，如本案例所示，mm/s的冲击速度导致脆性断裂主要以模式i发生。事实上，一旦发生冲击，裂纹开始从预裂纹尖端扩展，几乎以与水平方向近似的角度直线扩展，直到接近板边界。

　　图7

　　

　　例1中问题域的示意图(左)和自适应方法的初始配置(右)

　　图8

　　

　　例1中采用不同离散化方案的PD邻域表示。在标准方案中，精细网格和粗网格的邻域分别运行1793和441个节点。然而，在混合离散化中，计算只运行在粗网格内邻域的49个配置节点上

　　主要问题参数的取值与Dipasquale等人(2014)报道的相似，即考虑了平面应变条件，分别为N/(杨氏模量)、Kg/(质量密度)、J/(断裂能)和泊松比。冲击荷载是通过在沿两个预裂缝之间的板左侧部分的节点上施加mm/s的初始速度来模拟的。施加的速度沿水平方向定向，并在整个模拟过程中保持恒定。通过采用本节之前介绍的四种不同方法并假设mm和mm(即，，，其中为域粗部分使用的比率，而为域细部分使用的比率)来研究该问题。时间积分是通过利用velocity-Verlet算法进行的，如第3.1节所述。总时间长度为87.5，时间步长为。对于所提出的多自适应方案，即(自适应，混合)方法，在使用搭配方案离散化的区域部分，用于执行计算的单项基函数被取到二阶(参见(21))，这在二维中导致单项和搭配节点的数量分别为和。图8显示了考虑实例1求解所采用的不同离散化方案所得到的PD邻域。图中清晰地展示了使用本文方案相对于标准方案的优势，因为在考虑标准方案时，精细网格和粗网格内的PD邻域分别包含1793和441个节点，而在本文方法中，粗网格内的PD邻域仅包含49个搭配节点。(自适应，混合)和(自适应，标准)方法在仿真开始时的配置示意图如图7(右)所示。黑色用于表示域，其中包括预裂纹和边界。如前所述，在(自适应，混合)方法中，和分别使用标准方案和搭配方案进行离散化，而在(自适应，标准)方法中，两个领域都通过标准方案进行离散化。对于自适应方法中要利用的自适应算法，将因子(cf.(42))设为等于0.95。

　　图9显示了研究中考虑的四种模型得到的裂纹模式，以及两种自适应方法在四个不同时刻的域演化配置。除了(均匀，粗糙)模型外，所有模型都得到了相似的损伤轮廓图，因为在这种情况下，裂纹扩展角度和速度与其他方法模拟的结果略有不同。此外，在这种情况下，继发性骨折起源于与受冲击板相反的一侧。相反，其他模型没有报告这种现象。因此，(自适应，混合)方法的解与(均匀，精细)和(自适应，标准)方法的解非常吻合，从而证实了所提出的多自适应方案在研究动态裂纹扩展问题中的适用性。为了进一步验证所提方案，图9中考虑的四种不同模型在同一时刻得到的垂直方向速度场等值线图如图10所示。并通过实例验证了(自适应，混合)方案与(均匀，精细)和(自适应，标准)方案的结果吻合良好。由(均匀，粗糙)模型得到的解与其他方法的差异也反映在速度场中。

　　图9

　　

　　例1中不同模型得到的裂纹图;还包括(自适应的，标准的)和(自适应的，混合的)方法领域的发展

　　图10

　　

　　例1中采用(均匀，粗)、(均匀，细)、(自适应，标准)、(自适应，混合)四种方法获得的四个时刻速度场垂直方向等值线图;单位为[mm/s]

　　将(自适应，混合)模型的计算成本与(均匀，精细)和(自适应，标准)方法的计算成本进行比较，突出了(自适应，混合)模型的数值效率。表1报告了模拟开始和结束时每个域中(即自适应方法)的活动节点数量以及上述三种模型的运行时间。为了更好地观察所提出的多自适应方案的性能，还给出了加速。报告了两种不同并行配置(即12核和48核)的运行时间和相对加速，以便研究并行性对自适应方法在(统一的、精细的)模型上实现的加速的影响程度。如表1所示，当在12核上运行模拟时，两种自适应方案都比(均匀，精细)方法快1.86倍，而当模拟在48核上运行时，(自适应，混合)方案是计算效率最高的方法，因为它比(均匀，精细)模型快1.19倍，比(自适应，标准)方案快1.09倍。甚至

　　表1例1的运行时间。(毫米)

　　图11

　　

　　对比例1中(均匀、粗)、(均匀、细)、(自适应、标准)、(自适应、混合)四种方法在以预裂顶端为中心半径为10mm的圆形区域内体系的能量含量

　　虽然在使用12核进行模拟时，与使用两种自适应方法相关的计算资源优化更为明显，但与基于标准方案的自适应策略相比，使用更多核凸显了所提出方案的更高效率。考虑到表1中报告的结果并回顾图9和图10，可以得出结论，多自适应方案能够以显着较低的计算成本准确地再现参考解，即(均匀，精细)方法的解。

　　为了进一步验证多自适应方法的性能，图11显示了研究中考虑的四种不同模型在整个仿真过程中系统能量含量的变化。对于每种方法，能量含量都是在位于上部预裂纹尖端的圆形区域内计算的。从图11中可以看出，(自适应，混合)模型与(均匀，精细)和(自适应，标准)模型的结果非常吻合，表明本文方案能够以与上述模型相同的精度模拟裂纹形核和裂纹扩展速度，但在计算效率上优于上述模型，如表1所示。

　　在这项工作中提出的第二个例子涉及镁骨植入螺钉的生物降解。在过去的几十年里，这个问题引起了越来越多的关注，因为镁植入物结合了骨样机械和骨传导特性以及良好的生物相容性(Negrescu et al. 2020)。此外，由于它们在生物体中由于生物腐蚀而自然降解，它们可以用作临时的机械支撑，在周围组织愈合后，不需要在第二次手术干预中移除，但随着时间的推移，它们会自然溶解。不幸的是，迄今为止，由于其应用的高度复杂性，化学、机械和生物因素相互作用，以及与相关实验研究相关的相当大的成本，这种植入物的发展进展仍然受到很大的阻碍。因此，数值模拟越来越被视为克服这些限制的最有前途的工具(Gartzke et al. 2020;Ma et al. 2018)。尽管基于经典局部方法的成功发展，PD扩散求解器更适合于模拟生物降解过程，因为在这些腐蚀问题中，腐蚀前沿对应的浓度场实际上受到高度不连续的影响，从而使经典局部扩散模型难以充分处理这种现象。相反，PD模型可以方便地捕捉相边界(例如腐蚀界面)的演变，作为其解决方案的自然结果(Hermann et al. 2022)。尽管本文提出的方法的应用仅限于均匀腐蚀的建模(参见Jafarzadeh等人(2019b)了解均匀腐蚀建模的另一种方法)，但在过去十年中，PD已被广泛用于模拟不同类型的腐蚀现象，如点蚀(Chen和Bobaru 2015a;Jafarzadeh et al. 2019;Wang等人，2023)，带带状覆盖的点蚀(Jafarzadeh等人，2019c)，晶间腐蚀(Jafarzadeh等人，2018)，电蚀(Zhao等人，2021)，缝隙腐蚀(Jafarzadeh等人，2022)，应力辅助腐蚀(Jafarzadeh等人，2019b;Fan et al. 2022)和应力腐蚀开裂(Chen et al. 2021)。可以预见，该方法将扩展到其他腐蚀过程和损伤问题的模拟中。

　　图12

　　

　　示例2中问题域和初始配置的示意图表示

　　本例中考虑的三维问题包括一个金属试样，即如图12所示的Mg-10Gd植入螺钉，浸入充满液体电解质的球体中。该植入螺钉的特点是半径为1mm，高度为4mm，螺钉头为mm开槽。这里取球体半径为6mm，圆柱形区域的半径和高度分别为2.1 mm和6.2 mm，而固相(即完整金属)和液相(即电解质)的初始浓度分别为和。由于电解液对金属的腐蚀过程，固相减少，溶解的物质通过所谓的溶解通量逐渐扩散到液相中。为了使PD扩散模型(见2.2节)适用于腐蚀过程的模拟，有必要引入两个微扩散系数参数，即和微溶解性，这两个参数分别对应于连接液相材料点对的键和连接固相材料点与液相材料点的界面键。在固相中，不考虑材料点对之间的通量。另一个用于腐蚀过程建模的输入参数由阈值浓度表示。该参数用于确定整个溶解过程中固相的浓度。更具体地说，一旦固相中某一物质点的对应浓度低于该阈值，该物质点的相就由固相变为液相。由于这种相变机制，腐蚀前沿可以在模拟过程中自主进化。

　　图13

　　

　　例2中采用不同离散化方案的PD邻域表示。标准方案的计算在精细网格内的邻域运行超过7153个节点。在混合方法中，计算只运行在粗网格内的123个邻域的并置节点上

　　主要问题参数的值为、、、和(参见(12))。dirichlet型边界条件通过假设在距离球体表面有限距离内的所有节点上的约束浓度等于来实现，这里认为等于地平线(参见第2.2节)。该问题的研究采用了前面介绍的两种方法，即(均匀，精细)和(自适应，混合)模型，并假设mm和mm(即，，，其中为域粗部分使用的比例，而为域细部分使用的比例)。为了与4.1节示例1中所进行的讨论一致，所提出的多自适应方案在这里也被称为(自适应，混合)，尽管在这种情况下，如前所述，在腐蚀过程的模拟过程中不触发自适应算法。在所提出的方法中，图12中表示的和域的配置实际上是先验定义的，并且在整个模拟过程中保持固定。通过标准方案进行离散化的区域是为了只覆盖浓度场不连续的区域，即不同相特征的节点之间发生扩散的对应区域(相边界)。在使用搭配方案离散化的域部分(即域)中，用于执行计算的单项基函数被取到二阶(参见(21))，这在三维中导致了以下单项和搭配节点数量的规范，即分别为和。图13为例2求解采用不同离散化方案后得到的PD邻域。该图清楚地展示了使用本文方法相对于(均匀，精细)方法的优势，因为在考虑标准方案时，精细网格内的PD邻域由7153个节点组成，而在本文方法中，粗糙网格内的PD邻域仅包括123个节点。时间积分是通过利用3.1节中提到的显式前向欧拉算法来完成的(参见(17))。通过考虑时间步长为的总时间持续时间为14，解决了这个问题。

　　Mg-10Gd植入螺钉在液体电解质中浸泡14天后的PD腐蚀模拟结果如图14所示。该图显示了研究中考虑的两种模型在两个不同时刻得到的液相内扩散的固体金属和溶解物质的浓度分布图。为了清晰起见，固相和液相分别表示，只显示了液相球的一半。报告的结果表明，(自适应，混合)方案得到的溶液与参考方案非常吻合，因为在固相中可识别的腐蚀模式和溶解物质在液相中的扩散几乎完全符合(均匀，精细)方法模拟的结果。因此，结果证实了所提出的多自适应方案对大规模几何复杂问题中腐蚀过程建模的适用性。为了进一步验证所提出的方法，图15显示了上述两种模型的固相相对体积损失随时间的变化结果。对于计算成本问题，(自适应，混合)方案的解决方案在腐蚀的前2天内通过(均匀，精细)方法的解决方案进行验证。可以推断，所提出的方法重新获得了参考腐蚀速率。尽管如此，(自适应，混合)方案的模拟持续了14天的生物降解。结果显示，经过14天的腐蚀，超过一半的固相溶解了。

　　图14

　　

　　例2中采用(均匀，精细)和(自适应，混合)两种方法获得的两个时刻浓度场等值线图;单位为[mol/mm]。在(自适应，混合)方案中，相对于(均匀，精细)解决方案中相同区域的颜色，(自适应，混合)方案中的区域颜色较浅，这只是一种视觉效果，与用于离散该部分区域的较粗网格间距有关

　　通过与(均匀、精细)方法的计算成本比较，突出了(自适应、混合)方法的数值效率。表2报告了上述两种模型的每个域中活动节点的数量(即自适应方法的数量)和平均运行时间(每100个时间步)。为了更好地了解所提出的多自适应方案的性能，还给出了平均加速(每100个时间步长)。报告了48核并行情况下的运行时间和相对速度提升。如表2所示，(自适应的，混合的)方案比(统一的，精细的)方法快11.76倍。考虑到这些结果并回顾图14和图15，可以得出结论，所提出的多自适应方案能够以相当低的计算成本精确地再现(均匀，精细)方法的解。

　　图15

　　

　　在例2中，通过两种方法(均匀，精细)和(自适应，混合)获得固相随时间的相对体积损失

　　表2例2的运行时间。(毫米)

　　表1和表2的结果表明，在三维腐蚀实例中，(自适应、混合)方案获得的速度远远大于二维动态裂纹扩展实例。这与以下事实有关:在腐蚀问题中，与触发自适应算法相关的开销不存在，因为和域是先验定义的，并且在整个模拟过程中保持不变。

　　本文提出了一种多自适应方法来提高PD模型的数值性能，从而优化计算资源的使用。在提出的框架中，标准方案(通常用于实现PD模型的基于粒子的方法)的应用仅限于受不连续存在或具有复杂几何细节(即关键区域)和边界影响的局部区域，而最近提出的搭配方案则用于该领域的其余部分。

　　实施这种混合离散化的原因是，尽管配置方案在计算成本和空间积分计算方面具有优势，但该方案在靠近关键区域和靠近物体边界时可能会受到精度损失的影响，其中，由于能够方便地在非线性PD模型所描述的区域中计算所需的积分项，从而处理物质分离，因此已知标准方案的性能更好。

　　为了最大限度地提高PD模型的数值效率，本文将这种混合离散化方案与多网格方法结合在一起实现，其中不同间距的网格在同一域中耦合，从而将细网格间距的应用限制在所需域的一小部分(例如，关键区域和边界)。因此，在提出的多自适应框架内结合了混合离散化方案和多网格方法的优点，并通过两个实际数值实例评估了其处理二维和三维动态脆性断裂和腐蚀问题的能力。在Kalthoff-Winkler实验建模和镁骨植入螺钉生物降解的实例中，证明了在给定的精度要求下，所提出的多自适应方案能够再现PD模型得到的解，使用标准方案进行离散化，网格间距非常细，计算成本显著降低。结果还表明，在涉及数百万自由度的大规模、几何复杂的三维问题中，所开发方案的效率更为明显。

　　ccDownload: /内容/ pdf / 10.1007 / s10704 - 023 - 00709 - 8. - pdf