多尺度六面體單元網(wǎng)格連接方法研究

方錫武, 鄧正平, 蔣麒麟

(南京工程學(xué)院 機(jī)械工程學(xué)院,南京 211167)

1 引 言

工程領(lǐng)域經(jīng)常涉及局部-全局分析的多物理場和多尺度問題[1-3]，這類問題通常分析區(qū)域大且分析對象多。當(dāng)采用有限元方法分析時，其分析結(jié)果的精度與離散單元的精細(xì)程度直接相關(guān)，如果整個分析區(qū)域均采用精細(xì)網(wǎng)格，則解決問題所需的計算量巨大。故通常在感興趣的重要區(qū)域用精細(xì)網(wǎng)格，而不重要區(qū)域用大尺度粗糙網(wǎng)格，但這樣在粗糙網(wǎng)格和精細(xì)網(wǎng)格之間容易產(chǎn)生結(jié)點(diǎn)不匹配的單元或扭曲單元，導(dǎo)致過渡區(qū)域場量變化不連續(xù)或發(fā)生跳躍；另外，有限元分析不同物體對象之間發(fā)生的接觸碰撞等物理問題時，由于不同對象的網(wǎng)格單元結(jié)點(diǎn)在接觸界面很難匹配一致，也會導(dǎo)致界面兩側(cè)單元之間場量變化的不連續(xù)或跳躍現(xiàn)象。如何保證場量在不同尺度單元網(wǎng)格之間變化的連續(xù)性、一致性和各向同性是一個值得研究的課題。

目前解決此類問題的方法主要有兩類，一類是直接法，即對兩種尺度的網(wǎng)格各自采用有限元方法，而在網(wǎng)格的接觸界面構(gòu)建過渡模型方程或施加額外的約束條件方程。如Broughton等[4-7]提出的拉格朗日乘子法，F(xiàn)lemisch[8]采用的motar有限元方法，這些方法的本質(zhì)是在界面施加一個弱連續(xù)條件，要求界面的場量變化向量正交一個合適的拉格朗日乘子空間，從變分意義上加強(qiáng)界面連續(xù)條件來求解，此方法增加了系統(tǒng)額外的自由度，需要修改非正定系統(tǒng)矩陣。Loehnert等[9]用正投影法獲取接觸界面另一側(cè)網(wǎng)格單元的匹配信息，沿著界面用高斯積分來連接不同尺度的網(wǎng)格，此方法計算精度不高。Thevena等[10]采用近鄰結(jié)點(diǎn)插值法獲取兩網(wǎng)格結(jié)點(diǎn)之間的過渡信息，將尺度不同的網(wǎng)格連接起來，這個方法簡單且易執(zhí)行，但如果在界面附近，一個網(wǎng)格結(jié)點(diǎn)遠(yuǎn)離另一網(wǎng)格結(jié)點(diǎn)，則誤差不收斂。另一類是間接法，即將分析區(qū)域重建網(wǎng)格或特殊類型單元的網(wǎng)格，通過新的網(wǎng)格單元將非匹配結(jié)點(diǎn)網(wǎng)格轉(zhuǎn)換成匹配網(wǎng)格。Houzeaux等[11]提出在兩網(wǎng)格接觸界面區(qū)域構(gòu)建新的網(wǎng)格單元(即擴(kuò)展單元)，將兩種不同網(wǎng)格的單元結(jié)點(diǎn)連接起來，該方法需要在擴(kuò)展單元結(jié)點(diǎn)施加隱含的邊界條件來求解。Sohn等[12]采用多面體單元來完成結(jié)點(diǎn)不匹配的六面體單元網(wǎng)格之間的耦合作用，將六面體單元與其表面上新增加的額外結(jié)點(diǎn)一起考慮，重新劃分成一系列多面體單元，構(gòu)建多面體單元結(jié)點(diǎn)形函數(shù)，采用基于單元的光滑有限元方法實現(xiàn)網(wǎng)格連接。Sohn等[3,13,14]通過引入變結(jié)點(diǎn)單元(VNES)，直接增加新結(jié)點(diǎn)于接觸界面單元的邊上或者面上，原界面單元由變結(jié)點(diǎn)單元替換，基于移動最小二乘法或者點(diǎn)插值方法得到新增加結(jié)點(diǎn)的場量信息，使場量能在結(jié)點(diǎn)不匹配網(wǎng)格區(qū)域正確傳遞。此方法簡單，求解容易，連接效果良好，但在構(gòu)建形函數(shù)時要求多結(jié)點(diǎn)在單元邊或面上呈規(guī)律分布，而實際應(yīng)用中不同網(wǎng)格的單元結(jié)點(diǎn)很少呈規(guī)律分布，另外，移動最小二乘法構(gòu)建的形函數(shù)受權(quán)函數(shù)取值半徑的影響，方法難以適用于形狀復(fù)雜、尺度較大的分析區(qū)域。

