方向性熱膨脹系數可設計點陣結構均勻化性能分析

王 帥，牛 飛，牛 斌

(1.大連理工大學機械工程學院，遼寧 大連 116024;2.中國運載火箭技術研究院，北京 100076)

0 引言

隨著點陣結構的研究和發展，點陣結構的單胞數量越來越多，構型越來越復雜，對此類結構的分析往往很耗時甚至無法計算分析。通常的做法是采用近似模型，將宏觀結構中具有非均質特性的點陣結構單胞進行等效處理，以此來簡化宏觀模型的分析。常見的等效性能預測方法有：以變分原理為基礎的定界法，如Hashin[1]通過解析法得到材料等效參數的上限和下限范圍。基于夾雜理論的近似模型：自洽模型、廣義自洽模型、Mori-Tanaka(M-T)模型等，以及基于此構建的數值方法，如自洽法[2]、廣義自洽法[3]、M-T法[4]。考慮到實際計算模型的復雜性和預測方法的通用性，數值求解方法在預測結構等效性能方面比較流行，主要有代表體元法(RVE)和漸近均勻化方法(AH)。

代表體元法(RVE)[5-6]基于能量等效原理，選取周期性點陣結構(代表體元)作為研究對象，通過施加單位位移邊界條件(Dirichlet邊界)或者力邊界條件(Neumann邊界)，使得點陣結構應變能和相應的均質材料單胞應變能等效，從而求出等效性能。漸近均勻化(AH)[7～9]方法是另外一種預測周期性點陣結構等效性質的數值方法。漸近均勻化方法基于攝動理論和虛位移原理，具有嚴謹的數學理論推導，當單胞尺寸相對于宏觀結構尺寸非常小時，能夠給出精確的預測結果。近年來，許多學者將漸近均勻化與商業軟件結合，保留嚴格的數學理論邏輯，簡化程序開發流程，實現快速簡單計算，如程耿東等[10]提出一種新漸近均勻化實現方法，可以使用商業軟件輕松實現多種點陣結構等效彈性模量的計算，張永存等[11]提出了基于漸進均勻化方法預測周期性復合材料等效熱膨脹系數的新算法。

本文采用漸近均勻化方法獲得點陣單胞的等效彈性模量和熱膨脹系數，比較點陣結構等效分析與離散分析之間的差異。首先，本文針對一種單方向熱膨脹系數可設計的點陣單胞，用漸近均勻化方法，獲得點陣單胞的等效模量和熱膨脹系數屬性。然后，分別建立一個點陣結構離散模型和采用均勻化后等效模型，討論在不同幾何外形、邊界條件、載荷情況和點陣數目下均勻化方法預測的準確性。

1 漸近均勻化理論

對于周期性點陣結構，漸近均勻化方法可以預測其等效性能。點陣結構的等效彈性模量可以用應變能的形式表示為[10]。

(1)

其中，Y是點陣單胞的定義域，|Y|是點陣單胞的體積，D為彈性矩陣，ε0(kl)為單位應變，ε*(kl)為特征應變。

式(1)可以變形為

(2)

其中，x0(ij)是初始位移場，x*(ij)是特征位移場，f(kl)是初始位移場對應的節點反力，f*(kl)是特征位移場對應的節點反力。以二維為例，特征位移場x0(ij)可以表示為[10]:

(3)

f(kl)和f*(kl)可以分別表示為

f(kl)=Kx0(kl)

(4)

f*(kl)=Kx*(kl)

(5)

所以，只要獲得f(kl)和f*(kl)就可以計算有效彈性模量[10]。建立一個單胞有限元模型，對模型施加一個初始位移場x0(kl)進行一次有限元分析，提取節點的反力，即f(kl)。再將獲得的f(kl)施加在模型上，施加周期性邊界條件進行一次有限元分析，可以獲得特征位移場x*(kl)。最后，在模型上施加特征位移場x*(kl)進行一次有限元分析，獲得節點的反力，即f*(kl)。將獲得的數據代入到公式(2)中即可計算出有效彈性模量。

等效熱膨脹系數在等效彈性模量的基礎上進行計算。等效熱膨脹系數與等效彈性模量EH和熱彈性常數β有關[11]，關系如公式(6)所示，所以確定熱膨脹系數還需計算熱彈性常數β。

α=(EH)-1β

(6)

熱彈性常數β的表達式如式(7)所示：

(7)

其中，Rα為單位負溫升下的節點反力，Λ為特征位移場，RΛ為在特征位移場Λ下的節點反力。

計算等效熱膨脹系數的步驟如下[11]：首先，約束所有節點位移，施加單位負溫升，獲得節點反力Rα。然后將節點反力Rα施加在單胞的每一個節點上，施加周期性邊界條件，通過一次靜力分析，獲得特征位移場Λ。再將各節點施加特征位移場Λ，獲得節點反力RΛ。最終，將獲得的節點反力Rα和RΛ代入到公式(7)中計算出熱彈性常數β，再由公式(6)獲得單胞的等效熱膨脹系數。

