田 鈺, 和東宏, 馬 杭

(上海大學(xué)力學(xué)與工程科學(xué)學(xué)院, 上海 200444)

固體內(nèi)部的夾雜物(如粒子增強(qiáng)復(fù)合材料中的粒子)、第二相以及孔洞等對材料的性能有著重要的影響.對于這類材料的數(shù)值模擬, 為了精確地描述粒子的幾何形狀, 有限元需要采用大量的體積單元, 邊界元只需在粒子邊界, 即粒子與基體間的界面上剖分單元.然而, 為了能夠精確地模擬粒子的曲面邊界, 即使采用二次邊界單元, 仍然需要較多數(shù)量的單元, 這就使得解題的規(guī)模變得十分龐大, 給問題的求解帶來困難.近年來提出的等幾何單元[1-3]能夠準(zhǔn)確地模擬曲線和曲面邊界, 搭建了工程CAD 數(shù)據(jù)與科學(xué)計算之間的橋梁.等幾何單元的不足之處是需要配置控制點來規(guī)定對象的幾何形狀, 與定義物理量的配位點位置并不相同.Gao 等[4]提出了高次等參數(shù)閉合單元來模擬二維圓和三維球面, 使得采用具有少量節(jié)點的一個單元來離散一個粒子成為可能.閉合單元的不足之處是存在端點或者端線, 在這些位置上, 形函數(shù)及其模擬的場變量的導(dǎo)數(shù)存在跳躍, 影響了數(shù)值模擬的精度.

本工作首先基于Lagrange 插值多項式, 將其寫成單項式之和并進(jìn)一步寫成嵌套乘積的形式, 利用計算機(jī)程序自動計算嵌套乘積的系數(shù), 避免了高次單元形函數(shù)冗長推導(dǎo)的麻煩; 然后,在已有閉合單元的基礎(chǔ)上, 利用粒子幾何形狀的特性構(gòu)建了橢圓和橢球高次光滑邊界單元; 最后, 通過若干算例, 驗證了高次光滑邊界單元的精度與效率.

1 高次光滑邊界單元

1.1 形函數(shù)系數(shù)

邊界單元形函數(shù)的形成通常基于Lagrange 插值多項式.邊界上的幾何量或場變量可用節(jié)點上的數(shù)值近似表達(dá)為

式中:fk=f(ξk)為節(jié)點上的幾何量或場變量,ξk(k=1,2,··· ,N+1)為參數(shù)空間的節(jié)點局部坐標(biāo),N+1 為節(jié)點數(shù)目;(ξ)為N階Lagrange 多項式, 即

由于直接進(jìn)行Lagrange 多項式計算的效率并不高, 通常都是將其展開進(jìn)行計算, 這對于低次單元來說并無問題.而對于高次單元來說, 由于展開式的系數(shù)不僅取決于多項式的階數(shù), 也與節(jié)點的分布位置有關(guān), 使得多項式的展開成為一件冗長的工作.為了避免這一問題, 在式(2)中當(dāng)j ?=k時將下標(biāo)重新計數(shù), 并記所有的ξj為(j=1,2,··· ,N), 展開Lagrange 多項式為單項式之和的形式, 即

式中:dk為式(2)中的分母, 對于給定的節(jié)點數(shù)目和節(jié)點位置,dk均為常數(shù).式(3)中的系數(shù)可利用計算機(jī)程序自動加以計算, 其具體表達(dá)式為

這些系數(shù)可以預(yù)先算好備用.將式(3)進(jìn)一步寫成下列嵌套乘積的形式:

因為嵌套乘積的形式具有最高的計算效率[5], 容易寫出多項式的導(dǎo)函數(shù), 同樣可通過計算機(jī)程序自動加以計算:

這樣就避免了不同節(jié)點數(shù)目和不同節(jié)點分布時的高次單元形函數(shù)的一一推導(dǎo).

1.2 光滑橢圓單元

