求解固體力學問題的強-弱耦合形式單元微分法1)

胡 凱 高效偉 徐兵兵

(大連理工大學航空航天學院,工業裝備結構分析國家重點實驗室,遼寧大連 116024)

引言

在計算力學中,常用的數值算法基本上可以分為以下幾類:基于單元的網格類算法[1-2]、基于散點的無網格法[3-5]、基于邊界積分方程的邊界元法[6-7]、以及一些其他數值算法,例如比例邊界有限元法[8]、等幾何法[9]、有限塊法[10]等.這些算法又基本上可以分為兩類:一種是基于變分法或者加權余量法的弱形式算法,一種是基于配點法的強形式算法.基于加權余量法的弱形式算法通常有較高精度和更好的穩定性,如以上提到的廣泛使用的有限元法FEM.基于配點法的強形式算法具有更簡單的構造形式以及更快的計算效率,如無網格法中的有限點法[11]、徑向基點插值法[12-15]、自由單元法[16-17]等.基于分域算法的有限塊法以及等幾何配點法[18-20]也屬于強形式算法.但是配點法,特別是基于散點的無網格配點法的精度以及穩定性較差.

為了克服以上配點法的缺點,基于單元的強形式算法應運而生,如強形式有限元法[21-22]、節點梯度光滑有限元配點法[23]等.近年來,高效偉等[24-25]提出了一種求解二維及三維問題的新型強形式有限單元法——單元微分法EDM,并導出了二維及三維問題中形函數關于全局坐標的一階及二階偏導數的顯式表達式.目前,單元微分法已經用于求解靜力學問題、熱學問題[26-27]、動力學問題[28]、壓電問題[29]等多種單場及耦合場問題.與傳統配點法不同,單元微分法使用單元進行插值離散,因此具有更高的穩定性.與有限元法相比,由于不需要數值積分,因此具有更高的效率.當使用高階拉格朗日插值時,該算法具有較高的精度.但由于單元面中節點使用了應力平衡方程,對于一些模型中存在應力奇異的問題,如斷裂力學問題、V 型切口問題、間斷邊界問題等,單元微分法往往得到精度較差的結果.

為了解決此類問題,鄭永彤[30]等開發了弱形式的單元微分法,即對每個單元的內部點構造加權余量弱形式,并針對間斷邊界問題進行分析,得到了更精確的結果.但是由于每個單元都用到了數值積分,會帶來計算效率的降低.并且由于沒有改變界面點使用通量平衡這一方程,該算法在解決應力集中問題時仍存在一定問題.

為了解決以上問題,在本文中提出了將單元微分法與傳統有限元法相結合的耦合式算法.即對于模型的大部分,采用單元微分法,而對于模型中存在應力集中的部分以及采用強形式EDM 計算不準的區域,采用有限元法.通過這樣的處理,可以在保證計算效率的同時,提高算法的計算精度.本文將詳細描述單元微分法的基本原理,強-弱耦合算法的構造,以及該耦合算法在二維及三維力學問題中的應用,并對該耦合算法的精度,效率進行比較.

目前,求解固體力學問題中還未使用到強-弱耦合形式的單元微分法進行計算,因此對于該算法的研究非常有必要.

1 彈性力學基本方程

考慮一個處于平衡狀態的彈性體 Ω∈Rd,其計算域的邊界為 ?Ω=Γ=Γu∪Γt,其中 Γu為位移邊界,Γt為通量邊界,d為問題的維度.則彈性力學所滿足的控制微分方程如下:給定fi:Ω→Rd,gi:Γu→Rd,ti:Γt→Rd,尋找ui:→Rd,使得

其中,σij為柯西應力張量,fi為體積力,gi為給定位移,ti為模型邊界面力矢量,nj為 Ω 的外法線.在本文中,公式的重復指標代表求和.

應力張量和應變張量的關系可由以下廣義胡克定律給出

