在固定的独立覆盖网格中求解几何非线性问题

苏海东，董 鹏，颉志强

(1.长江科学院 材料与结构研究所，武汉 430010；

2.水利部水工程安全与病害防治工程技术研究中心，武汉 430010)

在计算力学中，当材料发生的位移或变形相对于材料自身的几何尺度不是很小时，就必须考虑由大位移或大变形造成的非线性，称为几何非线性问题。目前，相关的研究工作大都采用拉格朗日描述法或欧拉描述法[1](以下分别简称拉格朗日法和欧拉法)来求解，前者主要用于固体力学，后者主要用于流体力学。

拉格朗日法(又称为拉格朗日观点)，跟踪材料体上的物质点，计算网格附着于材料体而跟随物质点运动，其优点在于：可得到物质点的物理量时间历程，有利于求解随时间变化的材料非线性问题；
能准确描述材料体的边界随时间的运动；
控制方程相对简单。但其缺陷是特大变形引起网格扭曲，使计算精度急剧下降，甚至导致计算失败。常见的解决措施是重新划分网格，但有限元网格划分通常并不容易，而且需要通过插值运算在新、旧网格之间传递信息，易引入额外误差。

欧拉法(又称为欧拉观点)着眼于空间点，关注空间点上的物理量随时间的变化，材料体在固定于空间的网格中运动。其优点是计算网格始终保持不变，不存在特大变形下的网格扭曲，缺陷是：不易得到物质点的物理量时间历程；
难以准确描述材料边界随时间的运动；
控制方程存在迁移项(流体力学称为“对流项”)，造成处理上的困难，比如系数矩阵非对称，又比如基于欧拉观点的流体Navier-Stokes方程由对流项引起的非线性和计算稳定问题等。表1总结了拉格朗日法和欧拉法的特点，可以明显看出，2种方法优势互补。

表1 拉格朗日法和欧拉法的比较Table 1 Comparison between the Lagrangian method and the Eularian method

1964年Noh[2]提出任意拉格朗日-欧拉(Arbitrary Lagrange-Euler，ALE)方法，被认为吸收了拉格朗日法和欧拉法各自的优点：通过规定合适的网格运动形式(在ALE控制方程中考虑网格运动的影响)，既可以准确描述材料的移动边界，又能维持网格的合理形状。但要根据研究对象来设计网格运动方式需要经验和技巧，有时是非常困难的，会给使用者带来负担。

作为当今计算力学的研究热点，无网格法通过空间分布的离散点构造近似函数，从而避免了网格划分，其产生的一个重要背景就是针对有限元等网格类方法的网格扭曲问题[3]。无网格法中的物质点法(Material Point Method，MPM)[4]采用拉格朗日型的质点离散材料区域，用固定的欧拉背景网格求解动量方程，集合了2种描述的优势，避免了网格扭曲和迁移项。但无网格法的基础理论仍有待完善，计算量大的问题一直没有很好地解决。

笔者于2011年提出固定网格流形法[5]：利用数值流形方法的网格与材料分离的特性，采用在空间均匀布置的正三角形网格和一阶多项式覆盖函数，并考虑材料相对于空间固定网格的移动，初步实现了固定网格中的几何非线性分析。该方法同样集合了拉格朗日法和欧拉法的优势。然而，均匀网格不利于不同区域物理场的计算精度控制，而有限元网格通过结点相连的要求也给局部网格加密带来不便，同时，常规数值流形方法还存在高阶情况下的系统方程线性相关问题。

笔者后期提出独立覆盖流形法，具有网格任意连接和任意加密以及独立覆盖级数升阶的优点，且高阶情况的系统方程线性无关。本文在固定的独立覆盖网格中研究几何非线性问题，为下一步的计算精度控制和自适应求解以及固定网格中的多重非线性分析打下基础。为简便起见，本文仅研究单纯的几何非线性问题，不考虑由大应变引起的本构关系变化。

