熱力耦合問題數(shù)學(xué)均勻化方法的物理意義

朱曉鵬， 黃俊， 陳磊,*， 邢譽峰

(1. 安徽華電工程咨詢設(shè)計有限公司， 合肥 230022； 2. 北京航空航天大學(xué)航空科學(xué)與工程學(xué)院， 北京 100083；3. 北京航空航天大學(xué) 合肥創(chuàng)新研究院， 合肥 230012)

復(fù)合材料具有比強度高、比剛度大等優(yōu)點,廣泛應(yīng)用于航天、航空工業(yè)領(lǐng)域。眾所周知，對于很多復(fù)合材料的宏觀解，如低階頻率和模態(tài)，可以使用等應(yīng)變模型或等應(yīng)力模型[1]及其他均勻化方法[2]求解，但相對于宏觀應(yīng)力分析，細觀結(jié)構(gòu)分析要復(fù)雜很多。為了在計算精度和效率之間達到平衡，各種多尺度方法相繼被提出，如數(shù)學(xué)均勻化方法(MHM)[3-4]、廣義有限單元方法(GFEM)[5-6]、多尺度有限單元方法(MsFEM)[7-8]、異類多尺度方法(HMM)[9-10]及多尺度特征單元方法(MEM)[11-12]。其中，MHM是最具代表性的多尺度方法之一，因其具有嚴格的數(shù)學(xué)邏輯，且可以包含復(fù)合材料所有的細觀結(jié)構(gòu)信息,廣泛應(yīng)用于分析熱力耦合作用下的周期復(fù)合材料結(jié)構(gòu)響應(yīng)問題[13-18]。

在應(yīng)用MHM處理熱力耦合作用下的復(fù)合材料結(jié)構(gòu)響應(yīng)問題時，攝動階次的選擇是不可回避的問題，然而多數(shù)文獻只考慮了一階攝動項，這對于細觀結(jié)構(gòu)物理和力學(xué)信息的捕捉往往是不夠的[19]，在很多實際工程問題中[20-25]，二階攝動項對于計算精度的影響不可忽略，這就需要尋找高階MHM求解復(fù)合材料細觀結(jié)構(gòu)的物理和力學(xué)屬性。Han等[26-27]基于二階均勻化方法處理功能梯度材料的靜態(tài)熱力耦合問題，準確估計了有效彈性模量和應(yīng)力應(yīng)變場，并基于周期均勻化方法推導(dǎo)了用于求解復(fù)合材料彈性問題的二階攝動項格式。Guan[20]和Cui[22]等給出了一種二階多尺度方法用于預(yù)測周期復(fù)合材料的物理和力學(xué)屬性，并應(yīng)用于實際工程問題中。Allaire和Habibi[23]使用雙尺度二階漸進展開方法研究了具有輻射邊界條件下的熱傳導(dǎo)模型，并證明了均勻化過程中的雙尺度收斂。通過二階攝動項可以更加準確地捕捉到材料內(nèi)部物理和力學(xué)行為的細觀信息，這一結(jié)論得到越來越多學(xué)者的認可，然而關(guān)于二階攝動項的必要性問題基本上都是通過具體物理和力學(xué)問題的計算精度進行論證和說明，其結(jié)論往往不具有普遍性。

MHM本身是具有嚴密邏輯的數(shù)學(xué)方法，在處理力學(xué)領(lǐng)域具體工程問題時，如何從力學(xué)角度解釋MHM是應(yīng)用數(shù)學(xué)學(xué)者經(jīng)常關(guān)心的問題，一旦能夠研究清楚各階攝動項的物理意義，無疑將有助于從力學(xué)角度確定合適的展開項數(shù)或階次，避免一味追求高階次或盲目舍棄高于某階的項。Xing和Chen[28]在分析比較各種多尺度方法基礎(chǔ)之上，針對線彈性平面問題揭示了MHM前三階攝動項的物理意義，為該方法在數(shù)學(xué)和力學(xué)之間建立了聯(lián)系，并通過物理意義的分析得到了二階攝動項與一階攝動項同為攝動主項的結(jié)論。除此之外，鮮有其他關(guān)于MHM物理意義研究的報道。相對于線彈性問題MHM物理意義的研究，熱力耦合問題MHM各階攝動項中同時包含彈性影響函數(shù)和熱影響函數(shù)，以致出現(xiàn)更多、更復(fù)雜的虛擬載荷形式，其物理解釋也變得復(fù)雜，然而確定熱力耦合問題MHM的物理意義將有助于更好地理解和應(yīng)用該方法處理具體的力學(xué)和物理問題，并且從物理機理的角度研究合適的展開階次所得到的結(jié)論相對于通過具體力學(xué)和物理問題的研究所得到的結(jié)論更具有普遍性。

