无网格稳定配点法及其在弹性力学中的应用

王莉华， 刘义嘉， 钟 伟， 钱志浩

(同济大学 航空航天与力学学院，上海 200092)

1 引 言

无网格法[1-4]在数据输入时不需要提供单元连接信息，即节点之间不受网格结构限制，很大程度上节省了建模的时间和成本，而且可以在计算中根据需要改变节点的位置而不存在网格畸变问题，因此在大变形、高速碰撞、断裂破坏、金属成型以及微观粒子运动等复杂问题分析中具有明显优势，常应用到一些传统的数值计算方法(如有限元法和边界元法等)无法很好解决和尚未触及的领域。常用的无网格法主要分为基于伽辽金法的弱形式和基于配点法的强形式两类。伽辽金型无网格法主要包括扩散单元法DEM[5]、无网格伽辽金法EFG[6]、重构核粒子法RKPM[7]、hp云团法[8]、单位分解法PUM[9]、无网格局部彼得洛夫-伽辽金法MLPG[10]、径向点插值法RPIM[11]和光滑粒子伽辽金法SPG[12]等。配点型无网格法主要包括径向基函数配点法RBCM[3-15]、分区径向基函数配点法SRBCM[6-18]、局部径向基函数配点法[19]、径向点插值配点法RPICM[20]、重构核近似配点法RKCM[21,22]、有限点法FPM[23]和hp无网格云团法[24]等。

伽辽金型无网格法具有精度高和稳定性好的优点，但是由于无网格形函数通常是有理式，且形函数影响域与背景积分网格通常不重合，所以伽辽金型无网格法通常难以实现准确的数值积分。虽然稳定节点积分SCNI[25-27]等技术能够提高积分精度，但是积分过程较为复杂，尤其是构建高阶稳定节点积分显著增加了问题求解的复杂度。因此伽辽金型无网格法的效率较低，进行大规模计算非常耗时。配点型无网格法构建模式简单，不需要积分，所以效率非常高，而且结合一些特定的形函数还可以获得超收敛率[15]。然而，由于采用单个配点方程代表附近求解区域的信息而不对整个区域积分，其在求解一些复杂问题时精度和稳定性往往较差[3]。

本文针对伽辽金型和配点型无网格法的缺点，介绍一种改进的配点法,即无网格稳定配点法，该方法采用重构核函数作为近似函数，构建规则子域进行准确积分。该方法在构造近似场函数和数值积分中不需要背景网格，是一种真正的无网格法。其仅在规则子域积分，而不是像伽辽金法一样在整个求解区域积分，既保留了配点法效率高的特点，又具备伽辽金型无网格法精度高和稳定性好的特点，而且还兼具有限体积法满足局域离散方程守恒的特点。与伽辽金型无网格法相比，无网格稳定配点法的计算效率更高；与配点型无网格法相比，本文方法稳定性好且精度高。该方法兼具效率高、精度高、稳定性好和局域离散方程守恒的特点，可以应用于固体力学问题的求解，未来可将其进一步应用于流体力学和流固耦合问题分析。

2 重构核近似

重构核近似和移动最小二乘是无网格法中最常用的两类形函数，当采用多项式作为基函数时，移动最小二乘近似和重构核近似在本质上是等价的。在重构核近似中，假定近似函数可以表示为

(1)

式中c(x)是需要确定的未知矢量，HT(x-xI)是一组基矢量，其表达式如下，

HT(x-xI)=[1,x-xI,y-yI,…，

(x-xI)m(y-yI)n,…,

(y-yI)p]

(m+n≤p)(2)

(3)

式中rx=|x-xI|/θx,ry=|y-yI|/θy,θx和θy分别表示x和y方向影响域的大小，通常可取θx=θy。核函数可取如下三次B样条函数，

(4)

式中r=‖x-xI‖/θ。为了保证数值求解的精度和稳定性，形函数需满足如下的p阶一致性条件，

(m+n≤p)(5)

式中N为影响域内离散点的个数。式(5)可改写为

(6)

将式(1)代入式(6)可得未知矢量为