2.1 独立覆盖流形法简介

1991年，石根华先生[6-7]发明了数值流形方法(以下简称流形法)，将现代数学的流形思想引入数值计算：如图1(a)所示，一系列相互重叠的覆盖用于分割求解域，使复杂的物理场在每个覆盖简单化，以利于局部近似函数的逼近。

进一步提出基于独立覆盖的数值流形方法，简称独立覆盖流形法[10]，强调独立覆盖是计算分析的主要对象：如图1(c)所示，独立覆盖占据一个覆盖的主要面积，重点研究独立覆盖内的级数，而覆盖之间的重叠区域成为面积尽可能小(可只占覆盖面积的1%)的窄“条形”，仅起连接各覆盖级数的作用，并以有限单元的线性形函数作为单位分解函数。这样，由各覆盖分区的级数联合形成的整体近似函数，通过加权残值法(如伽辽金)逼近真实解，称为基于流形思想的“分区级数解”，可作为微分方程的最终解答。

图1 流形覆盖示意图Fig.1 Manifold covers

独立覆盖流形法的整体收敛是基于各覆盖自身的局部收敛而建立的[10]。以通常采用的多项式级数为例，相当于在各覆盖区域内将真实场函数展开为泰勒级数，或者说，各覆盖内的真实场函数(包括其导数，比如位移的导数是应变或应力)被多项式级数无限逼近。数学上的级数逼近理论保证了该方法具有严格的收敛性，且这种收敛性不受覆盖内实际占有的材料体(即真实的积分区域)形状的影响，也不受相邻覆盖在条形重叠区域是否错位连接的影响。

暂不考虑覆盖的条形重叠区域(将由程序自动添加，见下文)，新方法可应用2种覆盖网格方式：①对求解域进行任意形状的块体网格划分[11]，不同大小及形状的网格可随意连接；
②即本文所采用的规则网格方式，以图2(a)为例，用矩形网格(最好是正方形网格)覆盖求解域，通过切割操作将求解域分割成块体(即网格内的真实积分区域)。

基于独立覆盖流形法的收敛特性，材料不必占满所在网格的全部面积，真正有意义的是材料体所占据的积分区域，通常在材料体边界附近的网格内为任意形状。因此，网格不必与物理边界相匹配(本质边界条件可采用罚函数法施加)，这是数值流形方法及其发展而来的独立覆盖流形法的重要特性之一，突破了现有的大多数数值方法在网格划分上的限制，同时也是在固定网格中准确描述材料移动边界的基础。

当然，积分方法也与常见的映射到标准母单元的方式不同：以图2(a)右上角带有一条曲线边界的块体c为例，如图3(a)所示，将相邻顶点构成的各直线段依次与某一外点P0相连，形成有向单纯形(二维为三角形)[7]，在单纯形内可采用数值积分[11]或解析积分[12]，然后考虑各单纯形的正、负向后，累加得到整个块体上的积分；
图3(b)的曲线和直线所围区域，曲线边界可直接按精确的曲线方程参与积分计算[11,13]，实现材料边界的准确模拟；
然后如图2(b)所示，以块体e和f为例，在块体之间自动地添加条形，进行条形区域的积分，同时扣除重复的块体积分[11]。

大、小覆盖之间可任意错位相连，因而能对覆盖进行任意加密，比如图2(c)所示的对块体f进行细分和再细分。总之，覆盖网格具有任意形状(的积分区域)、任意(在条形处错位)连接以及任意加密这3个“任意”特性，配合计算精度控制，有望实现自适应分析甚至是无需人工参与的自动计算[14-15]。

图2 规则网格覆盖求解域Fig.2 Domain to be solved covered with regular meshes

图3 包含曲线边界的块体的单纯形划分Fig.3 A block with a curve boundary divided intoseveral simplexes

2.2 流形法计算格式及其改进