在此背景下，本文首先通過構(gòu)造MHM各階攝動項的全解耦格式，推導(dǎo)熱力耦合問題MHM高階攝動項的數(shù)學(xué)表達式及矩陣列式；其次將彈性影響函數(shù)和熱影響函數(shù)分別比擬為彈性虛擬位移和熱虛擬位移，利用熱力耦合問題MHM矩陣列式實現(xiàn)方法、量綱分析方法和虛擬載荷自平衡性質(zhì)，分別給出熱特征位移和彈性特征位移的空間矢量圖，再借助構(gòu)造有限元位移場和溫度場中廣義坐標的概念，給出了各階彈性虛擬位移、熱虛擬位移和攝動項的物理意義，并通過物理意義的分析得到了二階攝動項的必要性依據(jù)；最后通過數(shù)值算例和最小勢能泛函計算結(jié)果驗證了熱力耦合問題MHM有限元矩陣列式推導(dǎo)和物理意義分析的正確性。

1 熱力耦合問題MHM的數(shù)學(xué)表達式和有限元矩陣列式

基于微觀結(jié)構(gòu)的周期性和單胞域的一致性假設(shè)，均勻化理論可以解耦成為微觀結(jié)構(gòu)異質(zhì)邊界值的單胞問題和宏觀問題，對于結(jié)構(gòu)線性小變形問題，熱力耦合的應(yīng)力-應(yīng)變關(guān)系為

(1)

總應(yīng)變-位移和熱應(yīng)變分別表示為

(2)

(3)

式中：uε為真實位移；amn為熱膨脹系數(shù)；T為結(jié)構(gòu)溫差。

將式(2)和式(3)代入式(1)，并考慮彈性常數(shù)張量的對稱性可以得到

(4)

二維周期復(fù)合材料控制方程為一致橢圓偏微分方程，其具有Dirichlet邊界條件[29]：

(5)

(6)

表1 攝動位移及影響函數(shù)控制方程Table 1 Perturbation displacements and governing equation of influence functions

3) 三階及以上階次影響函數(shù)控制方程具有相同的遞推形式。

熱力耦合問題MHM的應(yīng)用分3步實施：①求解單胞問題，得到各階熱影響函數(shù)、彈性影響函數(shù)及均勻化材料參數(shù)EH和aH；②求解各階宏觀位移場導(dǎo)數(shù)和溫度場導(dǎo)數(shù)；③求解MHM真實位移和應(yīng)力。

表1中各階彈性影響函數(shù)和熱影響函數(shù)控制方程的矩陣列式可以歸納為

(7)

(8)

(9)

(10)

(11)

式中：

(12)

其中：s=3,4,5…；e為單胞內(nèi)的子單元數(shù)量；Ye表示單胞內(nèi)子單元的區(qū)域；式(8)和式(9)中Eε表示周期結(jié)構(gòu)內(nèi)不同材料的彈性常數(shù)矩陣，式(10)和式(11)中的彈性常數(shù)矩陣在文獻[28]附錄A中給出；B為應(yīng)變矩陣；χ為彈性影響函數(shù)矩陣；ψ為熱影響函數(shù)矩陣；N為形函數(shù)矩陣；B、χ、ψ、N的維度取決于使用有限元(FEM)計算時所選擇的單元階次。

以線性矩陣單元為例，N為2×8矩陣，形式如下：

(13)

式中：Ni(i=1,2,3,4)為單胞內(nèi)子單元4個節(jié)點的形函數(shù)。此時B為3×8矩陣，與標準有限元方法形式相同；一階影響函數(shù)χ1和ψ1分別為8×3矩陣和8×1矩陣，二階影響函數(shù)χ2和ψ2分別為8×6矩陣和8×2矩陣，三階影響函數(shù)χ3和ψ3分別為8×12矩陣和8×4矩陣，更高階次影響函數(shù)矩陣的維度以此類推。

(14)

結(jié)構(gòu)由于熱膨脹作用只產(chǎn)生線應(yīng)變，剪應(yīng)變?yōu)?，所以這種由于熱變形產(chǎn)生的應(yīng)變可以看作結(jié)構(gòu)的初始應(yīng)變ε0。對于二維問題，ε0的表達式為

(15)