c(x)=M-1(x)H(0)

(7)

式中M(x)称作矩量矩阵，表示为

(8)

将式(7)代回式(1)可得重构核近似的表达式为

(9)

其一阶和二阶导数为

(α=x,y)(10)

(α=x,y;β=x,y)(11)

(12)

(13)

3 重构核配点法

在重构核配点法RKCM(Reproducing Kernel Collocation Method)[21,22]中，采用重构核近似作为近似函数，基于配点法构建离散方程。不失一般性，考虑弹性力学问题如下，

σi j,j+bi=0inΩ

(14)

σi jnj=hionΠ

(15)

ui=gionΓ

(16)

式中σi j=Ci j k lεk l为应力张量，Ci j k l=λδi jδk l+μ(δi kδj l+δi lδj k)为弹性张量，其中拉梅常数λ=νE/(1+ν)(1-2ν),μ=E/2(1+ν)，E为弹性模量，ν为泊松比，εi j=u(i,j)=(ui,j+uj,i)/2为应变张量，ui为位移，Ω为求解域，Π和Γ分别是自然和本质边界，bi为体力，hi和gi为给定的面力和位移，nj为自然边界Π的外法线方向。

将位移用重构核近似函数离散为

(17)

式中a=[a1,…,aN]T,aI=[ax I,ay I]T

(18)

Φ=[Φ1,…,ΦN],ΦI=ΨII

(19)

(20)

式中

(21)

(22)

(23)

(24)

4 重构核粒子法RKPM

重构核粒子法RKPM(Reproducing Kernel Particle Method)[7]采用重构核近似作为近似函数，结合伽辽金法构建离散方程。在伽辽金法中，引入权函数wi，并在整个问题区域内对控制方程(14)积分得

(25)

采用分步积分并引入散度定理，式(25)可表示为

(26)

考虑自然边界条件方程(15)，并代入应变表达式得

(27)

代入近似函数(17)并定义权函数的近似函数为w≈Ψ，式(27)可改写为

(28)

(29)

(30)

(31)

本质边界条件(16)可以通过拉格朗日乘子法[6]和罚函数法[28]等方法施加，本文采用拉格朗日乘子法施加本质边界条件。

在数值计算中，无网格法可以借鉴有限元法的数值积分模式，将求解区域离散成相连而互不重叠的背景积分网格Ωl(l=1,2,…,Nl)，在背景积分网格上对刚度矩阵采用高斯积分得

(32)

(33)

同理，将自然边界离散成相连而互不重叠的背景积分网格Πb(b=1,2,…,Nb)，对载荷矩阵采用高斯积分得

(34)

(35)

5 无网格稳定配点法

在无网格稳定配点法[29]中，首先在求解域及其边界上布置若干个离散点，以每个离散点为中心构建一子域，域内子域通常可选用正方形，在子域内积分如下。

(36)

(37)

(38)

(39)

(40)

将重构核近似位移离散函数(17)代入式(36～38)，可得

(41)

(42)

无网格稳定配点法中，在规则的四边形子域内进行高斯积分，可表示为

(43)

(44)

(45)

无网格稳定配点法和典型的数值方法的比较列入表1，发现稳定配点法具有如下优点。

表1 典型数值方法比较Tab.1 Comparison of the numerical methods

(1) 构建近似函数和积分均不需要背景网格，是真正的无网格法。

(2) 在子域内容易实现高阶准确积分，积分精度高，效率高。

(3) 积分子域只和离散点的位置有关，在变形前后均保持规则形状。

(4) 满足局域离散方程守恒，这对于提高流体力学问题的数值计算精度和稳定性尤为重要。

另外，需要特别说明的是，无网格稳定配点法和无网格局部彼得洛夫-伽辽金法(MLPG)都是在子域内进行积分。不同的是，无网格稳定配点法采用的是强形式直接积分；而无网格局部彼得洛夫-伽辽金法采用的是弱形式，引入散度定理，由于采用子域积分，需采用特定的试函数使得域内子域自然边界条件积分为0以简化计算，而且其和边界交界的子域积分也需经过特殊处理。此外，无网格局部彼得洛夫-伽辽金法的子域在边界上可能会截断，影响计算精度，而无网格稳定配点法的子域可以超出计算域，更为自由。因此，无网格局部彼得洛夫-伽辽金法的计算要比无网格稳定配点法更为复杂，且求解精度受积分精度影响，而无网格稳定配点法更容易实现准确积分。