線單元的起點和終點都是單元的端點, 端點上布置了節(jié)點的單元是協(xié)調(diào)單元.文獻(xiàn)[4]所構(gòu)建的閉合單元是協(xié)調(diào)單元, 位于單元起點和終點的兩個節(jié)點的位置重合, 形成了閉合單元.由于閉合單元還存在一個位置重合的端點, 因此端點效應(yīng)依然存在.為了避免端點效應(yīng), 以節(jié)點數(shù)N=8 為例構(gòu)建的光滑橢圓單元如圖1 所示.利用橢圓的周期性, 在不增加單元實際節(jié)點數(shù)目的條件下, 增加插值節(jié)點的數(shù)目, 單元上有3 個相鄰節(jié)點被重復(fù)使用了2 次, 即節(jié)點1(9),2(10), 3(11)(見圖1(a)).插值節(jié)點的數(shù)目增加以后, 位于單元起點的節(jié)點1(9)與單元終點的節(jié)點3(11)同時也位于單元的內(nèi)部, 即節(jié)點1, 3 分別與節(jié)點9, 11 重合.在對應(yīng)的參數(shù)空間中, 位置重合的節(jié)點用空心圓“?”表示, 如圖1(b)所示.由于插值節(jié)點數(shù)目增加了2 個, 插值多項式的階數(shù)也隨之提高了二階, 即

圖1 光滑橢圓單元及參數(shù)空間(N =8)Fig.1 Smooth elliptical element with the parametric space (N =8)

光滑橢圓單元的形函數(shù)定義如下:在光滑單元中, 位于單元起點和終點的節(jié)點同時也位于單元的內(nèi)部, 從而能夠去除端點效應(yīng),保證幾何量或場變量的光滑性, 即導(dǎo)數(shù)的連續(xù)性.插值多項式階數(shù)的提高以及端點效應(yīng)的消除使得光滑單元的插值精度有所提高.采用3 個重合節(jié)點的目的是保持離散后的單元關(guān)于橢圓短半軸的對稱性(見圖1(a)).光滑橢圓單元的積分區(qū)間仍為[?1,+1], 與常規(guī)的閉合單元相同.

1.3 光滑橢球單元

節(jié)點數(shù)N= 14 的光滑橢球單元如圖2 所示, 其中經(jīng)線和緯線構(gòu)成橢球表面的曲面坐標(biāo),節(jié)點13 和14 為極點,x3為旋轉(zhuǎn)對稱軸.參數(shù)空間如圖2(b)所示, 其中局部坐標(biāo)ξ1和ξ2分別對應(yīng)于經(jīng)線和緯線, 極點在參數(shù)平面中表達(dá)為雙實線.分別記N1,N2為經(jīng)線和緯線上節(jié)點的數(shù)目, 由于存在2 個極點, 則N1?2,N2分別是橢球單元上緯線和經(jīng)線的數(shù)目, 它們與節(jié)點總數(shù)N的關(guān)系為

圖2 光滑橢球單元及參數(shù)空間(N =14)Fig.2 Smooth ellipsoidal element with the parametric space (N =14)

用k1,k2對經(jīng)線和緯線上的節(jié)點進(jìn)行局部編號, 則單元整體編號m與局部編號的關(guān)系為

兩個極點的整體編號分別為m=N ?1 和m=N.在圖2(b)所示平面中, 正方形陰影部分即文獻(xiàn)[4]所構(gòu)建閉合單元的參數(shù)平面.注意到, 經(jīng)線都是開放的半圓或半橢圓, 兩個極點是所有經(jīng)線的共同端點, 當(dāng)ξ2=±1 所對應(yīng)的經(jīng)線就是閉合單元的端線.端點和端線的存在破壞了幾何量或場變量導(dǎo)數(shù)的連續(xù)性, 降低了插值精度.解決辦法是利用橢球的周期性, 在不增加單元實際節(jié)點數(shù)目的條件下, 沿緯線和經(jīng)線兩個方向上分別增加插值節(jié)點的數(shù)目, 構(gòu)建光滑橢球單元.由于緯線都是圓, 節(jié)點增加的方式以及插值多項式的形成與二維光滑橢圓單元完全相同, 緯線方向的形函數(shù)(ξ)仍可用式(7)表示(分別用k2和N2替換式中的k和N).以赤道為例(如圖3(a)), 整體編號為8, 5, 6 的節(jié)點是沿緯線方向重復(fù)使用兩次的節(jié)點, 括號中的數(shù)字為節(jié)點的局部編號, 其中重復(fù)使用的節(jié)點都有兩個局部編號.在圖2(b)所示的參數(shù)平面中, 沿緯線方向位置重合的節(jié)點用空心圓“?”表示.