流形法采用分时步的计算格式，在t时刻计算完毕后更新材料体的边界，按更新后的构形计算t+Δt时刻的增量位移u，采用“惯性优势平衡方程式”[7]，即

(KL+Kg)u=Q+Fg-Fσ。

(1)

式中：KL为小变形情况下的线性刚度矩阵；
Q为荷载(本文考虑依赖于变形状态的所谓“跟随荷载”)；
当考虑动力荷载(不计阻尼)时，将惯性力分离出来，增加由此引起的刚度矩阵Kg及荷载Fg，即

(2)

(3)

Fσ称为初应力荷载,即

(4)

式中：ρ为密度(本文暂时忽略其变化)；
v为速度；
T为形函数矩阵；
BL为线性应变矩阵；
σ为线弹性应力;Ω为积分区域。每步计算完毕后的应力和速度按如下公式更新：

e=BLu,

(5)

t+Δtσ=tσ+De，

(6)

(7)

式中：相同变量以左上标区分其发生的构形；
e为线性应变增量；
D为线弹性材料矩阵。

与计算几何非线性问题常用的更新拉格朗日(Updated Lagrangian, UL)格式[1]相比，流形法计算格式的不同之处在于：

(1)缺少了几何刚度矩阵(也称为初应力矩阵)。在非线性分析中，刚度矩阵没必要很准确，只需保证正确的收敛方向。UL格式更接近于切向刚度，有利于收敛的稳定性，而流形法的刚度矩阵在初应力较大时有可能偏离方向而出现计算振荡，因此要引入惯性，使材料保持其原有的变形和运动方向。即使在静力分析中，也推荐采用动力松弛法[7]，即在式(7)中将速度tv乘上一个系数(如0.95)，由衰减的动力方式计算静力问题。

(2)没有非线性平衡迭代过程。在步长很小的情况下，不进行迭代计算而直接将失衡力计入下一步荷载是允许的。但需要控制步长，否则会引起计算结果的“漂移”。

(3)应变计算完全按照小变形公式，未考虑应变的二阶项，在计算带有旋转的大变形或大位移问题时，步长稍大就会产生所谓的“体积膨胀”。对此改进如下：式(6)中的线性应变增量e用Green-Lagrange应变增量ε代替，即

ε=e+η。

(8)

式中η为应变二阶项，其分量ηij用张量形式表示为[1]

(9)

式中：uk，i为k中的下标“i”表示位移分量uk对独立坐标xi求偏导数；
指标k重复出现2次，表示在该指标的取值范围内(如平面问题，k=1和2)遍历求和。

(4)其他如应力在不同构形中的转换、时间积分格式等差异还需进一步研究。

独立覆盖流形法沿用了经典流形法的数值计算格式，仅仅在数值解的空间离散上采用了“分区级数解”方案，表现为形函数T的不同形式：独立覆盖i中的位移u、速度v和应力σ等均为级数表达，比如多项式级数,即

(10)

式中：p为多项式的最高阶次；
ank为级数系数；
(x0，y0)取为独立覆盖的中心坐标；
lx和ly分别表示x和y这2个方向的尺度。而在2个覆盖之间的条形重叠区域内，有

V=φ1V1+φ2V2。

(11)

其中，单位分解函数采用一维有限元的形函数[1]，即

(12)

式中ξ为在条形厚度方向定义的局部坐标(以条形中点为原点)。

表1给出了拉格朗日法和欧拉法的优缺点。对于固体计算而言，拉格朗日法跟踪物质点、控制方程简单的优点很有吸引力，而引入欧拉法的固定网格以解决其网格扭曲问题的关键在于两点：①物质点相对于固定网格(或固定空间点)的移动；
②固定网格中的材料移动边界的准确模拟。本节主要讨论第①点。