6 数值算例

在数值计算中，采用如下的L2和H1误差来评估算法的精度

(46)

(47)

(i=1,2,…,Nd)(48)

式中RNd(i)∈(-1,1)为服从正态分布的随机数，Nd为随机数的个数。定义相对误差如下，

(49)

(50)

6.1 一维杆问题

考察如图2所示左端固定的一维杆问题，其控制方程和边界条件如下。

图1 稳定配点法中点的布置

图2 一维杆

EAu,x x+f(x)=0,u(0)=0,EAu,x(L)=0

(51)

图3和图4(相对误差定义见式(49，50))比较了采用重构核配点法、重构核粒子法和稳定配点法的精度、收敛性和稳定性，其中重构核粒子法和稳定配点法均采用两点高斯积分，结果表明稳定配点法的精度、收敛率和稳定性均优于重构核配点法和重构核粒子法。图5比较了三种方法的刚度矩阵条件数，其中稳定配点法的条件数最低，这也使得其具有更好的稳定性。从图6的计算时间比较可以看出，稳定配点法和重构核配点法的效率均高于重构核粒子法。

图3 一维杆问题的数值计算结果

图4 一维杆问题的稳定性分析

图5 一维杆问题的刚度矩阵条件数

图6 一维杆问题的计算时间

6.2 空心圆筒问题

如图7所示，无限长空心圆筒受内压为P0=500 N，外压为P1=500 N，圆筒内径为r0=4 m，外径为r1=10 m。弹性模量E=2×107Pa，泊松比ν=0.3。由于对称性，仅采用圆筒的1/4作为数值计算模型进行离散。其控制方程和边界条件可表示为

图7 空心圆筒问题

(52)

式中ti=σi jnj。此问题的解析解为

(53)

图8和图9的数值计算结果再次表明，稳定配点法的精度和稳定性均优于直接配点法(重构核配点法)和伽辽金型无网格法(重构核粒子法)，而且其对应的刚度矩阵条件数更低，如图10所示，这也有利于提高数值结果的稳定性。图11的计算时间比较表明，稳定配点法的效率略低于直接配点法，明显优于伽辽金型的重构核粒子法。

图8 空心圆筒问题的数值计算结果

图9 空心圆筒问题的稳定性分析

图10 空心圆筒问题的刚度矩阵条件数

图11 空心圆筒问题的计算时间

7 结 论

本文介绍了一种新的无网格法,即无网格稳定配点法，该方法不需要网格构建近似函数，也不需要背景网格进行积分，是一种真正的无网格法。由于子域都是规则的正方形，采用高斯积分可以实现准确积分，提高了算法的精度和稳定性，而且子域积分有利于降低离散矩阵的条件数，进一步提高了算法的稳定性。由于对应积分子域只和配点位置有关，和问题求解区域形状以及变形均无关，对于大变形和局部高梯度等复杂问题[18]同样适用，对于任何问题，积分子域都是规则子域而不会发生变形。相较于直接配点型无网格法，稳定配点法通过准确子域积分提高了精度和稳定性，积分非常高效，仍然保持配点型无网格法的高效性。相较于伽辽金型无网格法需要在全域内积分，而且难以实现准确积分，稳定配点法只需要在局部子域进行积分，而且通过传统高斯积分就可以实现高阶准确积分，所以比伽辽金型无网格法具有更高的效率。弹性力学算例的数值计算结果展示了该方法精度高、效率高和稳定性好的优点。由于稳定配点法的基本方法是子域法，与有限体积法一样满足局部离散方程守恒，因此其在流体力学和流固耦合分析中具有较大潜力，未来将进一步拓展这些方面的应用研究。