本文方法相較于上述兩種方法，不需要建立額外的約束條件和過渡網(wǎng)格，只需在六面體單元等參插值方法基礎(chǔ)上，構(gòu)建額外結(jié)點(diǎn)形函數(shù)和修正本身結(jié)點(diǎn)形函數(shù)，形函數(shù)將結(jié)點(diǎn)場值的影響域限制在基點(diǎn)和三個方向近鄰結(jié)點(diǎn)所確定的六面體區(qū)域內(nèi)，保證了傳統(tǒng)六面體單元有額外結(jié)點(diǎn)情形下的等參插值計算，從而實現(xiàn)兩接觸網(wǎng)格結(jié)點(diǎn)屬性在接觸界面的無縫連接，保證接觸界面場值變化的連續(xù)性、一致性和各向同性。通過兩個不同尺度的六面體單元的結(jié)點(diǎn)屬性在接觸界面的過渡實例和圓柱滾子軸承滾子與外圈接觸的多尺度分析實例，驗證了方法的正確性和有效性。

2 問題分析

有限元方法分析三維問題時，一般將研究區(qū)域劃分為四面體或者六面體單元(本文的研究對象)網(wǎng)格，為了保證每個單元內(nèi)場值變化的連續(xù)性，構(gòu)建單元結(jié)點(diǎn)形函數(shù)，單元上任一點(diǎn)的場值由單元結(jié)點(diǎn)的場值插值計算得到，由于網(wǎng)格相鄰單元結(jié)點(diǎn)的共有性和匹配性，保證了整個分析區(qū)域場值變化是連續(xù)性。但對于多尺度問題，如圖1(a)所示，尺度大的網(wǎng)格A和尺度小的網(wǎng)格B接觸，在其接觸界面，由于網(wǎng)格單元尺度的不同和網(wǎng)格單元相對位置的不確定性，導(dǎo)致單元結(jié)點(diǎn)不能匹配；A網(wǎng)格上、左和前的一個單元在與B網(wǎng)格接觸時，如圖1(b)所示，增加了額外的B網(wǎng)格單元的結(jié)點(diǎn)1和2，同時由于A和B網(wǎng)格單元邊相交而增加了交點(diǎn)結(jié)點(diǎn)3，4和5；同樣，B網(wǎng)格下、左和前的一個單元，也增加了兩網(wǎng)格單元邊的交點(diǎn)結(jié)點(diǎn)4和5，如圖1(c)所示。分析不同尺度的兩個有限元網(wǎng)格的接觸界面，其兩單元面的四邊形之間的相對位置有四種可能，① 一網(wǎng)格單元面完全在另一網(wǎng)格單元面內(nèi)，如圖1(d)所示； ② 一網(wǎng)格單元面有三個結(jié)點(diǎn)落在另一網(wǎng)格單元面內(nèi)部，如圖1(e)所示； ③ 一網(wǎng)格單元面有兩個結(jié)點(diǎn)落在另一網(wǎng)格單元面內(nèi)部，如圖1(f)所示； ④ 一網(wǎng)格單元面有一個結(jié)點(diǎn)落在另一網(wǎng)格單元面內(nèi)部，如圖1(g)所示。上述情況均會給相互接觸的六面體單元增加額外結(jié)點(diǎn)，可將六面體單元存在額外結(jié)點(diǎn)的情況歸為四類， ① 單元一面存在多個額外結(jié)點(diǎn)； ② 單元相對面存在多結(jié)點(diǎn)； ③ 單元相鄰面存在多結(jié)點(diǎn)； ④ 上述三種情況 共存。

圖1 不同尺度單元網(wǎng)格接觸分析

