节点梯度光滑有限元配点法1)

樊礼恒 王东东 刘宇翔 杜洪辉

(厦门大学土木工程系,福建厦门 361005)

(厦门市交通基础设施智能管养工程技术研究中心,福建厦门 361005)

引言

基于变分法或等效积分弱形式的伽辽金法和基于强形式的配点法是计算力学领域两种常用的数值计算方法[1].与伽辽金法相比,配点法具有构造简单、计算高效的特点,但是配点法需要计算数值离散形函数的高阶梯度,例如势问题或弹性力学问题等二阶问题需要用到形函数的二阶导数.同时配点法的刚度矩阵不对称,稳定性相对较差.有限元法是目前应用最为广泛的一种典型伽辽金法[1],其形函数定义在单元上,一般网格情况下仅具有C0连续性[2],在相邻单元边界处会出现梯度跳跃,因而不能直接应用于需要计算形函数高阶梯度的配点法[1].近年来,随着采用高阶光滑形函数的无网格法和等几何分析的快速发展[3-7],配点法逐渐得到了愈来愈广泛的应用,例如径向基函数配点法[8-13],移动最小二乘与再生核无网格配点法[14-18],等几何配点法[19-21]等.值得指出的是,与无网格和等几何形函数相比[1],有限元形函数形式简单,计算快捷.最近,高效伟等[22-23]利用有限元形函数在单元内部高阶连续的特点,通过在选定配点邻域内构建局部单元,提出了自由单元配点法.

就配点法的精度而言,一个突出的问题就是其计算精度依赖于形函数所采用的基函数阶次的奇偶性,即奇数阶次基函数的配点法收敛率较其基函数阶次下降一次[19,24].王东东等[24-26]系统研究了二阶和四阶问题无网格配点法奇数阶次基函数精度掉阶的原因.研究表明,配点法的精度与形函数及其梯度的一致性或完备性条件密切相关.因此,王东东等[24-25]通过引入无网格形函数的梯度光滑方法,建立了超收敛无网格配点法,并在不引入待定配点位置的情况下有效解决了奇数阶次基函数配点法的精度掉阶问题,实现了采用奇数次基函数配点法精度的超收敛提升.基于梯度光滑方法,邓立克等[27]通过线性基函数构造了二阶光滑梯度,使得采用线性基函数就可以进行四阶薄板问题的伽辽金无网格分析.然而,注意到上述工作中所讨论的形函数高阶光滑梯度构造的基础是形函数的标准一阶梯度,这对于本身具有高阶光滑特性的无网格形函数没有问题,但是由于有限元形函数的一阶梯度在单元节点和边界处并不存在,因此难以直接构造有限元形函数的高阶光滑梯度和发展相应的配点法.

为了构造有限元形函数的高阶光滑梯度,本文在无网格法的梯度光滑理论框架内[24,28],首先引入了有限元形函数在节点处和全域内的一阶光滑梯度定义,然后在此基础上建立了有限元形函数的二阶光滑梯度,进而采用有限元形函数的光滑梯度建立了一种新型的节点梯度光滑有限元配点法.与传统有限元形函数梯度不同,有限元形函数的一阶和二阶光滑梯度具有全域连续的特点.不失一般性,本文以线性形函数为例,详细研究了有限元形函数一阶和二阶光滑梯度的构造特点和一致性条件,并对节点梯度光滑有限元配点法进行了精度分析.分析表明,线性有限元形函数的二阶光滑梯度能够满足二阶一致性条件,因而对应的节点梯度光滑有限元配点法具有二次精度和超收敛特性.注意到基于梯度光滑理论[28],Liu 等[29]构造了有限元形函数的一阶光滑梯度,提出了伽辽金光滑有限元法.但本文所讨论的一阶与二阶有限元形函数光滑梯度,其光滑核函数选取和构造方式均与伽辽金光滑有限元法不同,同时着重研究的是配点法计算列式.文中通过典型算例系统验证了节点梯度光滑有限元配点法的精度和收敛性.

1 配点法控制方程

本文考虑稳态势问题和弹性力学问题等典型二阶问题,其控制方程的强形式可以表示为