如图4所示，在t时刻计算完毕后，根据位移级数解更新所有的材料点至t+Δt时刻的新位置，包括更新材料体的原边界(虚线)至新的构形边界(实线)。仍采用2.2节的拉格朗日型计算格式，考虑到积分区域可为任意形状的特性，在确定的形状下(即实线边界内)可以做到对式(1)至式(4)中的各矩阵和向量积分的准确计算，其中的速度v和应力σ的确定方式如下所述。

图4 逆向追踪示意图Fig.4 Schematic ofbackward tracking

拉格朗日方程关注的是物质点，关键是如何得到向量Fg和Fσ中关于物质点的速度v和应力σ。从有限元法来看，网格附着在材料体上并随材料移动，移动后的网格结点携带着物理量，通过插值方式可获得网格内任意点的物理量(不考虑网格扭曲造成的精度下降)。如果网格不动，那么固定结点上的物理量需要修正，以反映材料相对于固定网格(及其所代表的固定空间)的移动，而有限元法不易做到准确修正，特别是应力。

独立覆盖流形法采用的“级数解”具有公式化的表达，为物理量的准确修正创造了条件。本文借鉴流体计算中的半拉格朗日法[16]，提出对物质点进行“逆向追踪”的思路：如图4所示，当前构形中坐标为(t+Δtx,t+Δty)的空间点A被某物质点占据，可通过上一时步得到的位移级数解逆向求解该物质点在上一时步的位置A′(tx,ty)，即

t+Δtx=tx+ux(tx,ty) ,t+Δty=ty+uy(tx,ty) 。

(13)

式中ux、uy分别为x、y两个方向的位移级数，高阶情况下难以显式求解，可采用直接迭代法，以当前位置为初值，在控制好步长及整体计算稳定的情况下，一般经过几次迭代就可得到(tx,ty)的收敛解。以此坐标代入上一时步的应力和速度的级数表达式求出σ和v，作为当前位置A的物质点所携带的物理量初值，代入式(3)和式(4)中的Fg和Fσ参与积分，从而实现当前时步的拉格朗日方程的准确计算。

以欧拉观点对上述过程给出另一种解释：欧拉描述关注固定空间点上的物理量的变化，即关注不同时刻经过该空间点的不同物质点所携带的物理量的变化，而迁移项(对流项)所反映的正是这种由于物质从一点移动到另一点的空间不均匀性所引起的变化。通过“逆向追踪”手段反映这种变化，就可消除迁移项，实现材料相对于固定空间点移动的准确描述。

以上讨论的一般性情况，同样适合于各独立覆盖。不过，针对各物质点的操作，需要形成级数表达式，以便于在条形区域内“加权”，以及式(6)中的应力累积和式(7)中的速度更新：在“逆向追踪”后，根据各积分点的速度v和应力σ的修正数值，用最小二乘法得到修正后的级数表达式。另外，考虑应变二阶项计算应变及应力时，由于式(9)的复杂性，也要通过最小二乘法获得与位移同阶次的应力级数，然后再代入式(6)与tσ累积得到t+Δtσ的级数表达式。

对于变形过程中可能出现的复杂形状边界，本文采用3次样条曲线描述，按单值性分为y=f(x)和x=f(y)2种，并可设置多条曲线用于描述复杂边界。构形变化时，更新各样条曲线的控制点。

图5 正方形网格中的积分区域Fig.5 Integral region in asquare mesh

采用正方形网格，与每个时步更新后的边界轮廓进行切割操作。本文提出形成固定网格内的积分区域的新方法，如图5所示，简要说明如下：

对于某个正方形网格，边界轮廓(包括直线段或样条曲线)依次与网格的4条边求交，比如，得到图5中左、右两条边上的交点A、B，作为有效顶点；
通过判断，结点1和结点2在域内，也作为有效顶点；
各段边界线的起点或终点也有可能落在网格内，成为有效顶点；
按逆时针走向，从边1-2开始，在各边(及网格内部)依次连接有效顶点，形成积分区域1-2-B-A-1；
其他与边界不相交的正方形网格，如判断其形心在域内，则整个网格作为积分区域参与积分。

