田 鈺, 和東宏, 馬 杭

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

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

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

1 高次光滑邊界單元

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

邊界單元形函數(shù)的形成通?；贚agrange 插值多項(xiàng)式.邊界上的幾何量或場(chǎng)變量可用節(jié)點(diǎn)上的數(shù)值近似表達(dá)為

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

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

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

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

因?yàn)榍短壮朔e的形式具有最高的計(jì)算效率[5], 容易寫出多項(xiàng)式的導(dǎo)函數(shù), 同樣可通過(guò)計(jì)算機(jī)程序自動(dòng)加以計(jì)算:

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

1.2 光滑橢圓單元

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

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

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

1.3 光滑橢球單元

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

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

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

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

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

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

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

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

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

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

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

2 數(shù)值算例

2.1 光滑橢圓單元

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

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

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

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

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

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

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

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

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

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

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

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

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

圖9 圓球沿1/2 赤道和1/2 經(jīng)線的半徑相對(duì)誤差與曲面法線誤差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), 對(duì)旋轉(zhuǎn)橢球的Eshelby張量Sijkl進(jìn)行數(shù)值計(jì)算和比較.Eshelby張量的數(shù)值計(jì)算采用下式[8]進(jìn)行:

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

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

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

3 結(jié)束語(yǔ)

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