其中u(x)为场变量,例如弹性力学问题的位移或热传导问题的温度,L 为空间域Ω 内的偏微分算子,Bg和Bt为强制边界Γg和自然边界Γt上的算子,b为源项,g和t分别为Γg和Γt上给定的边界值.方便起见,这里给出二维情况下的各个算子的定义.例如,对于稳态势问题,L,Bg和Bt为

其中n表示在自然边界Γt的单位外法线向量.对于平面应力弹性力学问题,上述算子可表示为

式(4)和式(6)中D为平面应力弹性本构矩阵,I为二阶单位阵,E和ν 分别代表材料的杨氏模量和泊松比.

将式(7)代入控制方程(1),并令其离散节点处的残值为零,便可以得到配点法的离散方程

其中K,d和f分别表示刚度矩阵,节点系数向量和力向量

由式(10)可见,配点法需要计算形函数的二阶梯度,而标准有限元形函数在单元边界处不存在二阶梯度,因此无法直接用于配点法.此外,虽然通过特殊构造方式可以建立一些特殊单元的高阶光滑有限元形函数,但是目前仍然缺乏构造一般单元高阶光滑形函数的通用有效方法[1].本文则是通过引入无网格法中的梯度光滑技术,建立了有限元形函数的一阶和二阶光滑梯度,进而发展了相应的节点梯度光滑有限元配点法.

2 有限元形函数的光滑梯度构造

2.1 有限元形函数光滑梯度构造理论

考虑图1 所示的有限元网格,ΩL表示第L个单元对应的空间区域,节点xI对应的有限元形函数是NI(x),NI(x)的影响域为

式中MI表示共享节点xI的单元个数.方便起见,定义下面的节点组和单元组

根据广义梯度光滑理论[4,24],有限元形函数NI(x)的一阶光滑梯度的定义为

图1 有限元形函数影响域示意图Fig.1 Illustration of the influence domains of finite element shape functions

式中φ(x,y)为梯度光滑核函数,NI,i(x)为传统的标准一阶梯度或导数,即NI,i(x)=∂NI(x)/∂xi,xi为x的第i个坐标分量.由于有限元形函数的梯度在单元边界处没有定义,因此式(14)并不能直接用于有限元形函数的光滑梯度构造.

为了构造有限元形函数的光滑梯度,这里首先引入有限元形函数在节点处的梯度值.为此,在式(14)中选取如下的梯度光滑核函数

若进一步考虑常应变单元,例如一维两节点单元、二维三节点三角形单元和三维四节点四面体单元,则NI,i(x)在每个单元内均为常数,式(16)可以简化为

其中VL表示单元ΩL的长度、面积或体积.

在式(17)定义的有限元形函数的节点梯度基础上,进一步选择梯度光滑核函数φ(x,xJ)=NJ(x),则根据式(14)对应的离散形式,可得如下的有限元形函数一阶光滑梯度

接下来,对式(18)进行求导,并用式(18)定义的一阶光滑梯度代替传统一阶梯度,最后可得有限元形函数的二阶光滑梯度

至此,已经构建有限元形函数的二阶光滑梯度.下面通过一维和二维线性有限元形函数一阶和二阶光滑梯度的构造过程来说明有限元形函数光滑梯度构造理论.

2.2 线性有限元形函数光滑梯度构造

清晰起见,首先考虑图2 所示的一维均匀有限元网格,其中h表示单元的长度,在此情况下,ΩI=,因此VI=h,=2h.有限元节点xI的形函数NI(x)可表示为

其对应的一阶梯度或导数NI,x(x)为

由式(20)、式(21)和图2 可见,有限元形函数NI(x)为C0连续,其一阶梯度NI,x(x)在单元节点xI−1,xI及xI+1 处不连续,无法直接计算二阶梯度,因而无法用于式(8)定义的配点法求解.

另一方面,根据式(14),有限元形函数在节点处的一阶光滑梯度为

图2 一维有限元形函数一阶光滑梯度构造过程Fig.2 Construction of the first order smoothed gradient of 1D finite element shape function

在此基础上,利用式(18),可得一维线性有限元形函数的一阶光滑梯度