所以均勻化宏觀問題可以表達為

(16)

式中：f為外載荷。

在實際工程應(yīng)用中，一般將一個單胞或多個單胞作為一個宏觀單元處理，提高計算效率。

在得到宏觀位移及其各階導(dǎo)數(shù)、各階溫差場導(dǎo)數(shù)及攝動位移后，結(jié)構(gòu)真實細觀位移可以使用如下漸進展開表達式得到：

(17)

式中：?u0/?x、?2u0/?x2及?3u0/?x3分別與χ1、χ2和χ3的上標kl、klp和klpq的組合相對應(yīng)求得；?T/?x和?2T/?x2分別與ψ2、ψ3的上標klp和klpq的組合相對應(yīng)求得。

從以上求解公式可以看出，熱力耦合問題MHM解耦形式都可以使用有限元方法得到。由式(17)可以看出，MHM的顯著優(yōu)點為細觀尺度的解可以用宏觀尺度來表示，這意味著在得到結(jié)構(gòu)單胞的影響函數(shù)后，其計算精度取決于宏觀位移、宏觀位移場導(dǎo)數(shù)和宏觀溫度場導(dǎo)數(shù)的求解精度。

2 熱力耦合問題MHM的物理解釋

2.1 二維周期復(fù)合材料單胞結(jié)構(gòu)

對于單胞模型內(nèi)部一個線性子單元，用于求解一階彈性影響函數(shù)χ1的彈性自平衡虛擬載荷形式為

(18)

式中：ν為泊松比；l1和l2分別為矩形子單元沿著坐標軸x1和x2方向的長度。顯然，如果l1=l2或c=1，即正方形子單元，F(xiàn)1的第1列(與kl=11相對應(yīng))中每個節(jié)點沿著坐標軸x1和x2方向的比值為±1/ν，第2列(與kl=22相對應(yīng))中每個節(jié)點沿著坐標x1和x2方向的比值為±ν，第3列(與kl=12相對應(yīng))中每個節(jié)點沿著坐標x1和x2方向的比值為±1。

同理，對于單胞模型內(nèi)部一個線性子單元，用于求解一階彈性影響函數(shù)ψ1的熱自平衡虛擬載荷形式為

(1+ν)aε-(1+ν)caε(1+ν)aε

(1+ν)caε-(1+ν)aε(1+ν)caε]T

(19)

圖1為包含3塊夾雜的二維單胞細分網(wǎng)格模型。為了更加清晰地給出虛擬載荷矩陣中每一列矢量圖，在本文中選取夾雜和基體的彈性模量分別為50 GPa和60 GPa，熱膨脹系數(shù)分別為1.5×10-6/℃和1×10-6/℃，泊松比均為0.2。前三階虛擬載荷和影響函數(shù)的物理意義討論如下：

① 只要夾雜的材料性質(zhì)相同，在每一塊夾雜上相同位置節(jié)點上的矢量都是相同的，這個結(jié)論普遍成立，與單胞包含夾雜的數(shù)量無關(guān)。

③ψ1和χ1反映了單胞內(nèi)材料的不連續(xù)性，因為它們是由分布在夾雜和基體交界處節(jié)點的分

圖1 包含3塊夾雜的二維周期復(fù)合材料單胞結(jié)構(gòu)Fig.1 2D periodical composites unit cell structures with 3 inclusions

段線性載荷計算得到的。另一方面也反映了夾雜和基體之間的相互作用，是單胞內(nèi)材料不均勻性的最基本關(guān)系，因為在基體和夾雜的內(nèi)部節(jié)點是沒有載荷作用的。

④u1為一階攝動主項，本質(zhì)上是χ1矩陣的3列和ψ1的1列虛擬位移的線性組合，所以一般來說，如果基體和夾雜的彈性模量、熱膨脹系數(shù)相差不大時，使用MHM處理周期結(jié)構(gòu)復(fù)合材料熱力耦合問題攝動到一階的精度就足夠了。

圖2 一階虛擬載荷矢量圖Fig.2 First-order virtual load vector

圖3 二階虛擬載荷矢量圖Fig.3 Second-order virtual load vector

圖4 三階虛擬載荷矢量圖Fig.4 Third-order virtual load vector

綜合前三階虛擬載荷矢量示意圖可以得到以下結(jié)論：

