高次光滑边界单元的构建与性能

田 钰, 和东宏, 马 杭

(上海大学力学与工程科学学院, 上海 200444)

固体内部的夹杂物(如粒子增强复合材料中的粒子)、第二相以及孔洞等对材料的性能有着重要的影响.对于这类材料的数值模拟, 为了精确地描述粒子的几何形状, 有限元需要采用大量的体积单元, 边界元只需在粒子边界, 即粒子与基体间的界面上剖分单元.然而, 为了能够精确地模拟粒子的曲面边界, 即使采用二次边界单元, 仍然需要较多数量的单元, 这就使得解题的规模变得十分庞大, 给问题的求解带来困难.近年来提出的等几何单元[1-3]能够准确地模拟曲线和曲面边界, 搭建了工程CAD 数据与科学计算之间的桥梁.等几何单元的不足之处是需要配置控制点来规定对象的几何形状, 与定义物理量的配位点位置并不相同.Gao 等[4]提出了高次等参数闭合单元来模拟二维圆和三维球面, 使得采用具有少量节点的一个单元来离散一个粒子成为可能.闭合单元的不足之处是存在端点或者端线, 在这些位置上, 形函数及其模拟的场变量的导数存在跳跃, 影响了数值模拟的精度.

本工作首先基于Lagrange 插值多项式, 将其写成单项式之和并进一步写成嵌套乘积的形式, 利用计算机程序自动计算嵌套乘积的系数, 避免了高次单元形函数冗长推导的麻烦; 然后,在已有闭合单元的基础上, 利用粒子几何形状的特性构建了椭圆和椭球高次光滑边界单元; 最后, 通过若干算例, 验证了高次光滑边界单元的精度与效率.

1 高次光滑边界单元

1.1 形函数系数

边界单元形函数的形成通常基于Lagrange 插值多项式.边界上的几何量或场变量可用节点上的数值近似表达为

式中:fk=f(ξk)为节点上的几何量或场变量,ξk(k=1,2,··· ,N+1)为参数空间的节点局部坐标,N+1 为节点数目;(ξ)为N阶Lagrange 多项式, 即

由于直接进行Lagrange 多项式计算的效率并不高, 通常都是将其展开进行计算, 这对于低次单元来说并无问题.而对于高次单元来说, 由于展开式的系数不仅取决于多项式的阶数, 也与节点的分布位置有关, 使得多项式的展开成为一件冗长的工作.为了避免这一问题, 在式(2)中当j ̸=k时将下标重新计数, 并记所有的ξj为(j=1,2,··· ,N), 展开Lagrange 多项式为单项式之和的形式, 即

式中:dk为式(2)中的分母, 对于给定的节点数目和节点位置,dk均为常数.式(3)中的系数可利用计算机程序自动加以计算, 其具体表达式为

这些系数可以预先算好备用.将式(3)进一步写成下列嵌套乘积的形式:

因为嵌套乘积的形式具有最高的计算效率[5], 容易写出多项式的导函数, 同样可通过计算机程序自动加以计算:

这样就避免了不同节点数目和不同节点分布时的高次单元形函数的一一推导.

1.2 光滑椭圆单元

线单元的起点和终点都是单元的端点, 端点上布置了节点的单元是协调单元.文献[4]所构建的闭合单元是协调单元, 位于单元起点和终点的两个节点的位置重合, 形成了闭合单元.由于闭合单元还存在一个位置重合的端点, 因此端点效应依然存在.为了避免端点效应, 以节点数N=8 为例构建的光滑椭圆单元如图1 所示.利用椭圆的周期性, 在不增加单元实际节点数目的条件下, 增加插值节点的数目, 单元上有3 个相邻节点被重复使用了2 次, 即节点1(9),2(10), 3(11)(见图1(a)).插值节点的数目增加以后, 位于单元起点的节点1(9)与单元终点的节点3(11)同时也位于单元的内部, 即节点1, 3 分别与节点9, 11 重合.在对应的参数空间中, 位置重合的节点用空心圆“◦”表示, 如图1(b)所示.由于插值节点数目增加了2 个, 插值多项式的阶数也随之提高了二阶, 即

