基于赫林格-赖斯纳变分原理的一致高效无网格本质边界条件施加方法1)

吴俊超 吴新瑜 赵珧冰 王东东

*(华侨大学土木工程学院,福建省智慧基础设施与监测重点实验室,福建厦门 361021)

†(厦门大学土木工程系,福建省滨海土木工程数字仿真重点实验室,福建厦门 361005)

引言

无网格法[1-7]是指基于节点信息构造离散形函数的一类数值方法的总称.该类方法通常具有高阶光滑、全域协调的形函数,形函数构造过程不依赖于网格单元的拓扑信息,适用于大变形分析[8]、薄板壳高阶问题[9]及裂纹扩展模拟[10]等.然而,与传统有限元法相比,高阶连续光滑的特点导致无网格形函数在离散节点上通常不具有插值性,在求解过程中难以直接施加本质边界条件[11-12].其次,无网格形函数通常为有理式,并且形函数影响域高度重叠,导致形函数在背景积分单元上为分段的有理式.在伽辽金法的求解过程中,传统高斯积分法无法精确数值积分由形函数组成的刚度矩阵和力向量,导致伽辽金弱形式不满足积分约束条件或称变分一致性条件[13-14],无法保证计算精度和最优误差收敛率.

为了使无网格法能够直接施加本质边界条件,许多学者构造了诸多具有插值性的无网格近似方法,例如奇异权函数法[15]、插值最小二乘法[16-17]、复变量移动最小二乘法[18]、广义移动最小二乘法[19]、变换法[20]等.然而,这类方法不是建立在变分原理基础上,并不能保证节点之间位移边界条件施加精度和无网格法的变分一致性.此外,对于满足积分约束条件的无网格数值积分方法,例如稳定节点积分法[13]、一致性积分法[21-22]、变分一致积分法[23]、嵌套子域积分法[24]、再生光滑梯度积分法[25]等,在计算过程中采用形函数的光滑梯度替代传统无网格形函数的导数,在保证无网格法的计算精度和最优误差收敛率的同时提高了计算效率.但其本质边界条件仍需要具有变分一致性的方法进行施加[26-27].

在伽辽金无网格法中,拉格朗日乘子法[1]是常采用的一种变分一致本质边界条件施加方法.该方法需要在整体刚度矩阵上增加额外的自由度离散拉格朗日乘子.当采用满足变分一致无网格数值积分方法时,拉格朗日乘子的自由度需要与光滑梯度构造过程中的数值积分点相一致.但是,过多的拉格朗日乘子自由度将增加整体刚度矩阵奇异性,不适用于高阶基函数的无网格法.Lu等[28]根据拉格朗日乘子的物理意义,采用位移自由度离散拉格朗日乘子,无需增加额外自由度,并称之为修正变分原理法.但该方法施加本质边界条件过程中的修正变分项降低了刚度矩阵的正定性,计算精度不高.Zhu和Atluri[29]将罚函数法引入无网格法中施加本质边界条件,该方法格式简洁易实现,但其形式并不满足变分的一致性,求解精度依赖人工参数,稳定性较低.Fernández-Méndez和Huerta[11]采用尼兹法(Nitsche method)施加无网格法本质边界条件,该方法可视为修正变分原理与罚函数法的结合,修正变分项保证了变分一致性,而罚函数项恢复了刚度矩阵的正定性,提高了求解的稳定性,但该方法仍然需要选择合适的人工参数来保证数值解的精度.