其中,εkl為應變張量,Di jkl為四階彈性本構張量,其表達式分別為

式中,uk為位移矢量,μ 為剪切模量,ν 為泊松比.將式(3)代入式(2),并利用本構張量的對稱性,將所得結果代入式(1)中,可得均質問題的彈性力學控制方程

及其面力邊界條件

2 單元微分法的基本原理

為了通過數值方法求解以式(5)及式(6)為控制方程的彈性力學問題,第一步需要對計算域進行分片離散.和傳統基于網格的算法類似,計算域 Ω 被離散為一系列互不重疊的四邊形或六面體有限單元

此外,未知變量f的p階導數可以表示為以下形式

其中,Nα稱之為單元的形函數,m為單元 Ωe內的節點數量.在本文中,為了減小矩陣帶寬及計算量,式(8)～式(10)中變量N通常取3.

接下來的問題是要處理不規則計算域的問題.對于幾乎所有的物理問題,其計算域往往并不是規則的.因此,需要用到坐標轉換技術,即建立函數對整體坐標的導數與函數對參數坐標的導數之間的關系.考慮計算域 Ωe內的某連續函數f,其一階導數可通過以下方式獲得

對于強形式算法,需要用到函數對全局坐標的二階導數.為了獲得單元的二階導數,通過對式(12)進行再次微分,可得函數對全局坐標的二階導數的顯式表達式為

其中Jik為坐標變換的雅克比矩陣,其逆矩陣可通過下式計算

此外,式(13)中雅克比的逆矩陣關于局部坐標的導數為

將式(14)及式(15)代入式(12)及式(13)中,并考慮式(11),可得變量關于全局坐標的一階導數及二階導數為

從式(16)及式(17)中可以看出,通過計算形函數關于全局坐標的一階及二階導數,可得函數f關于全局坐標的導數.以上變量關于全局坐標的一階導數及二階導數在單元微分法及有限元法均有相關應用.針對二階導數計算式(17),其相關具體每一項的展開式可見EDM 的文獻[24-25].

3 強-弱耦合形式算法

現有的數值算法可以分為兩大類:基于強形式的配點類算法和基于弱形式的伽遼金類算法.單元微分法是一種強形式算法.和配點類無網格法不同的是,單元微分法基于類似有限元法中的單元,而非散點,因此算法的穩定性要高于傳統的無網格法.由于不需要進行積分,其易用性、效率等方面都優于傳統有限元法.但是作為一種強形式算法,針對具體問題要達到較高的計算精度,往往需要更多的網格,尤其是針對具有應力集中的問題.為了解決以上問題,可以嘗試將單元微分法和有限元法進行耦合計算,即在利用單元微分法優點的前提下,提高其在求解具體問題時的能力.如圖1 所示為一個具有切口的二維模型的網格圖,其在切口附近的單元使用有限元法,其他部分的單元采用單元微分法.

圖1 單元微分法和有限元法耦合形式的模型離散方案Fig.1 Model discretization scheme in coupling form of element differential method and finite element method

3.1 單元微分法

作為一種強形式的有限元法,單元微分法和其他強形式的無網格法的基本思想類似,通過將形函數關于全局坐標的一階導數及二階導數直接代入到控制方程及邊界條件中構造系統方程.和強形式無網格法不同的是,該算法基于網格而非散點,因此穩定性更好.

如圖2 中所示,在單元微分法中,將模型中的點分為兩類:處于單元內部的點,稱之為內部點;處于單元邊界的點,稱之為邊界點.對于內部點來說,其滿足具有二階導數形式的控制方程(式(5)),對于單元邊界點,其滿足具有一階導數形式的面力平衡方程(式(6)).

圖2 單元內部點及單元邊界點Fig.2 Element internal nodes and element interface nodes

對于內部點,將函數關于全局坐標二階導數的微分關系式(17)代入式(5)中,可得以下方程