1) 一階彈性影響函數(shù)χ1和一階熱影響函數(shù)ψ1反映的是單胞內(nèi)基體和夾雜之間的相互作用的結(jié)果，一階攝動項u1是χ1矩陣中3列向量和ψ1的線性組合，反映單胞細觀信息的位移場；包含一階攝動項的MHM計算精度對于基體和夾雜之間彈性模量和熱膨脹系數(shù)相差不大的情況下是足夠精確的。

3) 一階影響函數(shù)χ1和ψ1反映了基體內(nèi)部、夾雜內(nèi)部及基體和夾雜相互作用的結(jié)果，二階攝動項u2是χ2矩陣中6列向量線性組合與ψ2矩陣中2列向量線性組合的和，反映單胞內(nèi)材料更加準確的細觀信息，所以對于一般周期復(fù)合材料而言，包含二階攝動項的MHM的計算精度是足夠的，但對于單胞內(nèi)的基體和夾雜材料參數(shù)相差特別大的周期復(fù)合材料，極限情況下，如具有孔洞的周期結(jié)構(gòu)，將孔洞視為夾雜，則基體和夾雜的材料彈性模量比例無限大，或某些顆粒增強結(jié)構(gòu)，基體和夾雜的彈性模量比值無限小，此時包含二階攝動項的計算精度是否足夠需要專門研究。

2.2 周期復(fù)合材料桿單胞結(jié)構(gòu)

文獻[28]中給出了線性單元桿單胞彈性虛擬載荷和彈性特征位移分析結(jié)果，本文使用二次單元分析桿單胞的彈性影響函數(shù)和熱影響位移，以及與之相對應(yīng)的彈性虛擬載荷和熱虛擬載荷。

考慮包含一塊夾雜的桿單胞，如圖5所示，對于這樣簡單的問題可以推導(dǎo)出其解析解形式。使用3個二次桿單元計算該單胞模型，該情況下，m=k=l=p=q=1。E1和E2分別表示基體和夾雜的彈性模量，a1和a2分別為基體和夾雜的熱膨脹系數(shù)，L為單元長度，A為橫截面積。

圖5 包含一塊夾雜的周期復(fù)合材料桿單胞結(jié)構(gòu)Fig.5 Periodical composites rod unit cell structure with one inclusion

針對圖5中的一個子單元，各階影響函數(shù)控制方程如下：

1) 一階熱影響函數(shù)控制方程

(20)

2) 一階彈性影響移函數(shù)控制方程

(21)

(22)

(23)

分別求解方程可得到單胞結(jié)構(gòu)彈性影響函數(shù)和熱影響函數(shù)為

(24)

(25)

值得注意的是，當(dāng)相鄰2個單元的彈性模量、熱膨脹系數(shù)均相同時，即均質(zhì)材料，式(24)右端項為0，同時對于非均質(zhì)材料，若耦合參數(shù)相同(E2a2=E1a1)時，方程右端項也為0；對比式(25)，彈性自平衡虛擬載荷只有在不同彈性模量交界處才非0，即非均質(zhì)材料，而熱自平衡虛擬載荷只有在不相同耦合模量交界處才非0，基體和夾雜內(nèi)部沒有虛擬載荷，所以一階虛擬載荷反映的是基體和夾雜之間的相互作用，是反映單胞內(nèi)部結(jié)構(gòu)最基本的信息載體。

式(12)、式(24)、式(25)聯(lián)立，可得得到均勻化彈性模量和均勻化熱膨脹系數(shù)分別為

(26)

二階熱虛擬載荷和彈性虛擬載荷分別為

(27)

二階彈性影響函數(shù)和熱影響函數(shù)分別為

(28)

三階熱虛擬載荷和彈性虛擬載荷分別為

(29)

三階熱影響函數(shù)和彈性影響函數(shù)分別為

(30)

從二階虛擬載荷和三階虛擬載荷形式可以看出，當(dāng)彈性模量相同時，彈性虛擬載荷為0，當(dāng)耦合參數(shù)相同時，熱虛擬載荷為0，即在非均質(zhì)材料交界處，彈性虛擬載荷一定非0，而熱虛擬載荷可能為0，所以彈性問題只考慮不同材料之間彈性參數(shù)的異同，而熱問題需要同時考慮彈性參數(shù)和熱膨脹系數(shù)；另一方面，彈性虛擬載荷和熱虛擬載荷在單元內(nèi)部節(jié)點不為0，所以二階和三階虛擬載荷除了反映基體和夾雜之間的相互作用，也刻畫了基體和夾雜內(nèi)部更加詳細的細觀信息。

3 數(shù)值分析

3.1 二維周期復(fù)合材料單胞結(jié)構(gòu)

