方錫武, 鄧正平, 蔣麒麟

(南京工程學(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ù)，采用基于單元的光滑有限元方法實(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ī)律分布，而實(shí)際應(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)情形下的等參插值計算，從而實(shí)現(xiàn)兩接觸網(wǎng)格結(jié)點(diǎn)屬性在接觸界面的無縫連接，保證接觸界面場值變化的連續(xù)性、一致性和各向同性。通過兩個不同尺度的六面體單元的結(jié)點(diǎn)屬性在接觸界面的過渡實(shí)例和圓柱滾子軸承滾子與外圈接觸的多尺度分析實(shí)例，驗(yà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ù)，實(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)插值方法上， 因?yàn)樵搯卧薪Y(jié)點(diǎn)的形函數(shù)沿額外結(jié)點(diǎn)所在單元面與其相對單元面方向的乘積因子并沒有改變。對于四邊形單元多結(jié)點(diǎn)插值方法，文獻(xiàn)[16]針對上述的第①～③三種情況已有詳細(xì)論述，第⑤種情況實(shí)質(zhì)是第④種情況有額外結(jié)點(diǎn)在單元邊的特例，因此本文主要論述第④種情況的處理方法。

(5)

j′和k′的結(jié)點(diǎn)場值由5和6結(jié)點(diǎn)場值線性插值得到，j″和k″的場值由7和8結(jié)點(diǎn)的場值線性插值得到。此方法的實(shí)質(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 實(shí)例驗(yàn)證

6.1 兩不同尺度單元連接驗(yàn)證

圖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)下收斂的完備性條件，實(shí)現(xiàn)了結(jié)點(diǎn)屬性在不同尺度的有限元網(wǎng)格接觸界面的連續(xù)光滑傳遞，保證了整個分析區(qū)域場值分布的協(xié)調(diào)一致性。

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

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