一种无网格固体力学仿真方法、电子设备及存储介质与流程

本发明涉及力学仿真系统，尤其涉及一种无网格固体力学仿真方法、电子设备及存储介质。

背景技术：

结构应力分析在航空、航天、汽车等不同的工程领域都是十分重要和必要的。从产品的设计、制造到维护，结构应力分析起着至关重要的作用。由于其重要性，几十年来，许多研究人员一直致力于提高这一过程的准确性和效率。

目前，有限元法(finiteelementmethod)在结构应力分析中得到广泛的使用。该方法采用了基于单元的、局部的、多项式的、连续的试验和检验函数来进行仿真。由于试验和检验函数都是多项式的，因此，有限元法中的积分是很容易计算的，并且总体刚度矩阵的对称性和稀疏性使得有限元法适用于大型结构仿真。然而，有限元法的精度在很大程度上取决于网格的质量；而为了得到满意的解决方案，通常会花费很多精力在网格划分上。而且，即使在初始计算时采用高质量的网格结构，在结构的变形较大的情况下，有限元法会出现网格畸变，使得求解精度大大降低。另外，为了使有限元法能够研究裂纹的形成、结构的断裂和破碎，必须反复重新划分网格、删除单元或使用内聚力模型，但是这些额外的手段使得有限元法在计算极端问题时费时费力，结果也不够准确。

无网格方法是上世纪末发展起来的一类部分或完全消除网格的方法。无网格方法主要分为无网格弱形式方法和无网格粒子法两种。无网格伽辽金(elementfreegalerkin)和无网格局部彼得洛夫伽辽金方法(meshlesslocalpetrov-galerkin)是两种经典无网格弱形式方法。这两种无网格方法主要利用移动最小二乘(movingleastsquares)或径向基函数(radialbasisfunction)近似来推导基于节点的试验函数。使用这两种近似方法，很容易实现函数的高阶连续性。此外，由于独立的节点取代了网格结构，无网格方法可以方便地插入或删除额外节点。然而，另一方面，移动最小二乘法或径向基函数法得到的试验函数不再是多项式函数，并且非常复杂。因此，无网格伽辽金和无网格局部彼得洛夫伽辽金方法中的积分计算极其繁琐，精度较低，甚至会影响整个方法的精度和稳定性。当然，为了降低计算量，提高积分精度，必须采用一些特殊的、新型的数值积分方法，如一系列的新型节点积分法。

光滑粒子流体动力学方法(smoothedparticlehydrodynamics)是一种形式简单、计算量小的无网格粒子方法。然而，光滑粒子流体动力学方法是一种基于强形式，而不是弱形式的方法。这种方法的稳定性较差，容易发生拉伸失稳等问题。再比如，近场动力学(peridynamics)方法作为一种相对较新的方法，其理论基础不同于传统的连续介质力学。因此，科学家和工程师们几十年来积累的经典的本构模型和工程经验很难准确、直接地应用在近场动力学方法中。

技术实现要素：

为了克服现有技术的不足，本发明的目的之一在于提供一种无网格固体力学仿真方法，其能够提高计算效率与仿真精度，并有利于解决传统仿真方法难以进行断裂、破碎等极端问题仿真的难题。

本发明的目的之二在于提供一种电子设备，其能够提高计算效率与仿真精度，并有利于解决传统仿真方法难以进行断裂、破碎等极端问题仿真的难题。

本发明的目的之三在于提供一种计算机可读存储介质，其能够提高计算效率与仿真精度，并有利于解决传统仿真方法难以进行断裂、破碎等极端问题仿真的难题。

本发明的目的之一采用如下技术方案实现：

一种无网格固体力学仿真方法，包括以下步骤：

试验函数和检验函数的构造步骤：根据物体的结构中分布的节点，得出每个子域中由节点位移表示的试验函数和检验函数的表达式；

伽辽金弱形式的构造步骤：引入数值通量修正得到带有数值通量修正的伽辽金弱形式；

代数方程组生成和求解步骤：根据每个子域的试验函数和检验函数的表达式，以及带有数值通量修正值的伽辽金弱形式，构造出结构的总体刚度矩阵和代数方程组，并求解代数方程组得出每个节点的位移值；

结构变形和应力计算步骤：根据每个节点的位移值，求出该结构各个位置的变形以及应力，完成结构的仿真分析。