結(jié)構(gòu)大小為15 mm×15 mm，單胞內(nèi)含有4塊夾雜，如圖6所示?；w的彈性模量和熱膨脹系數(shù)分別為E1=2×109Pa,a1=3×10-6/℃，夾雜的彈性模量和熱膨脹系數(shù)分別為E2=6×1010Pa,a2=1×10-6/℃，基體和夾雜的泊松比均為ν=0.2，結(jié)構(gòu)溫差T=16(x2+y2)。

分別使用式(12)計算得到均勻化彈性模量和均勻化熱膨脹系數(shù)分別為EH=2.67×109Pa,aH=2.5×10-6/℃。圖7給出了結(jié)構(gòu)沿A、B、C、D上節(jié)點在x1方向的位移曲線。表2給出了結(jié)構(gòu)總勢能泛函Π計算結(jié)果，其中FEM為有限元解，作為標準檢驗MHM的求解精度，MHM1表示包含一階攝動項的計算結(jié)果，MHM2表示包含二階攝動項的計算結(jié)果，MHM3表示包含三階攝動項的計算結(jié)果。結(jié)合圖7和表2可以得到：①包含一階攝動位移的計算結(jié)果MHM1誤差較大；②MHM2和MHM3的精度明顯高于MHM1；③MHM3精度高于MHM2。

圖6 二維周期復(fù)合材料單胞結(jié)構(gòu)Fig.6 Unit cell of 2D periodical composite structure

3.2 二維周期復(fù)合材料多胞結(jié)構(gòu)

結(jié)構(gòu)大小為45 mm×45 mm，包含5×5個單胞，每個單胞內(nèi)部含有一塊夾雜，如圖8所示，材料參數(shù)及溫差和3.1節(jié)相同。影響函數(shù)和宏觀位移導(dǎo)數(shù)分別使用超單胞周期邊界條件和微分求積有限單元法[30]計算得到，均勻化彈性模量和均勻化熱膨脹系數(shù)分別為EH=2.42×109Pa,aH=2.63×10-6/℃。表3給出了結(jié)構(gòu)的總勢能泛函計算結(jié)果，F(xiàn)EM、MHM1、MHM2、MHM3的含義與3.1節(jié)相同。圖9給出了結(jié)構(gòu)沿A′、B′、C′、D′上節(jié)點x1方向的位移曲線。結(jié)合圖9和表3可以得到：①包含一階攝動位移的計算結(jié)果MHM1誤差較大；②MHM2和MHM3的計算精度明顯好于MHM1；③MHM2和MHM3的計算精度較高，且兩者相差不大；④二階攝動項對計算精度的影響較大，而三階攝動項對計算精度的影響可以忽略。

圖7 沿縱線A、B、C、D上節(jié)點位移曲線Fig.7 Nodal displacement curves along longitudinal lines A,B,C and D

表2 二維周期復(fù)合材料單胞結(jié)構(gòu)的勢能泛函
Table 2 Potential energy function for 2D periodical
composites unit cell

攝動位移Π/(10-6J)MHMFEMΠFEM-ΠMHMΠFEM×100%MHM1-2.47515.8MHM2-2.916-2.9390.78MHM3-2.9370.068

圖8 二維周期復(fù)合材料多胞結(jié)構(gòu)Fig.8 2D periodical composite multi-cell structure

表3 二維周期復(fù)合材料多胞結(jié)構(gòu)的勢能泛函Table 3 Potential energy function for 2D periodical composite multi-cell structure

圖9 沿縱線A′、B′、C′、D′上節(jié)點位移曲線Fig.9 Nodal displacement curves along longitudinal lines A′，B′，C′ and D′

4 結(jié) 論

本文通過構(gòu)造攝動項的全解耦格式推導(dǎo)了熱力耦合問題高階MHM的數(shù)學(xué)表達式和有限元矩陣列式，將彈性影響函數(shù)和熱影響函數(shù)分別比擬為彈性虛擬位移和熱虛擬位移，利用虛擬載荷的自平衡性質(zhì)、量綱分析和空間矢量分布揭示了各階攝動項的物理意義；數(shù)值算例分析證明了公式推導(dǎo)結(jié)果和物理意義分析的正確性。通過熱力耦合問題MHM各階攝動項物理意義的分析和數(shù)值計算結(jié)果，給出了二階攝動項對于周期復(fù)合材料計算精度的必要性依據(jù)，為熱力耦合問題MHM的應(yīng)用奠定了力學(xué)基礎(chǔ)。