田 鈺, 和東宏, 馬 杭

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

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

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

1 高次光滑邊界單元

1.1 形函數系數

邊界單元形函數的形成通?；贚agrange 插值多項式.邊界上的幾何量或場變量可用節(jié)點上的數值近似表達為

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

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

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

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

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

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

1.2 光滑橢圓單元

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

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

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

1.3 光滑橢球單元

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

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

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

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

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

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

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

光滑橢球單元的形函數由緯線方向和經線方向形函數的乘積來表達, 即

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

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

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

2 數值算例

2.1 光滑橢圓單元

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

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

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

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

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

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

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

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

2.2 光滑橢球單元及應用

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

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

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

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

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

在應用方面, 以解析解[9]為參照基準, 對旋轉橢球的Eshelby張量Sijkl進行數值計算和比較.Eshelby張量的數值計算采用下式[8]進行:

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

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

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

3 結束語

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