进一步地，还包括模拟步骤：根据一定时间顺序的多个载荷，依次在每个载荷下对结构执行试验函数和检验函数的构造步骤、伽辽金弱形式的构造步骤、代数方程组生成和求解步骤以及结构变形和应力计算步骤得出结构的各个位置的变形以及应力，并采用一定的判据来判断结构中相邻子域之间是否开裂，模拟结构在一定时间顺序的多个载荷作用下的断裂与破碎过程。

进一步地，每个子域内均包含一个节点。

进一步地，所述试验函数和检验函数的构造过程如下：

步骤s11：首先定义每个子域的试验函数的类型，并记为u^h，该试验函数的类型为n阶多项式函数；其中，n为大于或等于1的自然数，一般取n等于1；如公式(1)所示，将试验函数在节点p0处进行n阶泰勒展开：

其中，是位移u^h在节点p0处的值，是节点p0处的1到n阶位移梯度；

步骤s12：确定每个子域的支持域，并且定义加权离散的l²范数j，其中，j＝(aa+u0-um)^tw(aa+u0-um)(2)；

(xi,yi)是节点pi的坐标，是位移u^h在节点pi处的值，是权重在节点pi处的值(i＝1,2,3,…m)；

步骤s13：根据公式(2)以及范数j的驻值条件，将位移梯度a由节点位移值表示出来，也即是：

a＝cue(3)；

其中，c＝(a^twa)-¹a^tw[i1i2]，

步骤s14：根据公式(3)和公式(1)得出每个子域中试验函数的表达式；

步骤s15：根据检验函数和试验函数的表达式具有相同的形式得到检验函数的表达式。

进一步地，所述支持域的定义包括但不限于以下任意一种：定义每个节点的支持域为节点的所有邻近节点所在的子域、定义每个节点的支持域为以节点为圆心的特定大小的圆内部的所有节点所在的子域；并且，当相邻子域之间存在裂纹时，裂纹一侧的节点不在另一侧节点的支持域中。

进一步地，所述数值通量修正的类型包含但不仅限于内罚函数数值通量修正。

进一步地，所述代数方程组生成和求解步骤还包括：

步骤s21：将每个子域的试验函数和检验函数带入数值通量修正后的伽辽金弱形式得出节点刚度矩阵和内部边界刚度矩阵，节点刚度矩阵ke＝∫eb^tdbdω，内部边界刚度矩阵e是子域的公共边；

步骤s22：根据节点刚度矩阵和内部边界刚度矩阵装配成总体刚度矩阵k，从而得到矩阵形式下的代数方程组：kq＝q，其中q是节点位移向量，q是载荷向量；

步骤s23：根据总体刚度矩阵和载荷向量，求解代数方程组得出结构的节点位移向量，进而得出结构中每个位置的位移值。

本发明的目的之二采用如下技术方案实现：

一种电子设备，包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现如本发明目的之一采用的一种无网格固体力学仿真方法的步骤。

本发明的目的之三采用如下技术方案实现：

一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现如本发明目的之一采用的一种无网格固体力学仿真方法的步骤。

相比现有技术，本发明的有益效果在于：

本发明通过提供一种新型不连续的试验函数、检验函数及相应的弱形式，形成了一种新型无网格固体力学仿真方法，其能够提高计算效率与仿真精度，并有利于解决传统仿真方法难以进行断裂、破碎等极端问题仿真的难题。

附图说明

图1为本发明提供的结构与随机分布的节点示意图；

图2为本发明提供的结构与其子域划分示意图；

图3为本发明提供的节点p0的支持域示意图；

图4为本发明提供的不连续的形函数示意图；

图5为本发明提供的不连续的试验函数示意图；

图6为本发明提供的无网格固体力学仿真方法的流程图；

图7为本发明提供的另一种结构与随机分布的节点示意图；

图8为图7中的结构在x方向上的位移场示意图；

图9为图7中的结构在x方向上的应力场示意图；

图10为本发明提供的含有裂纹平板结构与均匀分布节点示意图；

图11为本发明提供的裂纹两侧节点p1，p5的支持域与自由边界示意图；

图12为图10中含有裂纹平板结构在x方向上的位移场示意图；

图13为图10中含有裂纹平板结构在y方向上的应力场示意图。

具体实施方式

下面，结合附图以及具体实施方式，对本发明做进一步描述，需要说明的是，在不相冲突的前提下，以下描述的各实施例之间或各技术特征之间可以任意组合形成新的实施例。