2 點陣模型等效屬性計算

本文使用的點陣結構單胞如圖1(a)所示，在雙材料金字塔結構[12]的基礎上，將底面改為板結構，提高承載性能，由兩個雙材料板桿金字塔構型對接形成本文點陣結構單胞。其中圖1(a)、圖1(b)中板材為材料1，桿件為材料2。

圖1 點陣結構單胞圖

參考文獻[13]中的推導過程，金字塔結構各邊與高度h、角度β之間滿足關系

(8)

分別對L1、L2求微分可得：

(9)

其中，dL1、dL2、dh、dβ分別是L1、L2、h、β的微分形式。

根據熱膨脹系數的定義，圖1(c)中淺色桿件(材料1)的熱膨脹系數α1=dL1/L1dt，深色桿件(材料1)的熱膨脹系數α2=dL2/L2dt，公式(9)變形為：

(10)

其中，dh/hdt為四棱錐沿高度h方向的熱膨脹系數αh(后面簡寫為等效熱膨脹系數)，對公式(10)約分，再將α1、α2、αh代入到公式(10)中，得到公式(11)。

(11)

(12)

從公式(12)中可以看出，選定兩種不同的材料或設計角度β的大小，都可以改變等效熱膨脹系數αh的大小，實現正、零和負熱膨脹系數的變化。

當固定角度β，選擇不同的材料組合來改變αh的數值。點陣結構單胞中材料1選擇為1Cr8Ni9，材料2選擇為4J33，材料屬性如表1中所示。將表1中的熱膨脹系數代入到公式(12)，當角β=45°，點陣結構單胞等效熱膨脹系數為0。然后固定角度β為45°不變，選擇不同的材料組合，獲得不同等效熱膨脹系數的點陣結構單胞。采用漸近均勻化方法，計算點陣結構單胞的等效材料屬性，結果如表2所示。表2中Ex、Ey、Ex表示點陣結構單胞XYZ三個方向的彈性模量，αx、αy、αz表示點陣結構單胞XYZ三個方向的熱膨脹系數。

表1 材料屬性[11]

表2 單胞等效材料屬性

由表2中第一行數據，點陣結構單胞在XY方向上的熱膨脹系數與材料1的熱膨脹系數相近，在Z方向上的熱膨脹系數比XY方向上的熱膨脹系數小兩個數量級，符合低熱膨脹的設計。對比3個點陣結構單胞的等效彈性模量可以看出，單胞在XY方向上的彈性模量相同，Z方向上的彈性模量較低。對比3個點陣結構單胞的等效熱膨脹系數可以看出，點陣結構單胞沿XY兩方向的熱膨脹系數與板材料(材料1)的熱膨脹系數相同，Z方向上的熱膨脹系數是可以設計的，可以實現由正到低再到負熱膨脹系數的變化。

3 點陣結構單胞均勻化性能驗證

點陣結構單胞均勻化之后可以看作是一個各向異性的等效材料。先采用簡單幾何的點陣結構為例，分別建立均勻化后的等效模型和點陣離散梁/板模型，考慮不同邊界條件、載荷情況和單胞數目，驗證漸近均勻化方法預測的等效彈性模量和熱膨脹系數的準確性。再使用復雜拓撲的點陣結構，同樣建立均勻化后的實體模型和點陣離散模型，探究幾何外形對漸近均勻化方法預測結果準確性的影響。

3.1 不同單胞數目

采用圖1(a)中的點陣結構單胞，填充一個10 mm×10 mm×14.14 mm的長方體區域，采用不同大小的點陣結構單胞進行填充，分別形成2×2×2、5×5×5和10×10×10點陣結構，建立點陣離散梁/板模型，如圖2所示。隨著單胞數的增多，桿徑和板厚相應變小，保證單胞的體分比不變。基于點陣結構均勻化等效性能，建立均勻化后的等效模型。載荷及邊界條件為：底面固定，頂面受到一個4 MPa的均布載荷，位移計算數據如表3所示。表3中數據為頂面上的所有節點位移的平均值。

圖2 點陣填充結構

表3 平壓下節點位移平均值

由表3中數據可以看出，隨著點陣結構單胞數目的增多，均勻化獲得的等效性能越來越能反映實際情況。當單胞數量為10×10×10時，離散分析和均勻化等效分析在主方向上的誤差在1%左右。之后的仿真分析就以10×10×10的點陣結構離散模型和均勻化后等效模型展開。

3.2 不同點陣邊界條件

均勻化后的等效結構不論如何截取都是同樣材料屬性的實體，但點陣結構卻不同。采用圖3(a)所示的3種方式截取點陣結構作為3種不同的單胞，再以此3種單胞陣列成周期性結構，如圖3(b)、(c)、(d)所示，這樣形成的結構在中心大部分區域并無差異，但在邊界上卻不同。通過仿真分析，比較不同截取方式對點陣性能的影響。位移計算數據在表4中列出。表4中數據為頂面上的所有節點位移的平均值。

圖3 點陣截取方法及陣列結構

