基于三维数值流形法的位移边界处理方法

杨石扣，艾华东

(江西理工大学 土木与测绘工程学院，江西 赣州 341000)

1 研究背景

数值流形法(Numerical Manifold Method,NMM)是Shi[1]在关键块体理论(Key Block Theory,KBT)和非连续变形分析(Discontinuous Deformation Analysis,DDA)的基础上提出的一种具有双重覆盖系统的数值计算方法，可以统一处理连续和非连续变形问题。有限元法处理位移边界条件较为方便，然而对于数值流形法[2]、广义有限元法[3]和无网格法[4]等数值分析方法而言，其处理方法往往存在各自缺陷[5]。目前，数值流形法和无网格法处理位移边界条件常采用罚函数法和Lagrange乘子法[4,6-7]。采用罚函数法得到的整体刚度矩阵具有对称、正定且实施方便等特点，且数学网格一般与边界不匹配，故而在数值流形法中应用较为广泛[5]，但罚函数法选择合适的罚值困难，取值过大则影响方程组工作性态，取值过小则计算误差较大。Zhang等[8]提出采用位移约束方程法施加本质边界条件，并应用于无网格法，Wu等[9]则提出采用BFCM(Boundary Flux Collocation Method)处理无网格法中的本质边界问题。蔡永昌等[4,10]在无网格法高斯数值积分点上对刚度矩阵和荷载向量进行修正，提出了一种准确施加位移边界和材料不连续条件的方法，具有实施简单、稳定和求解精度高等优点。Mihara等[11]将罚函数法应用于混合虚功原理，提出了一种混合型罚函数法(Hybrid-Type Penalty Method,HPM)，可方便地进行非线性分析。Kourepinis等[12]针对高阶数值流形法的边界协调问题，使用投影矩阵提出了一种边界处理方法，可方便地处理DDA、NMM和XFEM中的边界问题。Zheng等[13]和李伟等[14]提出了一种精确的边界施加方法，可方便地处理本质边界、材料接触边界和渗流边界。苏海东等[5]分析了独立覆盖流形法中本质边界条件不易施加的原因，并提出了改进边界覆盖函数和直接设定独立覆盖函数2种方法，可严格满足边界条件。对于地震响应和固体振动分析，Yang等[15]、Wei等[16]和Wang等[17]提出了相应的边界设置方法和控制方程。

目前大多数关于三维数值流形法的研究[18-25]，主要采用罚函数法施加位移边界条件，但对于沿着某一方向施加边界位移的情况未作具体讨论，且对于多步加载，计算误差累积效应将会更加明显。本文通过修正三维数值流形法控制方程中的位移边界部分，推导出含沿着某一方向施加边界位移条件和相应的分步加载情况的计算公式，以期扩大控制方程在位移边界处理上的适用范围和减少多步加载误差累积效应。通过2个典型算例的对比分析，验证了该方法的准确性。

2 三维数值流形法的基本概念

图1 数值流形法的覆盖系统和流形单元Fig.1 Covers and manifold elements for NumericalManifold Method

三维数值流形法在数学覆盖Mi上定义权函数φi(x)，权函数φi(x)满足式(1)所示的条件。

(1)

(2)

传递系数相应取值为

(3)

(4)

CT(x)=

(5)

式中n为多项式次数。

总位移函数u(x)可以写为

(7)

在实际应用中，数值计算均基于单元，故改写式(7)可得单元e内部任一点x的位移函数ue(x)为

(9)

(10)

其中di的形式为

由于本文采用六面体数学网格，故矩阵N取有限元六面体单元形函数，并采用高斯积分[25]。

3 位移边界处理方法

(12)

式中Ω是以边界Γt和Γu所形成的物理域。

(13)

将式(8)代入式(13)并考虑虚覆盖系数δdT的任意性，可得

Kd=R。

(14)

式中K和R分别为刚度矩阵和荷载向量，其形式分别为

(15)

(16)

式中:B为应变转换矩阵；D为弹性矩阵。

考虑到计算中加载步对位移影响的修正，由第k加载步至第k+1加载步的位移增量(见图2)为

图2 沿某方向tn施加边界位移示意图Fig.2 Sketch of boundary displacement applied alonga direction tn

(17)

其中，

(18)

(19)

将式(19)代入式(16)可得

(20)

4 数值算例和结果