一种无网格固体力学仿真方法，具体包括以下三个部分：

第一部分：构造试验函数和检验函数。

对于一个如图1所示的给定形状的固体结构来说，首先将若干个节点分布在结构的内部和边界上。如图2所示，利用这些节点，将结构划分为多个任意形状的子域，并且在每个子域内都包含一个节点。对于子域的划分，本发明中不做限定，只需要保证在每个子域内都包含一个节点即可，比如泰森多边形方法。

构造试验函数和检验函数就是基于这些分布的节点及其子域进行构造的，具体步骤如下：

步骤a11：首先为每个子域定义一个多项式的试验函数(也即是位移场，是指结构变形前到变形后的位移)u^h。该试验函数可以为一、二、三阶或更高阶的多项式函数。

比如公式(1)所述，将试验函数在节点p0处进行n阶泰勒展开：

其中，是位移u^h在节点p0处的值，是节点p0处的1到n阶位移梯度。

当n＝1，也即是试验函数为一阶函数时，p0所在子域e0的试验函数u^h可表示为：

其中，u^h(x,y)表示坐标点(x,y)处u^h值，(x0,y0)是点p0的坐标，是u^h在节点p0处的值，节点p0处的位移梯度

为了将试验函数用节点位移表示出来，需要将位移梯度用节点位移表示出来。其中，节点位移是指节点发生的位移，比如图2和图3中黑点处的位移。

为了将位移梯度用节点位移来表示，首先需要对节点的支持域进行设定。本发明中对于支持域的设定可采用多种定义，例如：定义每个节点的支持域为该节点的所有邻近节点所在的子域，如图3所示。

当然，也可以采用如下定义：定义每个节点的支持域为以每个节点为圆心的特定大小的圆内部的所有节点所在的子域。另外，如果相邻子域之间存在裂纹，裂纹一侧的节点不在另一侧节点的支持域中。

步骤a12：确定每个子域的支持域，并且定义加权离散的l²范数j，其中，j＝(aa+u0-um)^tw(aa+u0-um)(2)，

(xi,yi)是节点pi的坐标，是位移u^h在节点pi处的值，是权重在节点pi处的值(i＝1,2,3,…m)；

步骤a13：根据公式(2)以及范数j的驻值条件，将位移梯度a由节点位移值表示出来，也即是：

a＝cue(3)；

其中，c＝(a^twa)^-1a^tw[i1i2]，

步骤a14：根据公式(3)和公式(1)得出每个子域中试验函数的表达式u^h＝nue(4)。

其中，这里的n称为形函数。

规定检验函数和试验函数的表达式具有同样的形式，也即是试验函数和检验函数具有相同的形函数。因此，采用同样的方法可以得到检验函数v的表达式：v＝nve(5)。

因此，根据位移场和应变、应力场的关系，可以将应变场ε^h和应力场σ^h分别通过试验函数和检验函数表示出来，具体如下：

其中d为材料刚度矩阵，与结构的材料属性有关，可以为各向同性、正交各向异性或其他属性。

第二部分：引入数值通量修正伽辽金弱形式。

根据结构中每个子域的试验函数和检验函数的表达式，就可以得到整个结构的试验函数和检验函数的表达式。

每个子域的试验函数、检验函数在子域内部是连续的。而对于整个结构来说，由于每个子域在公共边界上有不同的位移场，因此，试验函数和检验函数在整个结构中是不连续的。

例如，图5是试验函数的示意图，可以看出整个结构的试验函数是不连续的；图4为一个节点的形函数的示意图，同样的，形函数也是不连续的。

本发明是应用于现有的无网格伽辽金的方法中的，而现有的伽辽金弱形式的公式如公式(6)：

但是，由于本发明给出的相邻子域的试验函数和检验函数是不连续的，若直接将本发明计算出每个子域的试验函数和检验函数运用到伽辽金弱形式的公式中时，将导致数学上的不一致性，最终引起计算结果的不稳定性。

为了解决由于试验函数和形函数的不连续性导致的计算结果的不稳定，本发明引入数值通量修正来对现有的伽辽金弱形式进行改进。

本发明中的数值通量修正可以有多种类型。由于数值通量修正限制了函数在内部边界上的不连续性，因此可通过引入数值通量修正来解决的试验函数不连续造成的计算结果不稳定的问题。