本文针对满足积分约束条件的伽辽金无网格法,提出了一种具有变分一致且高效的本质边界条件施加方法.该方法以赫林格-赖斯纳(Hellinger-Reissner,HR)变分原理[30]为基础,采用混合离散的方式近似弱形式中的位移和应力.其中,位移采用传统无网格形函数进行离散,而应力则在每个背景积分单元中近似为局部多项式.此时,近似应力中多项式的系数可在背景积分单元中采用位移自由度表示,无需额外增加整体自由度.不仅如此,应力的系数表达式中包含再生光滑梯度表达式,在再生光滑梯度的理论框架下,光滑梯度构造过程中可采用优化的数值积分点以减少全域无网格形函数的计算量,计算效率高.整体的离散平衡方程具有与尼兹法相类似的形式,其中尼兹法中的修正变分项分别引入了再生光滑梯度和无网格形函数离散应力和位移,整体求解过程无需计算传统无网格形函数的导数.而在赫林格-赖斯纳变分原理的框架下,弱形式中自带稳定项,无需增加额外的稳定项,且稳定项中不存在人工参数,消除了对人工参数的依赖性.文中通过二维弹性力学算例验证了所提方法的精度、效率和误差收敛性.

1 赫林格-赖斯纳变分原理

不失一般性,在求解域 Ω 内考虑如下弹性力学控制方程

式中,W为弹性体的余能密度函数,其与应力之间的关系为

其中Cijkl为四阶弹性张量.

对上式进行变分可得与之相对应的弱形式

2 位移离散与再生核近似

在赫林格-赖斯纳变分原理的弱形式中,位移与应力可分别采用不同的近似方案进行离散.这里,位移采用基于再生核近似[29]的无网格形函数进行离散.首先,在求解域 Ω 上布置一系列无网格节点.此时,位移分量ui的近似表达式可表示为

式中,diI为无网格节点xI上的位移分量节点系数,ΨI为与之相对应的无网格形函数,根据再生核近似理论[2],无网格形函数具有如下表达式

其中,p为基函数向量,以二维p阶基函数为例,此时基函数向量为

φ为核函数,其影响域也决定了无网格形函数的影响域.在二维情况下,核函数的影响域通常为圆形或矩形,本文采用基于三次样条函数、具有矩形影响域的核函数

式中,rx=|xI-x|/sx,ry=|yI-y|/sy,sx和sy分别为x和y方向的影响域尺寸.φ为三次样条函数[1]

无网格形函数表达式(8)中的c为修正函数向量,该表达式可通过满足再生条件[2]确定

将式无网格形函数表达式(8) 代入再生条件(13)中,可得

式中M称为矩量矩阵

将式(14)、式(15)代入式(8)中,可得到再生核无网格形函数的表达式

图1为二维二次基函数无网格形函数,从图中可以看出,无网格形函数在全域上连续光滑.但在无网格节点处,无论节点位于域内或边界处,形函数都不具有插值性,这是伽辽金无网格法难以施加本质边界的主要原因.同时,相较于有限元形函数,无网格形函数并不具有显示表达式,形函数计算复杂耗时,无网格形函数的导数更是如此.

图1 无网格形函数Fig.1 Meshfree shape functions

3 应力离散与再生光滑梯度近似

表述方便起见,引入沃伊特标记(Voigt notation),将近似应力 σh采用矩阵向量形式改写为

4 赫林格-赖斯纳变分原理下的本质边界条件施加方法

基于赫林格-赖斯纳变分原理的弱形式(4)中已考虑了本质边界条件.前两节已分别采用再生核无网格形函数和背景积分单元上的分段多项式近似式(4)中的位移和应力,进一步将近似的位移分量表达式(7)和应力分量表达式(17)代入弱形式(6)中,有

其中,δd为位移变分节点系数向量,f为外力向量.将式(19)代入式(33)中,式(33)等式左边可化简为

式中最后一个等式引入了关系式(20)～式(26)进行化简,详细的推导过程可参考文末附录A.将简化后的式(34)代入式(33),即可得到离散的平衡方程

从式(35)和式(36)可以看出,在赫林格-赖斯纳变分原理的框架下,本质边界条件的施加过程与尼兹法[11]具有相类似的格式.清晰起见,这里给出基于尼兹法的无网格法离散平衡方程

其中,K和f与式(36) 中的刚度矩阵K和力向量f相同.Kv和fv是基于变分原理施加本质边界条件产生的刚度和力向量,两者保证了边界条件施加方法的变分一致性,其具体表达式为

再者,式(40)中的Kp和fp对应于尼兹法中采用罚函数法施加稳定项的刚度矩阵和力向量