图1 光滑椭圆单元及参数空间(N =8)Fig.1 Smooth elliptical element with the parametric space (N =8)

光滑椭圆单元的形函数定义如下:在光滑单元中, 位于单元起点和终点的节点同时也位于单元的内部, 从而能够去除端点效应,保证几何量或场变量的光滑性, 即导数的连续性.插值多项式阶数的提高以及端点效应的消除使得光滑单元的插值精度有所提高.采用3 个重合节点的目的是保持离散后的单元关于椭圆短半轴的对称性(见图1(a)).光滑椭圆单元的积分区间仍为[−1,+1], 与常规的闭合单元相同.

1.3 光滑椭球单元

节点数N= 14 的光滑椭球单元如图2 所示, 其中经线和纬线构成椭球表面的曲面坐标,节点13 和14 为极点,x3为旋转对称轴.参数空间如图2(b)所示, 其中局部坐标ξ1和ξ2分别对应于经线和纬线, 极点在参数平面中表达为双实线.分别记N1,N2为经线和纬线上节点的数目, 由于存在2 个极点, 则N1−2,N2分别是椭球单元上纬线和经线的数目, 它们与节点总数N的关系为

图2 光滑椭球单元及参数空间(N =14)Fig.2 Smooth ellipsoidal element with the parametric space (N =14)

用k1,k2对经线和纬线上的节点进行局部编号, 则单元整体编号m与局部编号的关系为

两个极点的整体编号分别为m=N −1 和m=N.在图2(b)所示平面中, 正方形阴影部分即文献[4]所构建闭合单元的参数平面.注意到, 经线都是开放的半圆或半椭圆, 两个极点是所有经线的共同端点, 当ξ2=±1 所对应的经线就是闭合单元的端线.端点和端线的存在破坏了几何量或场变量导数的连续性, 降低了插值精度.解决办法是利用椭球的周期性, 在不增加单元实际节点数目的条件下, 沿纬线和经线两个方向上分别增加插值节点的数目, 构建光滑椭球单元.由于纬线都是圆, 节点增加的方式以及插值多项式的形成与二维光滑椭圆单元完全相同, 纬线方向的形函数(ξ)仍可用式(7)表示(分别用k2和N2替换式中的k和N).以赤道为例(如图3(a)), 整体编号为8, 5, 6 的节点是沿纬线方向重复使用两次的节点, 括号中的数字为节点的局部编号, 其中重复使用的节点都有两个局部编号.在图2(b)所示的参数平面中, 沿纬线方向位置重合的节点用空心圆“◦”表示.

图3(b)为互为镜像的两条经线及其节点的整体编号, 括号中的数字则为当前经线上节点的局部编号, 其中当前经线和镜像经线分别用实线和虚线表示, 点划线为椭球的旋转对称轴.对于经线方向, 在经线两端越过极点各增加s个插值节点(本工作取s=2), 使得两个极点都位于单元的内部, 不再是经线的端点.沿经线方向增加的节点用棱形符号“◇”表示(见图3(b)),其插值多项式为

图3 赤道与经线上节点的整体与局部编号Fig.3 Global and local numberings for nodes along the equator and the meridians

当前经线所增加的插值节点的位置越过了极点, 从而位于虚线表示的镜像经线之上, 括号中出现的0 甚至负数是为了保持节点局部编号的连续性.由于利用了椭球的周期性, 光滑椭球单元上的一些节点被重复地使用, 重复使用节点的整体编号添加“′”以示区别(见图2(b)).记经线方向的形函数为

光滑椭球单元的形函数由纬线方向和经线方向形函数的乘积来表达, 即

式中:k2=1,2,··· ,N2;m由式(10)定义;M(k2) 表示镜像函数, 具有

由此可知, 光滑椭球单元的经线数目或纬线的节点数目必须采用偶数.对应于两个极点的形函数定义为