有限元分析時，如果對接觸界面兩側(cè)的多結(jié)點(diǎn)單元仍采用傳統(tǒng)的等參插值方法，則在接觸界面的兩側(cè)單元會發(fā)生場值變化不連續(xù)的現(xiàn)象，如圖2所示，出現(xiàn)此問題的根源是傳統(tǒng)等參插值方法忽略了兩網(wǎng)格單元接觸時相互增加的額外結(jié)點(diǎn)帶來的屬性影響。如何處理額外結(jié)點(diǎn)對本單元帶來的影響，以保證場值在接觸界面的連續(xù)性，本文在傳統(tǒng)等參插值方法基礎(chǔ)上，通過構(gòu)建額外新增加結(jié)點(diǎn)的形函數(shù)和修正原六面體單元八個結(jié)點(diǎn)的形函數(shù)，實現(xiàn)了場值在單元界面的連續(xù)傳遞，保證了整個分析區(qū)域場值分布的協(xié)調(diào)一致性。

圖2 傳統(tǒng)等參插值方法分析不同尺度單元網(wǎng)格

3 傳統(tǒng)六面體單元等參插值方法分析

(1)

(2)

此等參插值方法滿足結(jié)點(diǎn)場值在六面體單元內(nèi)部和邊界上的連續(xù)性。

圖3 傳統(tǒng)等參插值方法結(jié)點(diǎn)形函數(shù)分析

分析立方體單元八個結(jié)點(diǎn)形函數(shù)的構(gòu)成(以 圖3(b)結(jié)點(diǎn)1的形函數(shù)N1為例，其他結(jié)點(diǎn)類似)，可表示為

(3)

(4)

4 多結(jié)點(diǎn)六面體單元結(jié)點(diǎn)形函數(shù)構(gòu)建

4.1 單元一個面上存在額外結(jié)點(diǎn)形函數(shù)構(gòu)建方法

當(dāng)一個單元面上存在額外結(jié)點(diǎn)，額外結(jié)點(diǎn)的分布有五種可能,① 單元面一邊存在額外結(jié)點(diǎn)； ② 單元面相對邊存在額外結(jié)點(diǎn)； ③ 單元面相鄰邊存在額外結(jié)點(diǎn)； ④ 單元面內(nèi)部存在額外結(jié)點(diǎn)； ⑤ 上述四種情況的混合。按照第3節(jié)三維等參插值向二維和一維轉(zhuǎn)化的思路，問題的焦點(diǎn)集中在四邊形多結(jié)點(diǎn)插值方法上， 因為該單元所有結(jié)點(diǎn)的形函數(shù)沿額外結(jié)點(diǎn)所在單元面與其相對單元面方向的乘積因子并沒有改變。對于四邊形單元多結(jié)點(diǎn)插值方法，文獻(xiàn)[16]針對上述的第①～③三種情況已有詳細(xì)論述，第⑤種情況實質(zhì)是第④種情況有額外結(jié)點(diǎn)在單元邊的特例，因此本文主要論述第④種情況的處理方法。

(5)

j′和k′的結(jié)點(diǎn)場值由5和6結(jié)點(diǎn)場值線性插值得到，j″和k″的場值由7和8結(jié)點(diǎn)的場值線性插值得到。此方法的實質(zhì)是將單元面5-6-7-8分解為5-j′-j″- 8,j′-k′-k″-j″和k′- 6-7-k″三個子四邊形，這些子四邊形中，5-j′-j″- 8的邊j′-j″有一個額外結(jié)點(diǎn)j；j′-k′-k″-j的邊j′-j″和k′-k″有額外結(jié)點(diǎn)j和k；k′- 6- 7-k″的邊k′-k″有額外結(jié)點(diǎn)k。

(6)

(7)

(8)

(9)

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

圖5 單元面5-6-7-8上各結(jié)點(diǎn)形函數(shù)等值線云圖和內(nèi)部等值面圖

圖6 另一種處理方法及兩種方法的單元等值線云圖比較

4.2 單元相對面上存在額外結(jié)點(diǎn)形函數(shù)構(gòu)建方法

4.3 單元相鄰面上存在額外結(jié)點(diǎn)形函數(shù)構(gòu)建方法

(10)

單元面m1-m2-m3-m4上四個結(jié)點(diǎn)與m點(diǎn)的ζ坐標(biāo)均相同，形函數(shù)的ζ向乘積因子也相同，如m1結(jié)點(diǎn)ζ向有效取值線段為m1- 5和m1- 1，則其形函數(shù)ζ向乘積因子為

(11)

單元面1-2-3-4上所有結(jié)點(diǎn)的ζ向乘積因子相同，修正為(以結(jié)點(diǎn)1為例)

(12)

圖7(b～e)為幾個結(jié)點(diǎn)和虛結(jié)點(diǎn)的等值線云圖和等值面圖，限于篇幅，沒有列出所有結(jié)點(diǎn)的等值線云圖。