圖3(b)為互為鏡像的兩條經(jīng)線及其節(jié)點的整體編號, 括號中的數(shù)字則為當(dāng)前經(jīng)線上節(jié)點的局部編號, 其中當(dāng)前經(jīng)線和鏡像經(jīng)線分別用實線和虛線表示, 點劃線為橢球的旋轉(zhuǎn)對稱軸.對于經(jīng)線方向, 在經(jīng)線兩端越過極點各增加s個插值節(jié)點(本工作取s=2), 使得兩個極點都位于單元的內(nèi)部, 不再是經(jīng)線的端點.沿經(jīng)線方向增加的節(jié)點用棱形符號“◇”表示(見圖3(b)),其插值多項式為

圖3 赤道與經(jīng)線上節(jié)點的整體與局部編號Fig.3 Global and local numberings for nodes along the equator and the meridians

當(dāng)前經(jīng)線所增加的插值節(jié)點的位置越過了極點, 從而位于虛線表示的鏡像經(jīng)線之上, 括號中出現(xiàn)的0 甚至負(fù)數(shù)是為了保持節(jié)點局部編號的連續(xù)性.由于利用了橢球的周期性, 光滑橢球單元上的一些節(jié)點被重復(fù)地使用, 重復(fù)使用節(jié)點的整體編號添加“′”以示區(qū)別(見圖2(b)).記經(jīng)線方向的形函數(shù)為

光滑橢球單元的形函數(shù)由緯線方向和經(jīng)線方向形函數(shù)的乘積來表達(dá), 即

式中:k2=1,2,··· ,N2;m由式(10)定義;M(k2) 表示鏡像函數(shù), 具有

由此可知, 光滑橢球單元的經(jīng)線數(shù)目或緯線的節(jié)點數(shù)目必須采用偶數(shù).對應(yīng)于兩個極點的形函數(shù)定義為

需要注意的是: ①單元極點處的形函數(shù)僅僅由經(jīng)線方向的形函數(shù)來定義, 所以位于兩個極點處的曲面的法線方向是不確定的; ②經(jīng)線方向增加的插值節(jié)點越過了極點位于其鏡像位置.圖2(b)所示參數(shù)平面中棱形符號“◇”所在區(qū)域上的法線方向與積分區(qū)間的法線方向是相反的, 但對單元積分并沒有任何影響, 因為光滑橢球單元沿緯線和經(jīng)線兩個方向上的積分區(qū)間均為[?1, +1], 即正方形陰影部分, 該積分區(qū)間與常規(guī)的閉合單元[4]完全一致.事實上, 由于本工作構(gòu)建的光滑橢球單元利用了橢球的周期性和對稱性, 在緯線和經(jīng)線兩個方向上增加了重復(fù)使用的插值節(jié)點, 提高了插值多項式的階數(shù), 提高了節(jié)點的利用率, 消除了端點效應(yīng), 從而能夠較大地提高插值的精度, 而單元的實際節(jié)點數(shù)目并沒有變化.

2 數(shù)值算例

2.1 光滑橢圓單元

橢圓的幾何參數(shù)如圖4 所示, 其中a和b分別為橢圓的長半軸和短半軸.分別采用二次單元、閉合單元、光滑單元進(jìn)行離散.用二次單元離散時, 節(jié)點總數(shù)與單元數(shù)之比為2∶1.采用高次單元進(jìn)行離散, 只需用一個單元.用3 種單元計算的圓面積(b/a= 1)相對誤差與節(jié)點總數(shù)的關(guān)系如圖5(a)所示.由圖5(a)可以看出, 兩種高次單元的計算精度都遠(yuǎn)高于二次單元, 而且隨著節(jié)點總數(shù)的增加, 高次單元相對誤差下降的速率也遠(yuǎn)高于二次單元.由于本工作構(gòu)建的光滑橢圓單元提高了插值多項式的階數(shù), 消除了端點效應(yīng), 光滑單元的計算精度比常規(guī)閉合單元大約提高了一個數(shù)量級.