在上式中,出現了形函數對全局坐標的二階導數.這部分公式的推導見第二章.

對于單元邊界點I,其通常與NI個單元相連.對于任意一個相連的單元e,單元邊界點I與單元e的SeI個單元面相關.根據所處位置不同,單元邊界點可分為兩類:模型內部點和模型邊界點.對于模型內部的單元界面點,其滿足通量平衡方程,即

對于模型邊界節點,其滿足的方程和模型內部界面點方程(式(19))類似,不同的是方程右端項為指定面力,即對于模型邊界點,有

對于一個EDM 單元,可依次根據式(18)～式(20)對單元的節點進行方程組集,形成單元系數矩陣,并可組裝成整體剛度矩陣.

3.2 有限單元法

有限元法是一種弱形式算法.在本文提到的強弱耦合形式的算法中,將采用伽遼金有限元法與單元微分法相結合的形式求解問題.在本文中,將使用加權余量法簡要推導相關方程.

對于單元中的每個配置點,權函數w取為計算點所在單元 Ωe的插值形函數N*.考慮彈性力學控制方程(5),可得到以下積分形式

對于上式,考慮分部積分,采用高斯散度定理,并考慮面力邊界條件(6)和形函數離散式(11),可得

其中 Γe為單元 Ωe的邊界.從上式可以看出,伽遼金有限元法中只出現了形函數對全局坐標的一階導數,可采用第2 節中推導的公式直接計算.對單元中每個點進行計算,可得有限元的單元剛度陣.

3.3 耦合形式算法

從上兩節中可以看出,單元微分法在計算形成系統方程組的過程中,采用與無網格法中的配點形式相似的方法,不需要對控制方程進行積分變換,但需要計算單元中函數對整體坐標的二階偏導數.而對于有限元法,通過分部積分,只需要計算單元中函數對整體坐標的一階偏導數.并且由于使用積分,計算精度較單元微分法高,但是FEM 計算需要用到高斯數值積分,計算量較大.通過該文中提到的耦合形式算法,可以較好地解決以上問題.

對于如圖1 所示的模型,其分為兩種域:單元微分域和有限元域,在不同域內使用不同的數值算法.但是需要注意的是,對于兩個不同域之間的界面節點(圖1 中紅色點),為了更方便進行計算,強制其滿足式(22)的伽遼金弱形式,即對于界面的單元微分單元是一種雜交元.這種耦合形式是簡單明了的,不需要對強-弱耦合界面做特殊操作,也不需要對總體系數矩陣做特殊變換.

下面給出具體的剛度矩陣K以及載荷向量b的表達式,如表1 中所示.

表1 兩種方法的具體剛度矩陣K 及載荷向量b 表達式Table 1 The specific stiffness matrix K and load vector b expressions of the two methods

通過以上形式,可得最終的系統方程組為

通過引入第一類邊界條件,求解以上方程組,可以求得節點的位移值,并可進一步計算應力.

4 數值算例

通過本文上述所推導出的公式,將其編寫成通用的Fortran 程序,并借助下面一些彈性力學的分析計算來驗證單元微分法與有限單元法耦合算法的正確性.

4.1 頂部受拉的V 型缺口平板

含有V 型缺口的平板問題是一個典型的具有應力奇異點的問題.對于強形式算法,求解此類問題所得結果往往精度較差.在使用本文所給出的耦合形式算法中,在切口附近區域采用有限元算法,其他部位采用單元微分法,可以求得滿意的結果.

如圖3 所示的模型為一矩形平板,長度L=100 mm,寬度D=50 mm,右側有一V 型缺口,缺口角度為30°,寬5 mm.模型下端固定,上端受y軸正向均布載荷.

圖3 頂部受拉的V 型缺口平板模型Fig.3 V-notch plate model with tension at the top