进一步将式(22)和式(23)代入式(19),可得一维线性有限元形函数的二阶光滑梯度

图2 和图3 给出了线性有限元形函数一阶和二阶光滑梯度的构造过程.由图2 和图3 可见,虽然线性有限元形函数的梯度在节点处不连续,而且二阶梯度不存在,但由本文所提有限元形函数梯度光滑理论构造所得的一阶和二阶光滑梯度仍然是连续函数,可以直接用于配点法求解.

对于二维和三维情况,考虑三节点三角形单元和四节点四面体单元,其对应的形函数为标准的线性形函数[1],这里不再赘述.二维和三维有限元形函数的一阶和二阶光滑梯度也可按照上述步骤进行构造.图4 给出了基于平面三角形线性单元构造一阶光滑梯度和和二阶光滑梯度的过程.与一维情况相同,线性三角形单元的形函数在节点和边界处的一阶梯度均不连续,不存在二阶梯度.另一方面,基于本文的有限元形函数光滑梯度构造理论,可以方便地构造图5 所示的一阶和二阶光滑梯度.此外,图2～图5 同时表明,与有限元形函数和其标准梯度不同,有限元形函数的一阶和二阶光滑梯度与形函数本身具有相同的连续性,但光滑梯度的影响域则大于形函数的影响域.

图3 一维有限元形函数二阶光滑梯度构造过程Fig.3 Construction of the second order smoothed gradient of 1D finite element shape function

图4 二维有限元形函数一阶光滑梯度构造过程Fig.4 Construction of the first order smoothed gradient of 2D finite element shape function

图5 二维有限元形函数二阶光滑梯度构造过程Fig.5 Construction of the second order smoothed gradient of 2D finite element shape function

3 节点梯度光滑有限元配点法及其精度分析

将式(18)和式(19)定义的有限元一阶和二阶节点光滑梯度代入配点法列式(8),便形成了节点梯度光滑有限元配点法.以势问题为例,节点梯度光滑有限元配点法离散方程刚度矩阵的各个元素为

节点梯度光滑有限元配点法力向量的各个元素与式(10)相同.为不失一般性,下面首先分析线性有限元形函数光滑梯度的一致性条件,然后以势问题为例对基于线性单元的节点梯度光滑有限元配点法的计算精度进行理论研究.

3.1 线性有限元形函数光滑梯度的一致性条件

对于传统线性有限元,其形函数及其一阶梯度的一致性或完备性条件为[1]

其中xI j为节点xI的第j个坐标分量,δi j为克罗内克符号.注意到有限元形函数的一致性条件取决于形函数的阶次.例如,线性有限元形函数仅满足式(26)～式(29)定义的线性一致性条件,由于其二阶梯度不存在,故不满足更高一次的二阶一致性条件.但是对于有限元形函数的光滑梯度,接下来证明其不仅满足线性一致性条件,在均布网格情况下还满足二阶一致性条件.

对于均匀有限元网格,由形函数的周期重复性可知,节点光滑梯度满足如下条件[24-25]

其中k=0,1 或2.首先考虑k=2 的情况,利用式(32),式(33)变为

再者,对于k=1,式(33)变为

式(34)～式(36),证明了线性有限元形函数的二阶光滑梯度满足二阶梯度一致性条件.与之形成鲜明对比的是,标准线性有限元形函数不存在二阶梯度.

3.2 节点梯度光滑有限元配点法精度分析

为不失一般性,以平面势问题为例,在均匀网格条件下对节点梯度光滑有限元配点法的精度进行解析研究.节点梯度光滑有限元配点法的精度可以方便地采用如下的截断误差进行度量[24,30]

其中uI为节点xI的解析常变量,即uI=u(xI),h表示单元的特征长度,k为精度的阶次.

在式(37)中引入uI的泰勒展开形式

同时注意到在均布网格条件下,有限元形函数的光滑二阶导数具有周期性和偶对称的特点,因此当n为奇数时,有

将式(31),式(33)～式(36)及式(40)代入式(39),可得

由式(41)可以看出,节点梯度光滑有限元配点法具有二阶精度,其L2误差和H1或能量误差均为二阶收敛.与线性有限元法相比,虽然两者L2误差收敛特性相同,但节点梯度光滑有限元配点法的H1或能量误差具有高一阶次的超收敛特性.