圖4 橢圓的幾何參數(shù)Fig.4 Geometrical parameters of an ellipse

定義位置誤差為離散橢圓邊線與理想橢圓邊線的距離.采用兩種8 節(jié)點高次單元計算沿橢圓周邊的相對位置誤差, 結(jié)果如圖5(b)所示.由圖5(b)同樣可以看出, 光滑單元的計算精度比常規(guī)閉合單元大約提高了一個數(shù)量級.從整體的效果來看, 無論對于閉合單元還是光滑單元, 單元中部(ξ=0 附近)的精度優(yōu)于靠近兩端的精度.

圖5 面積與位置相對誤差的對比Fig.5 Comparisons of relative errors of areas and positions

采用8 節(jié)點光滑單元分別計算無限域中圓孔和橢圓孔在遠(yuǎn)場x2方向單位載荷(σ0=1)作用下的場變量, 即位移和應(yīng)力.圖6(a)為第一象限邊界上的位移與應(yīng)力計算值與理論解的比較, 圖6(b)為靠近邊界(沿圖4(a)第一象限中的虛線,d/a=0.01)處的位移與應(yīng)力計算值與理論值的比較.當(dāng)計算點靠近邊界時, 產(chǎn)生的近奇異積分采用距離變換的數(shù)值積分技術(shù)加以解決[6-7].圖6 表明本工作構(gòu)建的光滑橢圓單元的計算值與理論解吻合良好.

圖6 無限域中圓孔和橢圓孔的場變量Fig.6 Field variables on boundary of a circle and close to boundary of an ellipse

采用光滑橢圓單元在不同條件下對無限域中橢圓孔前方x2方向的應(yīng)力分量的分布進(jìn)行了計算, 結(jié)果如圖7 所示.由圖7 可知, 采用光滑橢圓單元進(jìn)行離散, 當(dāng)橢圓的形狀參數(shù)b/a降低時, 應(yīng)適當(dāng)增加保持精度所需要的節(jié)點數(shù)目.總體而言, 對于光滑單元, 采用數(shù)量不多的節(jié)點即可獲得滿意的精度, 例如圓孔(b/a=1)僅用了4 個節(jié)點, 這是其他類型的單元所做不到的.

圖7 橢圓孔前方的應(yīng)力分布Fig.7 Stress distributions in domain ahead of the elliptical holes

2.2 光滑橢球單元及應(yīng)用

首先比較傳統(tǒng)閉合單元與本工作構(gòu)建的光滑單元幾何參數(shù)的精度.圖8 為圓球表面積與體積的計算相對誤差隨著高次單元節(jié)點總數(shù)的變化情況, 其中采用96 個單元290 個節(jié)點的8 節(jié)點二次單元[8]進(jìn)行計算得到的表面積與體積各表達(dá)為一個點, 給出的水平虛線是為了便于誤差的比較.由圖8 可知, 兩種高次單元的計算精度都遠(yuǎn)遠(yuǎn)高于二次單元.由于提高了插值多項式的階數(shù)以及消除了端點效應(yīng), 光滑單元的計算精度比閉合單元提高了一到兩個數(shù)量級, 而且隨著節(jié)點總數(shù)的增加, 光滑橢球單元相對誤差下降的速率也高于傳統(tǒng)的閉合單元.

圖8 圓球表面積與體積的相對誤差對比Fig.8 Comparisons of relative errors of surface area and volume of ball

采用26 個節(jié)點對兩種高次單元計算的圓球半徑相對誤差與曲面法線誤差沿1/2 赤道(ξ1=1)和1/2 經(jīng)線(ξ2=1)的分布進(jìn)行計算分析, 結(jié)果如圖9 所示, 其中1/2 赤道和1/2 經(jīng)線分別對應(yīng)于緯線與經(jīng)線方向上半個單元的范圍(圖2(b)中的陰影區(qū)域).曲面法線誤差為

