多尺度六面体单元网格连接方法研究

方锡武, 邓正平, 蒋麒麟

(南京工程学院 机械工程学院,南京 211167)

1 引 言

工程领域经常涉及局部-全局分析的多物理场和多尺度问题[1-3]，这类问题通常分析区域大且分析对象多。当采用有限元方法分析时，其分析结果的精度与离散单元的精细程度直接相关，如果整个分析区域均采用精细网格，则解决问题所需的计算量巨大。故通常在感兴趣的重要区域用精细网格，而不重要区域用大尺度粗糙网格，但这样在粗糙网格和精细网格之间容易产生结点不匹配的单元或扭曲单元，导致过渡区域场量变化不连续或发生跳跃；另外，有限元分析不同物体对象之间发生的接触碰撞等物理问题时，由于不同对象的网格单元结点在接触界面很难匹配一致，也会导致界面两侧单元之间场量变化的不连续或跳跃现象。如何保证场量在不同尺度单元网格之间变化的连续性、一致性和各向同性是一个值得研究的课题。

目前解决此类问题的方法主要有两类，一类是直接法，即对两种尺度的网格各自采用有限元方法，而在网格的接触界面构建过渡模型方程或施加额外的约束条件方程。如Broughton等[4-7]提出的拉格朗日乘子法，Flemisch[8]采用的motar有限元方法，这些方法的本质是在界面施加一个弱连续条件，要求界面的场量变化向量正交一个合适的拉格朗日乘子空间，从变分意义上加强界面连续条件来求解，此方法增加了系统额外的自由度，需要修改非正定系统矩阵。Loehnert等[9]用正投影法获取接触界面另一侧网格单元的匹配信息，沿着界面用高斯积分来连接不同尺度的网格，此方法计算精度不高。Thevena等[10]采用近邻结点插值法获取两网格结点之间的过渡信息，将尺度不同的网格连接起来，这个方法简单且易执行，但如果在界面附近，一个网格结点远离另一网格结点，则误差不收敛。另一类是间接法，即将分析区域重建网格或特殊类型单元的网格，通过新的网格单元将非匹配结点网格转换成匹配网格。Houzeaux等[11]提出在两网格接触界面区域构建新的网格单元(即扩展单元)，将两种不同网格的单元结点连接起来，该方法需要在扩展单元结点施加隐含的边界条件来求解。Sohn等[12]采用多面体单元来完成结点不匹配的六面体单元网格之间的耦合作用，将六面体单元与其表面上新增加的额外结点一起考虑，重新划分成一系列多面体单元，构建多面体单元结点形函数，采用基于单元的光滑有限元方法实现网格连接。Sohn等[3,13,14]通过引入变结点单元(VNES)，直接增加新结点于接触界面单元的边上或者面上，原界面单元由变结点单元替换，基于移动最小二乘法或者点插值方法得到新增加结点的场量信息，使场量能在结点不匹配网格区域正确传递。此方法简单，求解容易，连接效果良好，但在构建形函数时要求多结点在单元边或面上呈规律分布，而实际应用中不同网格的单元结点很少呈规律分布，另外，移动最小二乘法构建的形函数受权函数取值半径的影响，方法难以适用于形状复杂、尺度较大的分析区域。

本文方法相较于上述两种方法，不需要建立额外的约束条件和过渡网格，只需在六面体单元等参插值方法基础上，构建额外结点形函数和修正本身结点形函数，形函数将结点场值的影响域限制在基点和三个方向近邻结点所确定的六面体区域内，保证了传统六面体单元有额外结点情形下的等参插值计算，从而实现两接触网格结点属性在接触界面的无缝连接，保证接触界面场值变化的连续性、一致性和各向同性。通过两个不同尺度的六面体单元的结点属性在接触界面的过渡实例和圆柱滚子轴承滚子与外圈接触的多尺度分析实例，验证了方法的正确性和有效性。

2 问题分析