表4 平壓下節點位移平均值

橫向對比表4中數據，發現基本滿足上一節的結論，即點陣單胞數量越多，均勻化結果預測越準確。再縱向對比表4中數據發現，截取方式1和截取方式3，對點陣性能的影響較小，但截取方式2的離散計算與均勻化預測結果誤差大。原因是，截取方式2會產生不完整點陣單胞，極大的削弱點陣結構的剛度。

3.3 不同載荷情況

改變模型的載荷及邊界條件，探究不同載荷情況下均勻化方法預測的準確性。將3.1節中計算的載荷情況視為第一種載荷條件。第二種載荷條件：底面固定，頂面上一條邊線上受到一個均布載荷。因為邊界效應的影響，比較位移的時候不考慮施加載荷處的位移情況，比較下一層點陣結構和等效模型對應位置的位移。因為載荷的不對稱性，同一個面上位移響應也是不同的，所以不再對比面上的平均位移，而是比較線上平均位移，具體位置如圖4(a)所示。線上平均位移計算結果如表5所示。

圖4 載荷及邊界條件

對比表5中的數據，均勻化后實體模型與點陣離散模型的位移基本吻合，誤差在7%左右，相比于平壓情況誤差有所增加。在一些位移較小的地方，如表5中的X方向位移，相對誤差較大，但絕對誤差數量級一致。與表3中數據對比發現，載荷條件越復雜，漸近均勻化預測的誤差越大。

表5 側壓載荷下節點位移平均值

第三種載荷條件：施加一個400 ℃均勻溫度載荷。為了模擬一個自由膨脹的過程，邊界條件設置為約束底面和底面上兩邊的形式。如圖4(b)所示，約束底面Z方向的自由度，約束底面上一條X方向邊線Y方向自由度，約束底面上一條Y方向邊線X方向自由度。均勻溫度載荷下節點位移計算結果如表6。表6中數據為點陣結構頂面上所有節點位移平均值。

表6 均勻溫度載荷下節點位移平均值

從表6中數據發現，兩者的相對誤差很大。但是，因為使用的點陣結構單胞在Z方向上的理論熱膨脹系數為0，Z方向位移很小。將點陣結構單胞的兩種初始材料修改為AL-5A02和1Cr8Ni9，材料屬性如表1所示，獲得一個新的等效材料屬性，如表7，并進行如圖4(b)相同的仿真分析，均勻溫度載荷下節點位移計算結果如表8。

表7 單胞等效材料屬性

表8 均勻溫度載荷下節點位移平均值

由表8中數據可以看出，均勻化后的等效模型與點陣離散模型位移吻合良好。比較表6和表8中的數據可以發現，兩次的絕對誤差均在10-3左右。可以看出，表6中的結果因為本身熱膨脹系數較小，導致的結果相對誤差較大。

3.4 不同幾何外形

本文使用的點陣結構單胞沿高度方向的熱膨脹系數具有可設計性，采用表2中3種單胞，填充一簡單拓撲結構，如圖5所示。分別創建均勻化后的等效模型和點陣結構離散模型，約束底面Z方向自由度、后面Y方向自由度和左端面X方向自由度，施加溫度載荷，比較兩者的位移響應。位移計算結果如表9所示。表9中數據為點陣結構頂面上所有節點位移平均值。

圖5 點陣填充拓撲結構

表9中第一行數據中，實體模型與離散模型Z向位移相對誤差很大，與表6的情況相同，都是因為點陣結構Z方向的熱膨脹系數為0。通過對比表9中數據，均勻化后的實體模型和點陣結構離散模型位移吻合良好。通過比較Z向位移，發現不同點陣結構可以實現結構單方向正、負熱膨脹系數的設計。

表9 均勻溫度載荷下節點位移平均值

4 結論

以雙材料四棱錐結構為基礎，針對一種單方向熱膨脹系數可設計的點陣結構單胞，采用漸近均勻化方法計算點陣單胞的等效彈性模量和熱膨脹系數。基于此，針對不同的點陣填充結構分別創建均勻化后等效模型和離散模型，驗證在不同幾何、不同邊界和載荷、不同單胞數目下均勻化方法預測的等效性能的準確性。獲得如下結論：

1)通過設計不同材料組合和幾何構型，可以實現點陣結構熱膨脹系數由正到負的變化，并采用漸近均勻化方法驗證等效熱膨脹系數的變化；

2)基于漸近均勻化方法計算的性能與點陣離散結構的吻合情況受點陣單胞數量的影響，當點陣單胞個數足夠多時，在力載荷和溫度載荷的作用下基于均勻化和離散計算的位移吻合結果令人滿意。隨著邊界和載荷條件的復雜，相比于點陣離散分析結果，熱彈性漸近均勻化的計算結果誤差變大。

3)采用不同截取方式形成的單胞陣列而成的結構，中心大部分并無不同，但結構的邊界卻有所差異，在邊界上存在不完整點陣單胞，降低點陣結構的等效剛度，使均勻化實體模型與離散模型在邊界上的位移計算結果相差很大。