熱力耦合問題數學均勻化方法的物理意義

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

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

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

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

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

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

1 熱力耦合問題MHM的數學表達式和有限元矩陣列式

基于微觀結構的周期性和單胞域的一致性假設，均勻化理論可以解耦成為微觀結構異質邊界值的單胞問題和宏觀問題，對于結構線性小變形問題，熱力耦合的應力-應變關系為

(1)

總應變-位移和熱應變分別表示為

(2)

(3)

式中：uε為真實位移；amn為熱膨脹系數；T為結構溫差。

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

(4)

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

(5)

(6)

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

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

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

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

(7)

(8)

(9)

(10)

(11)

式中：

(12)

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

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

(13)

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

(14)

結構由于熱膨脹作用只產生線應變，剪應變為0，所以這種由于熱變形產生的應變可以看作結構的初始應變ε0。對于二維問題，ε0的表達式為

(15)

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

(16)

式中：f為外載荷。

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

在得到宏觀位移及其各階導數、各階溫差場導數及攝動位移后，結構真實細觀位移可以使用如下漸進展開表達式得到：

(17)

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

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

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

2.1 二維周期復合材料單胞結構

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

(18)

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

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

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

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

(19)

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

① 只要夾雜的材料性質相同，在每一塊夾雜上相同位置節點上的矢量都是相同的，這個結論普遍成立，與單胞包含夾雜的數量無關。

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

圖1 包含3塊夾雜的二維周期復合材料單胞結構Fig.1 2D periodical composites unit cell structures with 3 inclusions

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

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

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

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

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

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

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

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

2.2 周期復合材料桿單胞結構

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

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

圖5 包含一塊夾雜的周期復合材料桿單胞結構Fig.5 Periodical composites rod unit cell structure with one inclusion

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

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

(20)

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

(21)

(22)

(23)

分別求解方程可得到單胞結構彈性影響函數和熱影響函數為

(24)

(25)

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

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

(26)

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

(27)

二階彈性影響函數和熱影響函數分別為

(28)

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

(29)

三階熱影響函數和彈性影響函數分別為

(30)

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

3 數值分析

3.1 二維周期復合材料單胞結構

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

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

圖6 二維周期復合材料單胞結構Fig.6 Unit cell of 2D periodical composite structure

3.2 二維周期復合材料多胞結構

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

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

表2 二維周期復合材料單胞結構的勢能泛函
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 二維周期復合材料多胞結構Fig.8 2D periodical composite multi-cell structure

表3 二維周期復合材料多胞結構的勢能泛函Table 3 Potential energy function for 2D periodical composite multi-cell structure

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

4 結 論

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