圖7 單元相鄰面存在多結(jié)點(diǎn)及其結(jié)點(diǎn)形函數(shù)等值線云圖和內(nèi)部等值面圖

5 誤差分析及消除

圖8 誤差分析及消除

6 實例驗證

6.1 兩不同尺度單元連接驗證

圖9(a)為不同尺度的相互接觸的A和B兩個六面體單元，其在單元面5-6-7-8和9-10-11-12之間相互接觸，單元面5-6-7-8增加了額外結(jié)點(diǎn) 9～12，如圖9(b)所示，按4.2節(jié)方法，在單元面5-6-7-8上可構(gòu)建虛結(jié)點(diǎn)12′,12″,9′,9″,11′,11″,10′和10″；在單元B的9-10-11-12單元面，由于A單元面上9′-9″和11′-11″與其單元邊有交點(diǎn)，則增加了9?和11?兩個虛結(jié)點(diǎn)，如圖9(c)所示，此兩虛結(jié)點(diǎn)既屬于單元A的單元面,也屬于單元B的單元面，其結(jié)點(diǎn)場值由B單元面上的結(jié)點(diǎn)12,11和9，10的場值線性插值得到。在各自單元分別建立這些結(jié)點(diǎn)的形函數(shù)，將兩個單元的結(jié)點(diǎn)賦予場值，按照本文方法，其表面等值線云圖如圖9(d)所示，與圖2(a)的傳統(tǒng)等參插值方法比較，其邊界過渡處協(xié)調(diào)一致性良好。

圖9 兩不同尺度的接觸單元

6.2 圓柱滾子軸承接觸問題多尺度分析

圖10(a)是圓柱滾子軸承模型，當(dāng)滾子運(yùn)轉(zhuǎn)到最下端位置時，與內(nèi)外圈之間的接觸應(yīng)力高度集中，直接決定軸承的疲勞壽命和軸承極限轉(zhuǎn)速。對其進(jìn)行有限元分析時，由于滾子、內(nèi)圈和外圈是不同的零件，有限元網(wǎng)格尺度和結(jié)點(diǎn)的匹配性很難協(xié)調(diào)一致，運(yùn)用本文方法則可以解決此問題。圖10(b)是滾子與外圈接觸的簡化模型。外圈局部簡化為底面固定的平板，滾子簡化為半圓柱，半圓柱上矩形表面有均布壓力載荷作用。模型的平板長寬高為12 mm，8 mm和2 mm，半圓柱半徑為 3 mm，高度為12 mm；接觸采用Hertz理論模型，滾子與平板之間的摩擦系數(shù)為0.3，楊氏模量E=2.06×1011Pa，泊松比ν=0.3。圖10(c)是不同尺度兩網(wǎng)格有限元分析模型，分別在滾子上表面施加不同載荷，圖10(d)是第4組載荷作用下多尺度分析接觸應(yīng)力分布云圖。不同載荷下的最大接觸應(yīng)力與理論公式(13)[17]計算結(jié)果比較列入表1，本文方法所得分析結(jié)果與理論公式計算結(jié)果誤差在1%～3%，兩者結(jié)果基本一致(注：式中P為均布載荷(N/m2)，R為半圓柱半徑(m)，ν為泊松比，E為楊氏模量(N/m2)，σH為最大接觸應(yīng)力(N/m2))。

圖10 圓柱滾子軸承多尺度分析

(13)

表1 最大接觸應(yīng)力的仿真值與理論計算值比較

7 結(jié) 論

(1) 以傳統(tǒng)等參插值方法為基礎(chǔ)，研究了多結(jié)點(diǎn)六面體單元結(jié)點(diǎn)形函數(shù)構(gòu)建方法，形函數(shù)具備在整體坐標(biāo)下收斂的完備性條件，實現(xiàn)了結(jié)點(diǎn)屬性在不同尺度的有限元網(wǎng)格接觸界面的連續(xù)光滑傳遞，保證了整個分析區(qū)域場值分布的協(xié)調(diào)一致性。

(2) 方法應(yīng)用于圓柱滾子軸承滾子與外圈之間的接觸分析，通過與理論計算結(jié)果比較，誤差在5%之內(nèi)，驗證了方法正確性。

(3) 方法的不足之處是局限于8結(jié)點(diǎn)六面體單元等參插值方法和線彈性分析模型，對20結(jié)點(diǎn)的六面體單元則無能為力；同時對某些額外結(jié)點(diǎn)的處理需要增加虛結(jié)點(diǎn)，增大了問題的復(fù)雜性，這些問題需要進(jìn)一步研究和完善。