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

胡 凱 高效偉 徐兵兵

(大連理工大學(xué)航空航天學(xué)院,工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室,遼寧大連 116024)

引言

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

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

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

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

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

1 彈性力學(xué)基本方程

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

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

應(yīng)力張量和應(yīng)變張量的關(guān)系可由以下廣義胡克定律給出

其中,εkl為應(yīng)變張量,Di jkl為四階彈性本構(gòu)張量,其表達式分別為

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

及其面力邊界條件

2 單元微分法的基本原理

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

此外,未知變量f的p階導(dǎo)數(shù)可以表示為以下形式

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

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

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

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

此外,式(13)中雅克比的逆矩陣關(guān)于局部坐標的導(dǎo)數(shù)為

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

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

3 強-弱耦合形式算法

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

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

3.1 單元微分法

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

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

圖2 單元內(nèi)部點及單元邊界點Fig.2 Element internal nodes and element interface nodes

對于內(nèi)部點,將函數(shù)關(guān)于全局坐標二階導(dǎo)數(shù)的微分關(guān)系式(17)代入式(5)中,可得以下方程

在上式中,出現(xiàn)了形函數(shù)對全局坐標的二階導(dǎo)數(shù).這部分公式的推導(dǎo)見第二章.

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

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

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

3.2 有限單元法

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

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

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

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

3.3 耦合形式算法

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

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

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

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

通過以上形式,可得最終的系統(tǒng)方程組為

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

4 數(shù)值算例

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4.2 超燃沖壓發(fā)動機燃燒室模型

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

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

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

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

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

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

圖10 燃燒室模型的網(wǎng)格Fig.10 Grid of combustion chamber model

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

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

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

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

不僅如此,在保證計算精度的情況下,大大提高了求解的速度.為了更清晰地比較整體分析的速度,在程序中記錄了組集方程及求解所用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 結(jié)論

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

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

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

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

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