在材料边界的移动过程中，不可避免地会形成在网格内只占很小区域的所谓“小块”，如图6(a)所示的右上角网格。其中，在材料进入新网格时，一定会首先出现小块。设置一定的小块比例，比如，小块的面积为所有块体平均面积的5%以下。

小块的存在会影响方程性态，造成数值解的不稳定。考虑到独立覆盖流形法具有在不同大小的网格之间错位连接的特性，采用的处理方式如图6(b)所示，将小块所在的网格与其周边大网格合并(目前暂时取小块的长边所连接的大网格)，应力和速度采用大网格中的级数。

随着材料边界的移动，小块有两种变化的可能性：第1种如图6(c)所示，边界向下移动，小块完全消失；
第2种如图6(d)所示，边界向上移动，小块逐渐占据一定比例的正方形网格面积(比如大于平均块体面积的5%)，则成为独立网格中的积分区域参与计算，其初始的应力和速度继承原大网格中的级数。第2种情况同时说明了材料进入新网格时，如何通过小块将信息传递给新网格。

图6 小块的处理Fig.6 Processing of small blocks

通过以上的各种处理，在几何非线性问题的分步计算中，只要在各覆盖分区内采用足够的多项式级数逼近，就能保证每步计算的每个物质点(包括边界点)的位移精度，从而保证边界轮廓的准确更新、速度的准确更新和应力的准确累积，在各时步实现移动边界的准确描述。

图7 固定的正方形网格中的悬臂梁变形Fig.7 Deformations of cantilever beam in fixedsquare meshes

算例1：计算如图7所示的水平放置的悬臂梁，左端固定。该梁长10 m，截面的高度和宽度均为1 m，弹性模量E为3×105kN/m2，泊松比μ为0.2。在梁自由端的截面中点作用始终垂直向下的集中力P。在固定的正方形网格中计算其大变形，梁的上、下边缘分别采用3次样条曲线。每个独立覆盖采用3阶多项式级数，步长为6.25 kN(数值试验表明，计算精度随步长减小而提高)。典型时步的自由端中点位移及其与理论解的比较见表2，相对误差在1%左右。在计算中发现，若完全按照静力计算，则在很大变形时(P>1 000 kN)会出现计算振荡现象。除了引入初应力刚度矩阵解决此问题以外，本文采用动力松弛法求解静力问题，也能保证计算的稳定性。

表2 悬臂梁自由端的中点位移Table 2 Displacements of the mid-point at the free endof cantilever beam

算例2：图7中的正方形网格内的梁作为刚性杆件，弹性模量取大值以模拟刚性。杆件的左端中点为固定约束，从原位置开始旋转，初始角速度为1 rad/s，步长为0.01 s。各时步的自由端中点坐标见表3，采用1阶和4阶多项式级数的计算结果相同，但是否考虑二阶应变对计算结果有影响：不计二阶应变，则杆件会在旋转过程中逐渐“拉长”；
考虑二阶应变，得到与理论解几乎一致的结果。此算例验证了本文方法模拟刚体运动的准确性。

表3 匀速旋转杆件的自由端中点坐标Table 3 Coordinates of the mid-point at the freeend of rod with uniform rotation

算例3：为了更好地说明本文方法相对于常规拉格朗日法的优势，给出如图8(a)所示的半圆结构在垂直体力作用下的变形分析算例。该结构半径为4 m，底部法向约束，取弹性模量为100 kN/m2，泊松比为0.2。考虑对称性，取1/4圆建立计算模型如图8(b)所示，从圆弧的1/2处划分为两段，分别采用y=f(x)和x=f(y)2种样条曲线，左侧和底部为法向约束。采用了8×8个正方形网格覆盖其可能的变形区域，以及4阶多项式级数。

图8 在垂直体力作用下的半圆结构及其固定网格Fig.8 Semicircle in fixed meshes under verticalbody loads