式中:n和n0分別為曲面單位法線向量的計算值與理論值.由圖9 可知, 半徑相對誤差與曲面法線誤差的極大值分別出現(xiàn)在節(jié)點之間區(qū)域與節(jié)點位置上, 插值多項式階數(shù)的提高和端點效應(yīng)的消除, 使得光滑單元的計算精度比閉合單元提高了一到兩個數(shù)量級.從整體的效果來看,無論對于緯線還是經(jīng)線方向, 單元中部(ξ1=0 與ξ2=0 附近)的精度優(yōu)于靠近兩端的精度, 與二維的情況類似.由于緯線都是封閉的圓, 而經(jīng)線都是開放的半圓或半橢圓, 沿緯線的擬合效果(見圖9(b))整體上要由優(yōu)于沿經(jīng)線的擬合效果(見圖9(a)).

圖9 圓球沿1/2 赤道和1/2 經(jīng)線的半徑相對誤差與曲面法線誤差Fig.9 Relative errors of radius and errors of surface normal along 1/2 equator and 1/2 maridian of ball

在應(yīng)用方面, 以解析解[9]為參照基準(zhǔn), 對旋轉(zhuǎn)橢球的Eshelby張量Sijkl進(jìn)行數(shù)值計算和比較.Eshelby張量的數(shù)值計算采用下式[8]進(jìn)行:

式中:Γ為橢球邊界;為面力基本解導(dǎo)出函數(shù);δij為kronecke 符號;xl為場點與源點距離在坐標(biāo)軸上的投影;μ,v分別為材料剪切模量與泊松比.令x3為橢球的旋轉(zhuǎn)對稱軸(見圖2(a)),c為對稱軸的半軸長,a為橢球x1和x2方向的半軸長.采用52 個節(jié)點的光滑橢球單元, 計算了旋轉(zhuǎn)橢球的7 個不同的Eshelby 張量隨c/a的變化規(guī)律, 計算結(jié)果與解析解的對比如圖10(a)所示.可以看出, 二者吻合良好.同時采用96 單元290 節(jié)點的二次單元對同一問題進(jìn)行計算, 比較二次單元與光滑單元計算結(jié)果相對于解析解的絕對誤差, 結(jié)果如圖10(b)所示.可以看出, 光滑單元的計算精度比二次單元提高了2～3 個數(shù)量級.

圖10 橢球Eshelby 張量和絕對誤差的比較Fig.10 Comparisons of Eshelby tensors and absolute errors of ellipsoids

在數(shù)值計算時, 每個二次單元采用8×8 點的高斯積分, 總的高斯積分點數(shù)為8×8×96=6 144.對于光滑單元, 當(dāng)c/a接近于1 時采用高斯積分點數(shù)為20×20 = 400; 當(dāng)c/a較小或較大(分別對應(yīng)于扁橢圓和長橢圓)時, 需要將積分域劃分為4 個極坐標(biāo)描述的三角形子域并采用近奇異積分技術(shù)[6-7]進(jìn)行數(shù)值積分, 總的高斯積分點數(shù)為20×8×4=640, 其中20 和8 分別為子域的徑向和角度方向的積分點數(shù).兩種單元的高斯積分點總數(shù)相差超過9 倍, 說明本工作構(gòu)建的光滑單元不僅有很高的計算精度, 也有很高的計算效率.可以預(yù)期, 將光滑單元用于粒子增強(qiáng)材料的數(shù)值模擬, 能夠提高該類材料性能模擬計算的精度和效率.

3 結(jié)束語

本工作通過將Lagrange 插值多項式轉(zhuǎn)化為嵌套積的形式, 實現(xiàn)了高次邊界單元形函數(shù)系數(shù)的計算機(jī)自動生成.在已有閉合單元的基礎(chǔ)上, 充分利用粒子幾何形狀的特性, 在不增加實際節(jié)點數(shù)目的條件下擴(kuò)大了參與插值的節(jié)點數(shù)目, 構(gòu)建了高次光滑邊界單元, 提高了插值多項式的階數(shù), 消除了端點效應(yīng).數(shù)值算例表明, 與傳統(tǒng)的二次單元和閉合單元相比, 高次光滑邊界單元能夠較大地提高橢圓和橢球粒子數(shù)值模擬的計算精度與效率, 適用于粒子增強(qiáng)材料的數(shù)值模擬.