有限元方法分析三维问题时，一般将研究区域划分为四面体或者六面体单元(本文的研究对象)网格，为了保证每个单元内场值变化的连续性，构建单元结点形函数，单元上任一点的场值由单元结点的场值插值计算得到，由于网格相邻单元结点的共有性和匹配性，保证了整个分析区域场值变化是连续性。但对于多尺度问题，如图1(a)所示，尺度大的网格A和尺度小的网格B接触，在其接触界面，由于网格单元尺度的不同和网格单元相对位置的不确定性，导致单元结点不能匹配；A网格上、左和前的一个单元在与B网格接触时，如图1(b)所示，增加了额外的B网格单元的结点1和2，同时由于A和B网格单元边相交而增加了交点结点3，4和5；同样，B网格下、左和前的一个单元，也增加了两网格单元边的交点结点4和5，如图1(c)所示。分析不同尺度的两个有限元网格的接触界面，其两单元面的四边形之间的相对位置有四种可能，① 一网格单元面完全在另一网格单元面内，如图1(d)所示； ② 一网格单元面有三个结点落在另一网格单元面内部，如图1(e)所示； ③ 一网格单元面有两个结点落在另一网格单元面内部，如图1(f)所示； ④ 一网格单元面有一个结点落在另一网格单元面内部，如图1(g)所示。上述情况均会给相互接触的六面体单元增加额外结点，可将六面体单元存在额外结点的情况归为四类， ① 单元一面存在多个额外结点； ② 单元相对面存在多结点； ③ 单元相邻面存在多结点； ④ 上述三种情况 共存。

图1 不同尺度单元网格接触分析

有限元分析时，如果对接触界面两侧的多结点单元仍采用传统的等参插值方法，则在接触界面的两侧单元会发生场值变化不连续的现象，如图2所示，出现此问题的根源是传统等参插值方法忽略了两网格单元接触时相互增加的额外结点带来的属性影响。如何处理额外结点对本单元带来的影响，以保证场值在接触界面的连续性，本文在传统等参插值方法基础上，通过构建额外新增加结点的形函数和修正原六面体单元八个结点的形函数，实现了场值在单元界面的连续传递，保证了整个分析区域场值分布的协调一致性。

图2 传统等参插值方法分析不同尺度单元网格

3 传统六面体单元等参插值方法分析

(1)

(2)

此等参插值方法满足结点场值在六面体单元内部和边界上的连续性。

图3 传统等参插值方法结点形函数分析

分析立方体单元八个结点形函数的构成(以 图3(b)结点1的形函数N1为例，其他结点类似)，可表示为

(3)

(4)

4 多结点六面体单元结点形函数构建

4.1 单元一个面上存在额外结点形函数构建方法

当一个单元面上存在额外结点，额外结点的分布有五种可能,① 单元面一边存在额外结点； ② 单元面相对边存在额外结点； ③ 单元面相邻边存在额外结点； ④ 单元面内部存在额外结点； ⑤ 上述四种情况的混合。按照第3节三维等参插值向二维和一维转化的思路，问题的焦点集中在四边形多结点插值方法上， 因为该单元所有结点的形函数沿额外结点所在单元面与其相对单元面方向的乘积因子并没有改变。对于四边形单元多结点插值方法，文献[16]针对上述的第①～③三种情况已有详细论述，第⑤种情况实质是第④种情况有额外结点在单元边的特例，因此本文主要论述第④种情况的处理方法。

(5)

j′和k′的结点场值由5和6结点场值线性插值得到，j″和k″的场值由7和8结点的场值线性插值得到。此方法的实质是将单元面5-6-7-8分解为5-j′-j″- 8,j′-k′-k″-j″和k′- 6-7-k″三个子四边形，这些子四边形中，5-j′-j″- 8的边j′-j″有一个额外结点j；j′-k′-k″-j的边j′-j″和k′-k″有额外结点j和k；k′- 6- 7-k″的边k′-k″有额外结点k。

(6)

(7)

(8)

(9)

式中i=1,2，…，7,8,j,j′,j″,k,k′,k″,且m≠i。

图5 单元面5-6-7-8上各结点形函数等值线云图和内部等值面图