图9 在不同垂直体力作用下的结构变形Fig.9 Deformations of the semicircle under variedvertical body loads

不同体力荷载f作用下的结构变形如图9所示，结构顶部竖向位移(VB)和底部外边缘的水平向位移(UA)列入表4，并与大型有限元分析软件MARC的结果对比。MARC模型划分了540个有限元网格，考虑随变形变化的荷载选项。从表4可见，在较小荷载作用时，本文计算值与MARC结果接近。随着荷载增大，MARC的有限元网格因变形而扭曲，最终在f=17 kN/m3时出现计算不收敛的情况，提示矩阵奇异。而本文方法没有网格扭曲现象，可一直计算下去，体现出方法的优越性。

表4 A点水平位移与B点垂直位移Table 4 Horizontal displacements of point A andvertical displacements of point B

本文提出应用级数对物质点进行“逆向追踪”的思路，在空间固定的网格中采用拉格朗日控制方程求解几何非线性问题，集合了拉格朗日法的跟踪物质点、控制方程简单、边界描述准确以及欧拉法的网格无扭曲的优势，避开了两种方法相应的缺陷。相比之下，任意拉格朗日和欧拉法需要设计网格运动算法，该算法通用性不强且有可能很复杂。本文采用级数公式反映物质点的相对运动，相当于公式化的物质点法，求解精度与独立覆盖流形法的常规计算无异，而物质点法的不足是用移动的质点作为积分点，相比常规有限元法的精度降低。本文的研究表明，除了无网格法的途径之外，网格类的方法也可以很好地解决特大变形下的网格扭曲问题。

与物质点法类似，采用固定网格会涉及物质点在相邻网格之间的穿越，这要求物理量在网格之间具有一定的连续性，特别是应力，否则会引起局部计算的振荡。因此，为保证每步计算的收敛性(包括应力的连续性)，本文算例多采用3至4阶多项式级数。相对于笔者前期提出的固定网格流形法而言，本文基于独立覆盖的固定网格具有任意连接和任意加密的特性，且覆盖级数升阶方便。目前正在引入自适应分析技术，在一定的精度控制下，通过加密网格(也可配合级数升阶)，可以实现每步的应力“分区级数解”在一定程度上的整体连续甚至光滑。因此固定网格并非一成不变，可在各积分点处由“整体连续”的级数公式计算应力和速度初值，从而实现在各时步的不同网格之间准确传递信息，以替代通常基于网格的传递方式。

本文方法跟踪物质点且边界描述清晰准确，便于求解随时间变化的材料非线性和接触非线性问题：通过“逆向追踪”，使得当前位置的物质点携带材料状态信息，其中当前边界上的物质点还携带边界接触状况。因此，将来有可能在空间固定的背景网格中同时进行几何、材料和接触的多重非线性分析。

致谢：感谢石根华先生的指导。感谢中央级公益性科研院所基本科研业务费项目(CKSF2019394/GC)的资助。

猜你喜欢 拉格朗级数欧拉 无穷级数敛散性的判别方法探讨科技风(2022年26期)2022-10-1019.93万元起售，欧拉芭蕾猫上市车主之友(2022年4期)2022-08-27欧拉魔盒哈哈画报(2022年1期)2022-04-19精致背后的野性 欧拉好猫GT车迷(2022年1期)2022-03-29欧拉秀玛杂记文苑(2020年8期)2020-09-09这样的完美叫“自私”杂文选刊(2019年12期)2019-12-06二重Dirichlet级数在收敛半平面内的增长性知识文库(2019年4期)2019-10-20一个非终止7F6-级数求和公式的q-模拟华东师范大学学报（自然科学版）(2019年3期)2019-06-24拉格朗日的“自私”意林·少年版(2018年22期)2018-12-05拉格朗日的“自私”新青年(2018年8期)2018-08-18

相关热词搜索：求解 网格 几何