本实施例以内罚函数数值通量修正为例，来对现有的伽辽金弱形式进行改进。公式(7)表示带有内罚函数数值通量修正的伽辽金弱形式。该公式(7)中力边界条件和内部边界的位移连续条件以弱形式施加，位移边界条件需要以强形式施加。

其中γu为位移边界，γt为力边界；{}为平均算子，[]为跳跃算子；当e为内部边界(子域e1和e2的公共边)时，是从e1指向e2的单位法向向量，并且

从上可知，本发明中所提供的带有内罚函数数值通量修正的伽辽金弱形式的公式(7)与传统的伽辽金弱形式的公式(6)相比，其区别在于前者在公共边界上存在跳跃和平均项，进而可解决试验函数和检验函数的不连续性导致的计算结果的不一致性。

第三部分：将每个子域的试验函数和检验函数的表达式代入到带有数值通量修正的伽辽金弱形式中，构造出结构的总体刚度矩阵和代数方程组，并根据代数方程组求解出结构中每个节点的位移值，进而根据结构中每个节点的位移值计算得出该结构中每个位置的变形和应力，完成结构的仿真计算。

步骤b1：首先将每个子域中的试验函数和检验函数代入到带有数值通量修正的伽辽金弱形式中。也即是将每个子域的试验函数和检验函数代入到公式(7)，得出每个节点的节点刚度矩阵和内部边界刚度矩阵。

其中，节点刚度矩阵ke＝∫eb^tdbdω，内部边界刚度矩阵e是子域的公共边。

步骤b2：:通过每个节点的节点刚度矩阵和内部边界刚度矩阵装配得出结构的总体刚度矩阵k，得出矩阵形式下的代数方程组：kq＝q(8)。

其中，q是待求解的节点位移向量，q是载荷向量，总体刚度矩阵k由节点刚度矩阵和内部边界刚度矩阵装配得到。

步骤b3：根据结构的总体刚度矩阵和当前的载荷向量，求解公式(8)的代数方程组，得出结构中每个节点的位移向量，也即是每个节点的位移值。

步骤b4：根据每个节点的位移值得出结构的各个位置的位移，并且通过位移与应力的关系得出结构中每个位置的应变以及应力。

通过本发明可以得出结构的变形和受力的分析结果，将其应用于结构设计过程中，用来判断机械、土木等物体结构的设计是否安全合理。

其中，载荷是指结构在工作状态下所受的力，可以是单一载荷工况，可以是一个载荷历史。同样的结构在不同的载荷作用下会产生不同的变形与应力。本方法如果应用于结构在工作状态下受到的一个载荷历史，可以将该载荷历史分为m步，在每一个载荷步中求解变形与位移，以及模拟裂纹的产生与扩展。

本发明还结合具体的图示给出了对于结构的应变以及应力的计算过程。

实施方式一：

针对如图7所示的结构：

步骤c11：首先利用结构中随机分布的节点，将结构划分为若干子域。

步骤c12：确定每个节点的支持域，使用加权最小二乘法计算该结构中每个子域内试验函数和检验函数的表达式。

其中，试验函数：u^h＝nue，n为形函数，ue为节点位移。

检验函数：v＝nve，和试验函数具有相同的形函数n。

步骤c13：引入数值通量修正得到带有数值通量修正的伽辽金弱形式。

步骤c14：将每个子域的试验函数和检验函数代入含有数值通量修正的伽辽金弱形式。也即是公式(7)，然后使用高斯积分法进行数值积分得到节点刚度矩阵ke，内部边界刚度矩阵kh。由于试验函数和检验函数均为多项式函数，因此采用高斯积分法进行数值积分时能够快速、准确地计算得出节点刚度矩阵ke，内部边界刚度矩阵kh。

步骤c15：将所有的节点刚度矩阵ke，内部边界刚度矩阵kh装配出结构的总体刚度矩阵k，并根据给定的载荷向量得出代数方程组为：kq＝q。q为给定的载荷向量，q为节点位移向量。

通过强制施加位移边界条件，计算出节点位移向量q，并将其带回试验函数的表达式中，得到结构的位移场。

根据位移场和应力场的关系，得到结构的应力场。

例如：图7中的结构在x方向的位移场如图8所示、x方向的应力场如图9所示。上述仿真结果与理论情况下的精确值作对比，位移场的相对误差为0.051％，应力场的相对误差为2.064％。