需要注意的是: ①单元极点处的形函数仅仅由经线方向的形函数来定义, 所以位于两个极点处的曲面的法线方向是不确定的; ②经线方向增加的插值节点越过了极点位于其镜像位置.图2(b)所示参数平面中棱形符号“◇”所在区域上的法线方向与积分区间的法线方向是相反的, 但对单元积分并没有任何影响, 因为光滑椭球单元沿纬线和经线两个方向上的积分区间均为[−1, +1], 即正方形阴影部分, 该积分区间与常规的闭合单元[4]完全一致.事实上, 由于本工作构建的光滑椭球单元利用了椭球的周期性和对称性, 在纬线和经线两个方向上增加了重复使用的插值节点, 提高了插值多项式的阶数, 提高了节点的利用率, 消除了端点效应, 从而能够较大地提高插值的精度, 而单元的实际节点数目并没有变化.

2 数值算例

2.1 光滑椭圆单元

椭圆的几何参数如图4 所示, 其中a和b分别为椭圆的长半轴和短半轴.分别采用二次单元、闭合单元、光滑单元进行离散.用二次单元离散时, 节点总数与单元数之比为2∶1.采用高次单元进行离散, 只需用一个单元.用3 种单元计算的圆面积(b/a= 1)相对误差与节点总数的关系如图5(a)所示.由图5(a)可以看出, 两种高次单元的计算精度都远高于二次单元, 而且随着节点总数的增加, 高次单元相对误差下降的速率也远高于二次单元.由于本工作构建的光滑椭圆单元提高了插值多项式的阶数, 消除了端点效应, 光滑单元的计算精度比常规闭合单元大约提高了一个数量级.

图4 椭圆的几何参数Fig.4 Geometrical parameters of an ellipse

定义位置误差为离散椭圆边线与理想椭圆边线的距离.采用两种8 节点高次单元计算沿椭圆周边的相对位置误差, 结果如图5(b)所示.由图5(b)同样可以看出, 光滑单元的计算精度比常规闭合单元大约提高了一个数量级.从整体的效果来看, 无论对于闭合单元还是光滑单元, 单元中部(ξ=0 附近)的精度优于靠近两端的精度.

图5 面积与位置相对误差的对比Fig.5 Comparisons of relative errors of areas and positions

采用8 节点光滑单元分别计算无限域中圆孔和椭圆孔在远场x2方向单位载荷(σ0=1)作用下的场变量, 即位移和应力.图6(a)为第一象限边界上的位移与应力计算值与理论解的比较, 图6(b)为靠近边界(沿图4(a)第一象限中的虚线,d/a=0.01)处的位移与应力计算值与理论值的比较.当计算点靠近边界时, 产生的近奇异积分采用距离变换的数值积分技术加以解决[6-7].图6 表明本工作构建的光滑椭圆单元的计算值与理论解吻合良好.

图6 无限域中圆孔和椭圆孔的场变量Fig.6 Field variables on boundary of a circle and close to boundary of an ellipse

采用光滑椭圆单元在不同条件下对无限域中椭圆孔前方x2方向的应力分量的分布进行了计算, 结果如图7 所示.由图7 可知, 采用光滑椭圆单元进行离散, 当椭圆的形状参数b/a降低时, 应适当增加保持精度所需要的节点数目.总体而言, 对于光滑单元, 采用数量不多的节点即可获得满意的精度, 例如圆孔(b/a=1)仅用了4 个节点, 这是其他类型的单元所做不到的.

图7 椭圆孔前方的应力分布Fig.7 Stress distributions in domain ahead of the elliptical holes

2.2 光滑椭球单元及应用

首先比较传统闭合单元与本工作构建的光滑单元几何参数的精度.图8 为圆球表面积与体积的计算相对误差随着高次单元节点总数的变化情况, 其中采用96 个单元290 个节点的8 节点二次单元[8]进行计算得到的表面积与体积各表达为一个点, 给出的水平虚线是为了便于误差的比较.由图8 可知, 两种高次单元的计算精度都远远高于二次单元.由于提高了插值多项式的阶数以及消除了端点效应, 光滑单元的计算精度比闭合单元提高了一到两个数量级, 而且随着节点总数的增加, 光滑椭球单元相对误差下降的速率也高于传统的闭合单元.