式中 1为单位对角矩阵

该部分通过设置合适的罚因子 α[31]来保证刚度矩阵的正定性,提高求解的稳定性.

对照本文所提方法和尼兹法的无网格离散方程,即式(35)和式(40),可见式(35)中的对应于尼兹法中利用变分原理施加本质边界条件,但其应变离散过程采用式(25)中的光滑梯度替代传统无网格形函数梯度,而刚度矩阵的正定性则是通过施加得以保证.与尼兹法相比,本文所提方法的稳定项表达式相对较为复杂,但在求解全过程中,无需计算复杂耗时的无网格形函数梯度,整体过程计算效率优于尼兹法,且稳定项也不包括任何人工参数.

在引入数值积分求解式(35)的过程中,为保证离散弱形式的变分一致性,数值积分在求解的各个过程也需保持一致.如图2 所示,所提方法需要两套数值积分点,和中的数值积分方案需要与,,f,和保持一致,为了保证理论误差收敛率,该数值积分方案的积分精度在文献[25]有详细讨论.另一方面,变分一致性要求式(20)中G与式(35)中的K采用相同数值积分方案.在所提方法中,由于再生光滑梯度已满足积分约束条件,再生光滑梯度也为低阶多项式,此时数值积分方案采用传统有限元中对应阶次的高斯积分即可保证数值精度.

图2 背景积分单元示意图和优化的数值积分方案Fig.2 Illustration of background integration cells and optimized quadrature rules

此外,构造再生光滑梯度需要计算积分点处无网格形函数.为了减小形函数的计算量、提高计算效率,这里采用无网格再生光滑梯度积分法[25],利用积分点在背景积分单元间的共享特性,从全域上优化整体求解过程中积分点的数量,提高计算效率.本文采用二维二次、三次基函数的优化数值积分方案,具体的数值积分点的位置与权重详见附录B.

5 数值算例

本节首先通过分片试验,验证了采用传统高斯积分法和再生光滑梯度积分法的不同本质边界条件施加方法下是否满足积分约束条件.然后,通过典型的悬臂梁问题和带孔无限大平板问题来验证采用二次和三次基函数时,所提算法的计算精度和效率.在精度分析中,采用如下的位移误差(L2-Err)和能量误差(He-Err)

表述清晰起见,下面的讨论中用“GI”和“RKGSI”表示高斯积分法和再生光滑梯度积分法,用“penalty”、“LM”和“Nitsche”表示罚函数法、拉格朗日乘子法和尼兹法这三种施加本质边界条件的方法.例如,本文所提的基于赫林格-赖斯纳变分原理的本质边界条件施加方法简写为“RKGSI-HR”.为了保证拉格朗日乘子法求解的稳定性,拉格朗日乘子均采用线性有限元形函数进行离散.当采用二次基函数时,GI 的域内积分采用13 点高斯积分,边界积分采用3 点高斯积分;当采用三次基函数时,GI 采用16 点高斯积分,边界积分采用5 点高斯积分.

5.1 分片试验

首先,采用线性、二次和三次弹性力学分片试验验证不同无网格本质边界条件施加方案的变分一致性.分片试验的求解域为边长等于1 的正方形,求解域的四边施加本质边界条件,边界条件的预设值根据如下所示的精确解施加

其中,n=1,2,3 代表线性、二次和三次分片试验.分片试验采用如图3 所示的非均布节点离散求解域,针对二次基函数的无网格近似,采用线性和二次分片试验进行测试,核函数相对影响域为2.5;而三次基函数的无网格近似采用二次和三次分片试验,核函数的相对影响域为3.5.

图3 分片试验无网格离散模型Fig.3 Meshfree discretization for the patch test

表1、表2 分别为具有二次、三次基函数的无网格法分片试验结果,从表中可以看出,传统高斯积分法由于不满足积分约束条件,即使采用高阶高斯积分也不能通过分片试验.采用满足积分约束条件的再生光滑梯度积分法时,罚函数法不具有变分一致性,无法通过分片试验.而拉格朗日乘子法由于其拉格朗日乘子采用线性形函数进行离散,无法与再生光滑梯度相匹配,也不能通过分片试验.当采用尼兹法和本文所提基于赫林格-赖斯纳变分原理的本质边界施加方法时,RKGSI 可通过分片试验.