4 算例

本节通过一维、二维及三维算例验证节点梯度光滑有限元配点法的收敛性和精度,其中FEM 和FEC 分别表示传统的伽辽金有限元法和本文所提的节点梯度光滑有限元配点法.为了和传统有限元法进行比较,算例中采用如下误差形式

其中nd表示算例的空间维度,σij和εi j分别为应力和应变分量.

4.1 一维弹性杆问题

考虑两端固定的一维弹性杆问题,其几何和材料参数为:杆长为L=1,抗拉刚度EA=1,体力b(x)=−24x2.该问题的解析解为u(x)=2x4.计算中采用图6 和图7 所示包含20,40 和80 个单元的均布和非均布有限元网格进行分析.

图6 一维弹性杆问题均布有限元离散模型Fig.6 Uniform finite element meshes for the 1D elastic rod problem

图7 一维弹性杆问题非均布有限元离散模型Fig.7 Non-uniform finite element meshes for the 1D elastic rod problem

图8 和图9 给出了一维弹性杆问题的L2和能量误差收敛情况.结果显示,节点梯度光滑有限元配点法的位移误差收敛率与有限元法相同,精度略低,但是其能量误差收敛率比有限元法高一阶,对应的精度也明显优于有限元法,很好地验证了节点梯度光滑有限元配点法的精度和收敛率.图10 为采用20 个均布和非均布单元得到的应力结果,可见节点梯度光滑有限元配点法(FEC)的应力精度显著高于传统有限元法(FEM).

4.2 二维矩形区域势问题

考虑如图11 所示的正方形势问题区域,其边长L=1,对应势问题的解析解为

图8 采用均布网格的一维弹性杆问题收敛率对比Fig.8 Convergence comparison for the 1D elastic rod problem using uniform meshes

图9 采用非均布网格的一维弹性杆问题收敛率对比Fig.9 Convergence comparison for the 1D elastic rod problem using non-uniform meshes

图10 一维弹性问题应力解对比Fig.10 Comparison of stress solutions for the 1D elastic rod problem

图11 方形区域势问题模型Fig.11 Description of the square region potential problem

图11 中所示自然和强制边界条件均由式(45)给定.图12 为该问题采用的有限元离散模型,对应的节点数和平面三角形单元数分别为289,1089,4225 和512,2048,8192.

图13 给出了该势问题的L2和H1误差收敛结果,其中节点梯度光滑有限元配点法的L2和H1误差收敛率分别为2 和1.9,与上节的理论收敛率基本一致,尤其是H1误差的收敛率和精度都明显高于有限元法,具有超收敛特性.图14 为采用512 个三角形单元计算得到的梯度解的相对误差对比,即,其中图14(a)和图14(b)分别表示有限元法和节点梯度光滑有限元法的误差分布图.由图14 的对比结果清晰可见,传统线性三角形单元的梯度解为常数且在单元边界处不连续,因此节点梯度光滑有限元配点法的梯度解误差远小于传统有限元法,表明了节点梯度光滑有限元配点法梯度解的精度优势.

4.3 悬臂梁问题

图15 所示为一平面应力悬臂梁问题,其左端为强制位移边界,右端受竖向力P=1000 作用,几何与材料参数为:长度L=4,梁截面高H=1,宽W=1,杨氏模量E=2.0×107,泊松比ν=0.3,惯性矩I=WH3/12.该问题的解析解为

图12 方形区域势问题的有限元离散模型Fig.12 Finite element discretizations for the square region potential problem

图13 方形区域势问题收敛率对比Fig.13 Convergence comparison for the square region potential problem

图14 方形区域势问题梯度解误差对比Fig.14 Error comparison of the gradient solution for the square region potential problem

图15 悬臂梁问题Fig.15 Description of the cantilever beam problem