在本問題中,結構的彈性模量為200 GPa,泊松比為0.3,模型上端的載荷參數為P=100 N.在該算例中,分別利用EDM,ANSYS,EDM_FEM 耦合的方法進行計算.在使用EDM 計算時,采用二階拉格朗日單元,單元總數為880 個,總節點數為3653 個節點,如圖3 所示.同時,有限元采用商業軟件ANSYS進行計算,所用單元類型為8 節點單元,網格密度和EDM 一樣,單元總數880 個,節點總數為2773 個.

首先比較位移計算結果.經過計算,三種不同算法在相同網格密度下所得的y方向最大位移如表2所示.

從表2 中可以看出,傳統的EDM 在計算該類問題時,所得位移結果誤差較大,相較于FEM 密網格所得參考值的相對誤差為3.37%.而對于文本所提出的耦合算法EDM_FEM,計算所得最大位移相對FEM 密網格所得參考值的誤差僅為0.045%,可以證明該耦合算法在求解這類問題時的精度較為可靠.

表2 三種方法所得y 方向最大位移(mm)結果對比Table 2 Comparison of maximum displacement (mm)in y direction obtained by three methods

需要注意的是,針對該類問題,即使是有限元法,在不同網格密度下也會得到相差較大的結果.因此,為了進行結果的比較,采用足夠密的有限元網格計算作為參考.在本文中采用ANSYS 進行計算,所用有限元網格單元數量為81 600 個.

為了比較位移計算結果,選取模型頂部y=L和中部y=L/2 處節點的位移進行比較,三種方法所得結果繪制于圖4 和圖5.從圖中看出,FEM 和EFM_FEM 所得結果十分吻合,并和參考值吻合較好.可見模型頂部的y方向位移能夠得到比較接近實際的結果,耦合的求解方法大大提升了計算的準確性.

圖4 頂部節點y 方向的位移Fig.4 Displacement of top node in y direction

此外,從圖5 可以看出,相比于傳統的EDM 形式,在應力集中區域,耦合的求解方法可以得到更加精確的結果.同時,繪制了y=L/2 線上的馮·米塞斯應力云圖,如圖6 所示.

圖5 y=L/2 線上y 方向的位移Fig.5 Displacement in y direction on y=L/2 line

圖6 y=L/2 線上馮·米塞斯應力Fig.6 von-Mises stress on y=L/2 line

從圖6 中能夠非常容易地看到,在V 形切口附近,模型的應力變化十分劇烈.實際上在切口尖端所得應力值應為無窮大,即應力奇異.但是在進行數值計算時,由于數值離散誤差,該點附近應力變化會被磨平.在計算該類問題時,強形式算法雖然能得到較大的應力結果,但是該值的可信度不高,如應力強度因子計算值并不精確.

下面接著給出三種方法所計算的全域應力云圖,如圖7 所示.從圖中可以看出,相較于傳統的EDM,本文所提耦合形式算法計算所得尖端應力場更光滑.

圖7 兩種方法的馮·米塞斯應力云圖Fig.7 von-Mises stress cloud maps of the two methods

就該問題而言,對于應力集中點附近的計算結果,有限元法(F E M)和本文提出的耦合算法(EDM_FEM)所得馮·米塞斯應力結果更接近于參考值.相較于傳統的EDM 方法,本文提出的耦合算法(EDM_FEM)能夠較大地提高算法的計算精度,得到理想的計算結果.

4.2 超燃沖壓發動機燃燒室模型

作為一種強形式算法,單元微分法在形成系數矩陣時不需要進行積分計算,因此在對大規模問題進行計算時,可以節省大量的計算時間.但在用EDM 求解大規模問題時,在模型較關鍵部位需要用到更密的網格,才能得到較精確結果.在利用本文提出的耦合算法求解大規模問題時,可在模型的大部分采用EDM 進行計算,在關鍵部件采用FEM 進行計算.通過該方法,可以在得到較精確的結果的同時,提高問題的計算效率.