表1 二次基函数无网格法分片试验结果Table 1 The results of patch test with quadratic basis functions

表2 三次基函数无网格法分片试验结果Table 2 The results of patch test with cubic basis functions

5.2 悬臂梁问题

考虑图4 所示悬臂梁问题,悬臂梁求解域的长和宽分别为L=48,D=12,悬臂梁杨氏模量E=3×106、泊松比 ν=0.3.悬臂梁的右端沿y轴正方向施加荷载P,而左端为固定支座.根据平面应力假设和圣维南原理,该问题的解析解[32]为

图4 悬臂梁问题模型Fig.4 Description of the cantilever beam problem

与之相应的应力分量为

为保证一致性,该问题的本质边界条件和自然边界条件分别根据解析解式(47)和式(48)进行施加.

悬臂梁求解域采用图5 所示的四个疏密程度不同的节点离散研究计算误差收敛特性,算例采用具有二次基函数的无网格近似,相应的核函数相对影响域为2.5.图6为该问题的位移误差和能量误差收敛率结果,结果表明所提RKGSI-HR 方法的精度与RKGSI-Nitsche 相当,均达到理论收敛率,但RKGSIHR 方法无需人工参数.传统高斯积分法、RKGSILM和RKGSI-penalty 均不具有变分一致性,不能达到理论误差收敛率.

图5 悬臂梁问题节点离散Fig.5 Meshfree discretizations of the cantilever beam problem

图6 悬臂梁问题误差对比Fig.6 Error comparison for the cantilever beam problem

图7为该问题的效率对比,其中图7(a)为整体求解过程中的效率分析.从图7(a) 中可以看出,RKGSI 整体的计算效率要优于GI,GI 所用CPU 耗时大概为RKGSI 的2.5 倍.这是由于RKGSI 整体无需计算传统无网格形函数导数,同时采用优化的数值积分点减小形函数的计算量,两方面共同作用的结果.此外,从该图可以看出拉格朗日乘子法施加边界条件时,整体的计算效率要略高于其他方法,这是由于拉格朗日乘子法需要增加自由度离散拉格朗日乘子,导致整体线性代数的计算量增加.为了更好说明本文所提RKGSI-HR 的计算效率,图7(b)对比了不同本质边界方法施加过程的效率对比.该图对比了施加过程计算形函数及梯度和组装相应的刚度矩阵和力向量所用时间,从图中可以看出,罚函数法和拉格朗日乘子法由于计算过程仅需要无网格形函数本身,两种方法的计算形函数耗时相同且优于其余两种方法.相比于传统尼兹法,基于赫林格-赖斯纳变分原理施加本质边界条件过程也无需计算形函数梯度,但需要构造相应的光滑梯度,此部分所用施加大约为罚函数法和拉格朗日乘子法的1.6 倍.而在组装力向量方面,赫林格-赖斯纳变分原理法和尼兹法具有相当的效率.拉格朗日乘子法由于采用有限元形函数离散拉格朗日乘子,其计算效率最高.需要指出的是,罚函数法和拉格朗日乘子法并不能保证理论误差收敛率,且从整体看拉格朗日乘子法的计算效率并不优于其他方法.综合本算例的精度分析和效率分析可得,RKGSI-HR 不仅保证计算精度和最优理论收敛性,并且相较于传统尼兹法具有更高的计算效率.

图7 悬臂梁问题效率对比Fig.7 Efficiency comparison for the cantilever beam problem

5.3 带孔无限大平板问题

考虑图8 所示带孔无限大平板问题,其中板的中心有一半径为a=1 的小孔,平板的无穷远处沿x轴方向施加均布荷载T=1000.平板材料参数为:杨氏模量E=3×106、泊松比 ν=0.3.根据Michell解[32]可得该问题的解析解为