图8 圆球表面积与体积的相对误差对比Fig.8 Comparisons of relative errors of surface area and volume of ball

采用26 个节点对两种高次单元计算的圆球半径相对误差与曲面法线误差沿1/2 赤道(ξ1=1)和1/2 经线(ξ2=1)的分布进行计算分析, 结果如图9 所示, 其中1/2 赤道和1/2 经线分别对应于纬线与经线方向上半个单元的范围(图2(b)中的阴影区域).曲面法线误差为

式中:n和n0分别为曲面单位法线向量的计算值与理论值.由图9 可知, 半径相对误差与曲面法线误差的极大值分别出现在节点之间区域与节点位置上, 插值多项式阶数的提高和端点效应的消除, 使得光滑单元的计算精度比闭合单元提高了一到两个数量级.从整体的效果来看,无论对于纬线还是经线方向, 单元中部(ξ1=0 与ξ2=0 附近)的精度优于靠近两端的精度, 与二维的情况类似.由于纬线都是封闭的圆, 而经线都是开放的半圆或半椭圆, 沿纬线的拟合效果(见图9(b))整体上要由优于沿经线的拟合效果(见图9(a)).

图9 圆球沿1/2 赤道和1/2 经线的半径相对误差与曲面法线误差Fig.9 Relative errors of radius and errors of surface normal along 1/2 equator and 1/2 maridian of ball

在应用方面, 以解析解[9]为参照基准, 对旋转椭球的Eshelby张量Sijkl进行数值计算和比较.Eshelby张量的数值计算采用下式[8]进行:

式中:Γ为椭球边界;为面力基本解导出函数;δij为kronecke 符号;xl为场点与源点距离在坐标轴上的投影;µ,v分别为材料剪切模量与泊松比.令x3为椭球的旋转对称轴(见图2(a)),c为对称轴的半轴长,a为椭球x1和x2方向的半轴长.采用52 个节点的光滑椭球单元, 计算了旋转椭球的7 个不同的Eshelby 张量随c/a的变化规律, 计算结果与解析解的对比如图10(a)所示.可以看出, 二者吻合良好.同时采用96 单元290 节点的二次单元对同一问题进行计算, 比较二次单元与光滑单元计算结果相对于解析解的绝对误差, 结果如图10(b)所示.可以看出, 光滑单元的计算精度比二次单元提高了2∼3 个数量级.

图10 椭球Eshelby 张量和绝对误差的比较Fig.10 Comparisons of Eshelby tensors and absolute errors of ellipsoids

在数值计算时, 每个二次单元采用8×8 点的高斯积分, 总的高斯积分点数为8×8×96=6 144.对于光滑单元, 当c/a接近于1 时采用高斯积分点数为20×20 = 400; 当c/a较小或较大(分别对应于扁椭圆和长椭圆)时, 需要将积分域划分为4 个极坐标描述的三角形子域并采用近奇异积分技术[6-7]进行数值积分, 总的高斯积分点数为20×8×4=640, 其中20 和8 分别为子域的径向和角度方向的积分点数.两种单元的高斯积分点总数相差超过9 倍, 说明本工作构建的光滑单元不仅有很高的计算精度, 也有很高的计算效率.可以预期, 将光滑单元用于粒子增强材料的数值模拟, 能够提高该类材料性能模拟计算的精度和效率.

3 结束语

本工作通过将Lagrange 插值多项式转化为嵌套积的形式, 实现了高次边界单元形函数系数的计算机自动生成.在已有闭合单元的基础上, 充分利用粒子几何形状的特性, 在不增加实际节点数目的条件下扩大了参与插值的节点数目, 构建了高次光滑边界单元, 提高了插值多项式的阶数, 消除了端点效应.数值算例表明, 与传统的二次单元和闭合单元相比, 高次光滑边界单元能够较大地提高椭圆和椭球粒子数值模拟的计算精度与效率, 适用于粒子增强材料的数值模拟.