该悬臂梁问题的有限元离散模型如图16 所示,三种有限元离散对应的节点数225,833 和3201,三角形单元数为384,1536 和6144.图17 为该算例的位移和能量误差收敛对比结果.与一维弹性杆问题的结果一致,传统线性有限元法的位移和能量误差收敛率分别为2 和1,而节点梯度光滑有限元配点法的位移误差收率也为2,但其能量误差收敛率为2,比传统有限元法高一阶,具有更好的收敛特性.此外,图17也表明节点梯度光滑有限元配点法的能量误差明显小于传统有限元法.图18 给出了采用384 个单元的剪应力解的相对误差云图,结果进一步显示出节点梯度光滑有限元配点法的高精度应力计算特点.

图16 悬臂梁问题有限元离散模型Fig.16 Finite element discretizations for the cantilever beam problem

图17 悬臂梁误差收敛率对比Fig.17 Convergence comparison for the cantilever beam problem

图18 悬臂梁问题应力解误差对比Fig.18 Error comparison of the stress solution for the cantilever beam problem

4.4 三维厚壁圆筒势问题

为了验证节点梯度光滑有限元配点法的三维适用性,最后一个算例是图19 所示的厚壁圆筒势问题,相应的几何参数为:内径为ri=1,外径为ro=2,高度H=1.该问题的解析解为

利用该问题的对称性,计算中取四分之一模型进行离散分析.图20 为该三维势问题的有限元离散情况,对应的节点数为225,1377 和9537,四节点四面体单元数为768,6144 和49 152.由于空间区域的构造特点,该问题的有限元网格具有非均布特性.图21 为采用图20 中三种有限元离散网格计算得到的L2和H1误差收敛率结果.与一维和二维算例的结果类似,有限元法和节点梯度光滑有限元配点法的L2误差收敛率都为2,节点梯度光滑有限元配点法的H1误差收敛率受非均布有限元网格影响,略低于理论值2,但其H1误差的收敛率和精度仍然明显高于传统有限元法,验证了节点梯度光滑有限元配点法求解三维问题的精度和有效性.图22 给出了768个单元对应的梯度解的相对误差对比.结果再次证明,即使对于非均布有限元离散,节点梯度光滑有限元配点法仍然在梯度解方面具有突出的精度优势.

图19 三维厚壁圆筒势问题模型Fig.19 Description of the 3D hollow cylinder potential problem

图20 三维厚壁圆筒势问题有限元离散模型Fig.20 Non-uniform finite element discretizations for the 3D hollow cylinder potential problem

图21 三维厚壁圆筒势问题收敛率对比Fig.21 Convergence comparison for the 3D hollow cylinder potential problem

图22 三维厚壁圆筒梯度解误差对比Fig.22 Error comparison of the gradient solution for the 3D hollow cylinder potential problem

5 结论

传统有限元形函数在离散区域内一般仅有C0连续性,在单元边界和节点处上不存在一阶和二阶梯度,因而不能直接用于配点法.本文提出了有限元形函数的光滑梯度构造理论,建立了具有与有限元形函数同样连续性的一阶和二阶光滑梯度,进而发展了节点梯度光滑有限元配点法.有限元形函数的一阶光滑梯度依赖于在广义梯度光滑理论框架内的一阶光滑梯度节点值的定义,以及选择有限元形函数作为核函数对一阶光滑梯度节点值进行梯度光滑.对有限元一阶光滑梯度进行求导,并用一阶光滑梯度替换有限元形函数的标准梯度,即可得到有限元形函数的二阶光滑梯度.文中详细证明了线性有限元形函数的光滑梯度不仅满足常规的一阶梯度一致性条件,而且在均布网格条件下还满足二阶梯度一致性条件.与之对应的是,线性有限元形函数的标准一阶梯度只满足一阶梯度一致性条件.正是由于有限元形函数的光滑梯度满足二阶梯度一致性条件,理论分析表明,基于有限元形函数光滑梯度的节点梯度光滑有限元配点法具有二阶精度,即L2误差和H1或能量误差的收敛率均为2,与传统伽辽金有限元法相比呈现超收敛特性.文中通过典型算例系统验证了节点梯度光滑有限元配点法在梯度或应力求解方面的超收敛特性和显著精度优势.本文的节点梯度光滑有限元配点法采用了线性有限元形函数,对于高次有限元形函数及高阶光滑梯度还需要进一步研究.