其中 κ和µ为

与之相对应的应力表达式为

根据带孔无限大平板的对称性,如图8 所示,取边长为b=5 的四分之一方形域进行研究,求解域的左端与下端(图中红线边界)约束法向位移,求解域上端、右端和圆孔边界(图中蓝线边界)根据解析应力表达式(51)施加自然边界条件.该问题的误差收敛率测试采用图9 所示的112 个、403 个、1525 个、5929 个节点进行离散,无网格近似采用三次基函数,核函数的相对影响域为3.5.

图8 带孔无限大平板问题模型Fig.8 Description of the plate with a hole problem

图9 带孔无限大平板问题无网格离散Fig.9 Meshfree discretizations of the plate with a hole problem

图10为位移和能量误差收敛的结果,结果再次验证了所提赫林格-赖斯纳变分原理法能够保证理论误差收敛率.同样地,具有变分一致性的RKGSINitsche 也达到了理论误差收敛率,而传统高斯积分法即使在节点较为稀疏时,也难以保证收敛效果.该问题中,虽然RKGSI-LM 不能满足变分一致性,但拉格朗日乘子法保持着较高的计算精度,其误差并未影响该方法达到理论误差收敛率.图11为该问题的效率对比,与同样具有变分一致性的RKGSINitsche 方法相比,RKGSI-HR 无需计算传统形函数导数,具有更高的计算效率.最后,图12为该问题的应力云图,图中清晰表明具有RKGSI 的计算精度要高于GI,GI 的应力云图均出现不同程度的振荡.

图10 带孔无限大平板问题误差对比Fig.10 Error comparison for the plate with a hole problem

图11 带孔无限大平板问题效率对比Fig.11 Efficiency comparison for the plate with a hole problem

图12 带孔无限大平板问题 σxx应力云图Fig.12 Contour plot of stress σxxfor the plate with a hole problem

6 结论

本文以赫林格-赖斯纳变分原理为基础,提出了一种满足积分约束条件的变分一致高效本质边界条件施加方法.该方法采用混合离散近似赫林格-赖斯纳变分原理弱形式中位移和应力,其中位移采用传统无网格形函数进行离散,而应力则在每个背景积分单元上近似为对应阶次的多项式.在赫林格-赖斯纳变分原理的框架下,该方法的离散平衡方程具有与传统尼兹法相类似的格式,可视为与再生光滑梯度积分法相配套的新型尼兹法.与传统尼兹法相比,所提方法的修正变分项采用无网格形函数和再生光滑梯度进行混合离散,在保证了变分一致性的同时避免了复杂耗时的形函数导数计算,明显提高了计算效率;而对应于尼兹法中的稳定项则直接源于赫林格-赖斯纳变分原理弱形式,无需额外增加稳定项,更重要的是稳定项中不包含任何人工参数,有效消除了尼兹法中的人工参数依赖性问题.文中通过典型算例系统地验证了所提基于赫林格-赖斯纳变分原理施加本质边界条件方法的变分一致性、计算精度和计算效率.

附录A

本附录中详细推导了第4 小节式(34)中最后一个等式的推导过程,该等式为

引入关系式(20)～式(26),式(A1)中等式的左端各项简化过程如下所示

综合式(A2)～式(A6),可得等式(A1)成立.需要注意的是,从式(A4)的化简过程可以看出,的表达式虽然直观上看不对称,但从原表达式中可知其实是对称矩阵.

附录 B

本附录中列出了本文所提方法在数值实现过程中所采用的积分点位置和权重,其中附表B1和附表B2 分别为二次基函数情况和三次基函数情况的优化数值积分方案.表中“”为KIJ,G所采用的数值积分点,所采用的数值积分点.并且 ξ,η和γ为三角形参数空间坐标,w和wB分别为三角形域内积分权重和边界积分权重.

表B1 二次基函数无网格法优化的数值积分方案Table B1 The optimized quadrature rules for quadratic basis function

表B2 三次基函数无网格法优化的数值积分方案Table B2 The optimized quadrature rules for cubic basis function