图6 另一种处理方法及两种方法的单元等值线云图比较

4.2 单元相对面上存在额外结点形函数构建方法

4.3 单元相邻面上存在额外结点形函数构建方法

(10)

单元面m1-m2-m3-m4上四个结点与m点的ζ坐标均相同，形函数的ζ向乘积因子也相同，如m1结点ζ向有效取值线段为m1- 5和m1- 1，则其形函数ζ向乘积因子为

(11)

单元面1-2-3-4上所有结点的ζ向乘积因子相同，修正为(以结点1为例)

(12)

图7(b～e)为几个结点和虚结点的等值线云图和等值面图，限于篇幅，没有列出所有结点的等值线云图。

图7 单元相邻面存在多结点及其结点形函数等值线云图和内部等值面图

5 误差分析及消除

图8 误差分析及消除

6 实例验证

6.1 两不同尺度单元连接验证

图9(a)为不同尺度的相互接触的A和B两个六面体单元，其在单元面5-6-7-8和9-10-11-12之间相互接触，单元面5-6-7-8增加了额外结点 9～12，如图9(b)所示，按4.2节方法，在单元面5-6-7-8上可构建虚结点12′,12″,9′,9″,11′,11″,10′和10″；在单元B的9-10-11-12单元面，由于A单元面上9′-9″和11′-11″与其单元边有交点，则增加了9‴和11‴两个虚结点，如图9(c)所示，此两虚结点既属于单元A的单元面,也属于单元B的单元面，其结点场值由B单元面上的结点12,11和9，10的场值线性插值得到。在各自单元分别建立这些结点的形函数，将两个单元的结点赋予场值，按照本文方法，其表面等值线云图如图9(d)所示，与图2(a)的传统等参插值方法比较，其边界过渡处协调一致性良好。

图9 两不同尺度的接触单元

6.2 圆柱滚子轴承接触问题多尺度分析

图10(a)是圆柱滚子轴承模型，当滚子运转到最下端位置时，与内外圈之间的接触应力高度集中，直接决定轴承的疲劳寿命和轴承极限转速。对其进行有限元分析时，由于滚子、内圈和外圈是不同的零件，有限元网格尺度和结点的匹配性很难协调一致，运用本文方法则可以解决此问题。图10(b)是滚子与外圈接触的简化模型。外圈局部简化为底面固定的平板，滚子简化为半圆柱，半圆柱上矩形表面有均布压力载荷作用。模型的平板长宽高为12 mm，8 mm和2 mm，半圆柱半径为 3 mm，高度为12 mm；接触采用Hertz理论模型，滚子与平板之间的摩擦系数为0.3，杨氏模量E=2.06×1011Pa，泊松比ν=0.3。图10(c)是不同尺度两网格有限元分析模型，分别在滚子上表面施加不同载荷，图10(d)是第4组载荷作用下多尺度分析接触应力分布云图。不同载荷下的最大接触应力与理论公式(13)[17]计算结果比较列入表1，本文方法所得分析结果与理论公式计算结果误差在1%～3%，两者结果基本一致(注：式中P为均布载荷(N/m2)，R为半圆柱半径(m)，ν为泊松比，E为杨氏模量(N/m2)，σH为最大接触应力(N/m2))。

图10 圆柱滚子轴承多尺度分析

(13)

表1 最大接触应力的仿真值与理论计算值比较

7 结 论

(1) 以传统等参插值方法为基础，研究了多结点六面体单元结点形函数构建方法，形函数具备在整体坐标下收敛的完备性条件，实现了结点属性在不同尺度的有限元网格接触界面的连续光滑传递，保证了整个分析区域场值分布的协调一致性。

(2) 方法应用于圆柱滚子轴承滚子与外圈之间的接触分析，通过与理论计算结果比较，误差在5%之内，验证了方法正确性。

(3) 方法的不足之处是局限于8结点六面体单元等参插值方法和线弹性分析模型，对20结点的六面体单元则无能为力；同时对某些额外结点的处理需要增加虚结点，增大了问题的复杂性，这些问题需要进一步研究和完善。