根据第3节的方法和

公式，编写了相应的C++程序。为了验证程序和文中所述方法的可行性和正确性，选取2个典型算例进行讨论。

4.1 矩形薄板受分布力荷载作用

算例1为一个宽度L=10 m、高H=20 m和厚W=1 m的薄板[26]，如图3所示，左边及下边受连杆支承，右边及上边分别受均布压力q1=20 kN/m2及q2=10 kN/m2作用，不计体力。其位移场u和v的解析解[26]分别为：

图3 矩形薄板受分布力荷载作用Fig.3 Rectangular thin plate subjected todistributed load

(21)

计算时取薄板弹性模量E=300 MPa，泊松比μ=0.20。采用10×1×20网格数，则初始物理覆盖数为462，流形单元数为200，如图4所示。在x=0.001 m、沿z方向和z=10 m、沿x方向分别布置监测点，计算荷载q1和q2一步加载完成，且取罚函数s=1 000E，并将计算结果和解析解绘于图5。比较后可知，文中方法的计算结果与解析解相吻合，误差较小，均控制在0.1%以内。

图4 数值模型网格划分Fig.4 Meshing of numerical model

图5 沿z和x方向的水平位移和竖直位移Fig.5 Horizontal and vertical displacements alongz direction and x direction

为了说明文中方法的适应性，将计算模型绕x、y和z轴分别旋转α=30°、β=45°和γ=30°，新模型如图4(b)所示。考虑两种罚函数s=1 000E和s=10E，在式(20)中按照考虑位移边界误差修正(M)和不考虑位移边界误差修正(NM)位移边界误差修正2种情况，将计算荷载q1和q2加载步分为1、2、5和10等情况分别进行计算。在x=0.001 m、沿z方向和z=10 m、沿x方向分别布置监测点，并将计算结果绘于图6。比较后可知，罚函数较大的计算结果，位移误差较小；考虑位移边界误差修正比不考虑位移边界误差修正的计算精度高；随着分步加载数目的增多，累积计算误差逐渐增大，对不考虑位移边界误差修正的误差累积效应更加明显，对考虑位移边界误差修正的影响不大。

图6 不同加载步数、不同罚函数和是否考虑位移边界误差修正时的位移误差Fig.6 Displacement errors in cases of different loading steps and penalty functions and whether the modification of displacement boundary error is considered

4.2 矩形薄板受位移荷载作用

算例2为一个宽度a=5 m、高度b=5 m和厚度W=1 m的薄板，如图7所示，左右两边及下边均被固定，上边的位移给定为u=0，v=-η(1-x2/a2)，不计体力。当a=b时，其位移场v的解析解[26]为

图7 矩形薄板受位移荷载作用Fig.7 Rectangular thin plate subjected to displacement load

(22)

计算时取薄板弹性模量E= 300 MPa，泊松比μ=0.20，罚函数s=1 000E。采用20×1×10(见图8)、40×1×20和80×1×40共3种网格数分别进行计算。

图8 数值模型网格划分Fig.8 Meshing of numerical model

监测点坐标为(-2.5，0，2.5)，将位移v的计算结果绘于图9。比较后可知，文中方法的计算结果与解析解相吻合，误差控制在1.21%以内；随着计算网格的加密，计算误差逐渐减小，计算精度逐渐提高。

图9 监测点在不同网格密度下的位移误差Fig.9 Displacement errors of a monitoring point under different grid densities

5 结 论

(1)本文提出了一种基于三维数值流形法的位移边界处理方法，并编写了相应的C++程序。此法通过修正传统三维数值流形法控制方程中的位移边界部分，推导出含沿着某一方向施加边界位移条件和相应的分步加载情况的计算公式，扩大了控制方程在位移边界处理上的适用范围和减少多步加载误差累积效应。

(2)2个算例的计算结果表明：文中方法的计算结果与解析解相吻合，修正公式适用于不同方向施加边界位移的情况，具有较强的适应性；罚函数较大的计算结果，位移误差较小；考虑位移边界误差修正比不考虑位移边界误差修正的计算精度高；随着分步加载数目的增多，累积计算误差逐渐增大，对不考虑位移边界误差修正的误差累积效应更加明显，对考虑位移边界误差修正的影响不大；随着计算网格的加密，计算误差逐渐减小，计算精度逐渐提高。