实施方式二：

针对如图10所示的含裂纹的平板结构：

步骤c21：依据同样的划分子域的方法将该结构划分为若干子域。

步骤c22：确定每个节点的支持域。由于该结构中存在裂纹，因此，在确定支持域时，定义将裂纹一侧的节点不会出现在另一侧节点的支持域中。如图11所示，节点p1的支持域中不含p5，节点p5的支持域中不含p1。

同时，将与裂纹重合的边界设为自由边界。如图11所示，边界γ26，γ15和γ47为自由边界。

步骤c23：使用加权最小二乘法计算每个子域内试验函数的表达式u^h＝nue，n为形函数，ue为节点位移。

检验函数v＝nve，和试验函数具有相同的形函数。

步骤c24：将每个子域的试验函数和检验函数代入含有数值通量修正的伽辽金弱形式，即公式(7)之中，使用高斯积分法进行数值积分得到节点刚度矩阵ke，内部边界刚度矩阵kh。

将所有节点刚度矩阵ke，内部边界刚度矩阵kh装配出总体刚度矩阵k，并当前的载荷向量得出代数方程组为：kq＝q。q为当前的载荷向量，q为节点位移向量。

通过强制施加位移边界条件，计算出节点位移向量q，并将其带回试验函数的表达式中，得到结构的位移场。

根据位移场和应力场的关系，得到结构的应力场。

图10的含裂纹的平板结构中x方向的位移场如图12所示、y方向的应力场如图13所示。本发明使用j积分计算出的应力强度因子与精确值误差为0.83％。

另外，结构在工作状态下，不仅有可能受到静态载荷，也可能经历一个载荷历史，这里的载荷历史包括按照时间顺序的一系列载荷。这时可以将该载荷历史分为m步，在每一个载荷步中使用该方法。然后根据计算得出的结构的各个位置的变形以及应力，采用一定的判据来判断结构中的相邻子域之间是否开裂。这样，就可以模拟该结构在一定载荷历史作用下的断裂与破碎过程。

其中，对于静态载荷，是指结构在工作状态下只受到一个载荷时，通过本发明所提供的一种无网格固体力学仿真方法得出结构中各个位置的变形以及应力，实现结构的仿真。

对于载荷历史，是指结构在工作状态下按照一定时间顺序受到的一系列载荷。根据一定时间顺序的多个载荷，依次对每个载荷执行本发明提供的一种无网格固体力学仿真方法的步骤得出结构中各个位置的变形以及应力，采用一定的判据来判断结构在对应载荷下相邻子域之间是否开裂，这样就可以模拟该结构在该载荷历史作用下的断裂与破碎过程，如图6所示。

本发明相对于现有的技术来说，其具有以下技术效果：

1、本发明提供了一种多项式的、不连续的试验函数和检验函数。

本发明通过一种全新的基于点的插值方式将结构划分为若干子域，然后计算每个子域的试验函数和检验函数，进而可得出的试验函数、检验函数在每个子域内均为多项式函数；而由于试验函数、检验函数在相邻子域的公共边界上存在不同的值，因此，在不同子域之间，每个子域的试验函数和检验函数是不连续的。

2、本发明还通过对现有的伽辽金弱形式进行了改进，通过引入数值通量修正来解决由于本发明提供的试验函数和检验函数的不连续性，而导致计算结果不一致的问题。

3、在计算结构的刚度矩阵时，由于本发明中试验函数和检验函数的多项式性质，采用高斯积分法或解析法可快速、简单、精确地进行积分计算。

4、由于本发明中总体刚度矩阵的对称、稀疏、正定性，本发明可适用于大规模结构仿真。

5、由于试验函数和检验函数的不连续性，本发明可以十分简单地切断节点之间的联系，并可以结合一定的开裂判据，简单高效地模拟裂纹的生成，扩展以及结构的断裂和破碎。

本发明还提供了一种电子设备，其包括存储器、处理器以及存储在存储器上并可在处理上运行的计算机程序，所述处理器执行所述程序时实现如文中所述的一种无网格物理力学仿真方法的步骤。

本发明还提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现如文中所述的一种无网格物理力学仿真方法的步骤。

上述实施方式仅为本发明的优选实施方式，不能以此来限定本发明保护的范围，本领域的技术人员在本发明的基础上所做的任何非实质性的变化及替换均属于本发明所要求保护的范围。