本算例分析一個較為復雜的三維問題,模型取自超燃沖壓發動機中的燃燒室.由于模型對稱,僅采用一半的模型進行分析,其結構如圖8 所示,并在圖中給出了模型的關鍵尺寸參數.

圖8 燃燒室模型及其重要尺寸(單位:mm)Fig.8 Combustion chamber model and its important parameters(unit:mm)

在本算例中,如圖9(a)中標號①所示為燃燒室外部框架完全固定,燃燒室內壁面受0.3 MPa 的壓力,如圖9(b)中標號②所示.結構的對稱面(圖9 中標號③所示)只固定法向的位移,兩個切向的位移不受約束.

圖9 燃燒室模型的邊界條件Fig.9 Boundary conditions of combustion chamber model

該結構所用材料的彈性模量為200 GPa,泊松比為0.3,網格劃分情況如圖10 所示.可以發現模型的尾部壁面較薄,在使用純EDM 求解時需用到大量網格,因此需在求解時將該區域利用有限元求解作為補充.由于采用兩個數值方法的耦合求解,因此在圖形中標注出不同算法的求解區域,綠色區域為FEM求解區域,灰色區域為EDM 求解區域.為便于比較,三種方法均采用如圖10 所示的同一套網格.

圖10 燃燒室模型的網格Fig.10 Grid of combustion chamber model

網格劃分為70 913 個三維27 節點單元,總節點數為641 577 個,對于此模型分別采用EDM,FEM以及EDM_FEM 的方法求解,并對如圖9(b)中AB線上的節點的總位移和馮·米塞斯應力進行比較,以驗證該數值求解方法的準確性.所得比較結果如圖11 和圖12 中所示.

圖11 三種方法在AB 線上的總位移Fig.11 Total displacement of three methods on AB line

圖12 三種方法在AB 線上的馮·米塞斯應力Fig.12 von-Mises stress of three methods on AB line

從圖11 和圖12 中可以看出在模型尾部,傳統的EDM 方法計算精度較差,而強-弱耦合的求解方法卻可以得到較好的計算精度,所得指定路徑上的位移和應力結果和傳統有限元法十分吻合.

不僅如此,在保證計算精度的情況下,大大提高了求解的速度.為了更清晰地比較整體分析的速度,在程序中記錄了組集方程及求解所用CPU 時間,并列于表3.其中CPU 采用i9-9900 k 3.60 GHz,RAM 為128 GB.

表3 三種方法組集方程及求解所用時間對比Table 3 Time comparison of three methods

由于三種方法求解方程的時間相近,不進行比較.可見強形式的EDM 方法由于不需要對控制方程進行積分,極大地降低了形成單元剛度矩陣并組裝整體剛度矩陣所用的時間.在該方面,EDM 組集方程所用的時間不到有限元法的1/13.

5 結論

單元微分法是一種強形式有限元法,其在計算時不需要用到積分操作,具有簡單高效等優點.但是其在處理存在應力奇異問題時表現不佳,因此,本文提出了將單元微分法與伽遼金有限元法相結合的強-弱耦合算法.針對傳統單元微分法本文介紹了單元微分法和傳統有限單元法相結合的耦合式算法,并給出算例驗證其精度和計算效率,可得出下列結論.

(1)本文采用EDM,FEM 及耦合的求解方法進行比較分析,通過比較模型特定部位的應力和位移數據,可較為清晰地發現耦合求解方法的適用性更強.

(2)由于不需要數值積分,強形式的單元微分法在組集系統方程中所消耗的時間大為減少,有利于降低總的求解時間.

(3)由于采用單元微分法與傳統有限元法耦合的方法,使能夠在滿足計算精度的同時顯著減少計算所消耗的時間.

需要說明的是,配點法在處理力邊界時比較麻煩,而且精度不高.因此可以進一步考慮在模型邊界部分采用有限元離散,提高計算精度.