基于等參梯度單元的功能梯度材料結構分析

張 龍，邱榮凱，程 俊，劉秉斌

（中國空氣動力研究與發展中心設備設計及測試技術研究所，綿陽 621000）

功能梯度材料可設計用于抵抗高溫度梯度、減小熱應力以及提高界面連接強度［1?2］，在航空航天熱結構領域得到了廣泛的應用研究，如高超飛行器的梯度密度一體化防隔熱材料［3］、沖壓發動機梯度隔熱層［4］和航空發動機渦輪葉片熱障涂層［5］等。近年來，隨著增材制造技術、化學氣相沉積等先進制造工藝不斷發展，使得功能梯度材料的可設計性越來越強、制備產業越來越成熟。與此同時，功能梯度材料特殊的力學行為也給結構分析與設計帶來了新的挑戰［6］。

伏培林等［7］采用解析法，開展了功能梯度懸臂梁的自由振動分析。針對歐拉?伯努利梁、瑞利梁和鐵木辛柯梁，采用分離變量法推導了功能梯度懸臂梁自由振動的解析解，進而給出了不同跨深比和材料性能梯度變化的懸臂梁前三階固有頻率。許琦等［8］提出了一種考慮橫法向熱應變的三參數Reddy 型高階功能梯度梁理論，可提高熱力響應的分析精度。李華東等［9］研究了四邊簡支具有功能梯度芯材的夾層板在分布載荷作用下的彎曲問題，基于Reissner 假設，根據功能梯度材料的本構方程得出了應力、位移及內力的表達式。高晟耀等［10］基于一階剪切變形理論推導了中等厚度功能梯度球環結構公式，利用里茲法求解得到功能梯度球環結構固有頻率，分析了一般邊界條件下中等厚度功能梯度球環結構的自由振動特性。高英山等［11］基于一階剪切變形假設和哈密頓原理開展了碳納米管增強功能梯度復合板非線性建模與仿真，討論了碳納米管體積分數、分布方式、結構寬厚比和載荷對功能梯度復合板的影響。蘇盛開等［12］采用經典歐拉梁理論和高階三角剪切變形理論，研究了多孔功能梯度梁的熱力耦合屈曲行為，討論了材料非均勻參數、孔隙率和長細比等參數對屈曲臨界溫度的影響。

以上研究均是基于解析或者半解析方法，研究對象為簡單的梁、板、殼等結構，受載形式也較為簡單。當結構或載荷復雜度增加時，有限元等數值計算方法是強有力的分析工具。但是由于功能梯度材料宏觀上的材料性質不均勻性，采用分層法［13］等常規有限元開展分析，勢必需要采用極高的網格密度以表征材料梯度，導致計算量巨大。Santare等［14］和Kim 等［15］提出了等參梯度有限元法，其單元本身可以反映材料力學性能梯度。Burlayenko等［16］和Xiong 等［17］分別在此基礎上提出了改進的梯度有限元法，可用于功能梯度材料的力?熱耦合分析和磁?電?彈性耦合分析。黃立新等［18］采用功能梯度材料平板，研究指出等參梯度單元相比常規單元具有更高的計算精度，并且隨著彈性模量梯度增大，梯度單元的優勢越大。韓國凱等［19］基于梯度有限元法思路，利用材料用戶子程序UMAT 定義材料梯度，開展了防隔熱一體化復合材料的整體性能優化設計。

本文發展了平面四節點與八節點等參梯度單元，從收斂性、計算精度和計算效率等方面驗證了梯度單元相對于常規單元的優勢。最后，基于ABAQUS 平臺開發等參梯度單元 UEL（User?defined element）子程序，開展了功能梯度材料開孔結構和懸臂梁的計算對比分析。

1 功能梯度材料模型

功能梯度材料由兩種或兩種以上材料混合制備而成，并且材料組分在空間上漸進均勻變化，不存在含量百分比階躍的材料界面，消除了傳統復合材料組分界面處力學性能突變這一缺陷，從而能更好地發揮復合材料的優勢。如圖1 所示，在x＜0區域內為材料A；在x＞L區域內為材料B；在0≤x≤L區域內為功能梯度材料，材料A 和材料B 的含量百分比漸進均勻變化。

圖1 功能梯度材料示意圖Fig.1 Sketch of functional graded materials

在x=0 邊界處，材料性能與材料A 相同；在x=L邊界處，材料性能與材料B 相同；在0＜x＜L區域內，根據組分含量百分比的變化程度不同，材料性能可為坐標x的線性函數或指數函數，即

式中：φ（x）、φA、φB分別為梯度材料、材料A、材料B的性能指標，可以是密度、彈性模量和泊松比等。

2 等參梯度單元

2.1 有限元基本理論方程

在有限元理論中，單元內任意一點處的位移可以表示為

式中：Ni為形函數；uei為編號為i的節點處位移；n為該單元的節點數量。

確定了單元位移后，利用幾何方程和物理方程，可將單元的應力應變表示為

式中：B為應變矩陣；D為材料剛度矩陣。

根據虛功原理，可推導有限單元剛度矩陣方程如下［20］

式中：fe為單元節點載荷；ke為單元剛度矩陣，可由式（7）得到

式中Ωe表示單元的控制域。將以上單元水平內的推導方程擴展到整個結構，即可得到總體的求解方程。

2.2 平面四邊形八節點與四節點等參單元

由于實際結構形狀復雜，實際單元可以為任意形狀。通過坐標變換，可以將總體坐標系xy下任意形狀的四邊形單元變換成局部坐標系ξη下的正方形單元。以平面八節點單元為例，如圖2所示。

圖2 坐標變化示意圖Fig.2 Sketch of coordinate change

采用與式（3）中相同的形函數Ni作為坐標變換插值函數，則單元幾何形狀和單元內場函數采用了相同數目的節點參數及相同的插值函數進行變換，即為等參變換，響應的單元被稱為等參單元。因此，坐標變換也可表示為

基于以上二維平面問題的等參變換，采用局部坐標系ξη表示式（7）中的單元剛度矩陣為

式中：t為平面單元厚度；J為坐標變換雅各比矩陣，可由式（11）得到

2.3 功能梯度材料的等參梯度單元

式（12）中的D為材料剛度矩陣，各項同性材料在平面應力狀態下的材料剛度矩陣為

式中E和ν分別為材料的彈性模量和泊松比，將E替換為E(1-ν2)、ν替換為ν(1-ν)即可得到平面應變狀態下的材料剛度矩陣。在單元內部，功能梯度材料的材料性能指標也由形函數Ni插值得到，即

將式（15～16）代入式（14）中得到

再將式（17）代入式（12），即可得到功能梯度材料等參梯度單元的單元剛度矩陣。

2.4 數值處理方法

本文采用高斯積分方法求解式（12）中的單元剛度矩陣，則

式中：m為每個積分方向上的積分點個數；ξi、ηj為積分點坐標；Wi、Wj為積分點處相應的權系數；f(ξi，ηj)為式（19）函數在積分點坐標ξi、ηj處的取值。

對于二維平面四邊形四節點和八節點單元，其內部高斯積分點的分布如圖3 所示。四節點單元內部共4 個高斯積分點，即每個積分方向2 個積分點；八節點單元內部共9 個高斯積分點，即每個積分方向3 個積分點。積分點坐標和權系數如表1 所示。

圖3 單元積分點Fig.3 Integration points of the element

表1 積分點坐標與權系數Table 1 Integration point coordinates and weight coeffi?cient

求得單元剛度矩陣ke后，可利用ABAQUS 進行二次開發功能，將等參梯度單元模型編寫為單元用戶子程序UEL，應用于功能梯度材料結構的數值模擬。單元用戶子程序UEL 的介紹及處理流程見參考文獻［18，21］。

2.5 后處理云圖顯示方法

由于UEL 單元只參與ABAQUS 計算，將其單元剛度矩陣用于集成總體剛度矩陣，而單元形狀、應力云圖無法直接在ABAQUS 的后處理中顯示。因此必須采用其他手段進行后處理查看云圖。常用的方法是將UEL 計算結果儲存到文本文件中，然后利用Matlab 等工具讀取文本數據繪制云圖。在將UEL 單元與ABAQUS 自帶單元的計算結果進行對比時，Matlab 等工具繪制云圖時采用的插值方法必須與ABAQUS 繪制云圖的插值方法一致，以避免插值方法不一致引入的誤差。而ABAQUS 后處理中不同類型單元所采用的插值方法不同，因此使用Matlab 等工具繪制云圖時也必須采用多種插值方法，過程比較復雜。

本文的后處理云圖顯示方法為：

（1）將UEL 單元計算結果儲存到文本文件中。

（2）將計算模型中的UEL 單元替換為與UEL節點、積分點一致的ABAQUS 自帶單元，得到新的計算模型。

（3）在新的計算模型中采用用戶子程序SIGI?NI，讀取文本文件中與該積分點坐標相匹配的UEL 單元計算結果。

（4）對新的計算模型施加“零載荷”，進行一次“假計算”，則后處理的顯示結果即為SIGINI 讀取的UEL 單元計算結果。

其中SIGINI 為ABAQUS 用戶子程序［21］，可用于定義計算模型的初始應力場。在本文中用于將UEL 單元的計算結果賦予與其節點、積分點一致的ABAQUS 自帶單元，然后進行“零載荷假計算”，從而生成的應力云圖即為UEL 單元的應力云圖。

3 單元驗證

3.1 模型說明與解析解

采用本文等參梯度單元對圖4 所示的功能梯度材料正方形平板開展有限元分析，驗證單元的收斂性，并與采用常規單元分層法［13］的計算結果進行對比，驗證本文等參梯度單元在計算精度和效率上的優勢。

圖4 計算模型示意圖（單位：mm）Fig.4 Sketch of computed model(Unit:mm)

正方形平板為材料A 和材料B 混合制備的功能梯度材料，邊長100 mm，厚度取單位厚度。考慮功能梯度材料的彈性模量如式（1）線性變化或式（2）指數變化的兩種情況，泊松比為常量0.3。令材料A 的彈性模量為100 GPa、材料B 的彈性模量為800 GPa。此處的材料性能數據用于本節等參梯度單元的驗證，不具有實際物理意義。

正方形平板受到的約束條件為在（0，0）點受x向和y向約束，在y=0 邊受y向約束。開展3 種載荷條件下的計算分析。

工況1：位移載荷條件，在y=100 邊施加位移Δy=1.0。

工況2：應力載荷條件，在y=100 邊施加均布拉力f=100 N/mm。

工況3：彎曲載荷條件，在y=100 邊施加線性分布力f（x）=100-2x（N/mm）；f（x）的正負分別表示指向y軸正向或負向。

當功能梯度材料的彈性模量線性變化時，3 種工況下y=0 坐標處y向應力的解析解分別如式（20～22）所示，推導過程見參考文獻［15］。

3.2 收斂性

采用本文的等參梯度單元計算工況1 下線性功能梯度材料、指數型功能梯度材料的應力分布分別如圖5、6 所示。圖中S為應力，S22表示y向應力分量。圖5 中平面四節點單元與平面八節點單元的計算結果采用的單元數目依次為1×1、3×3、6×6、10×10。從圖5 可以看出，隨著單元密度增加，應力峰值和云圖分布均保持穩定不變，表示已得到收斂結果。圖6 中平面四節點單元與平面八節點單元的計算結果采用的單元數目依次為1×1、3×3、10×10、30×30。從圖6 可以看出，隨著單元密度增加，應力峰值和云圖分布均發生變化，并最后趨于穩定，得到收斂結果。

圖5 線性功能梯度材料應力云圖（工況1）Fig.5 Stress contour of linearly functional graded material(Case 1)

圖6 指數型功能梯度材料應力云圖（工況1）Fig.6 Stress contour of exponentially functional graded ma?terial(Case 1)

從圖5、6 可以看出，對于線性功能梯度材料，無論是采用平面四節點還是八節點單元，均只需要采用1 個單元，即可得到收斂結果，而指數型功能梯度材料性能更復雜。采用平面四節點單元時，將10×10 網格進一步細化后，應力峰值變化量小于0.3%，在工程應用中可忽略，即采用平面四節點單元需要10×10 個單元得到收斂解。此外，采用平面八節點單元時，將3×3 網格進一步細化后，應力峰值變化量小于0.2%，在工程應用中可忽略，即采用平面四節點單元僅需要3×3 個單元得到收斂解。圖6（e）中四節點單元10×10 網格的節點數量為11×11=121，圖6（d）中八節點單元3×3 網格的節點數量為7×7=49，而圖6（d）的計算精度高于圖6（e），可見相比于四節點梯度單元，八節點梯度單元可以采用更少數量的單元或節點，而得到更高精度的計算結果。

將本文的等參梯度單元用于工況2 和工況3下功能梯度材料的計算，亦得到收斂解。為簡潔起見，工況2 和工況3 的結果不再詳細展示，僅給出單元數目為10×10 的計算結果，在3.3 節中可查看。

從以上分析可以看出，采用本文的平面四節點或八節點等參梯度單元均可得到功能梯度材料的收斂解，但八節點單元比四節點單元需要更少的單元數/節點數即可得到收斂結果。

3.3 計算精度和效率

本節給出了等參梯度單元與常規單元的計算結果，并與解析結果進行對比，驗證本文方法在計算精度和效率上的優勢。等參梯度單元內部的材料呈梯度特性，而常規單元內部材料為均勻材料。采用常規單元模擬功能材料，需要采用分層法［13］，對于不同層的單元賦予不同的材料屬性以表征材料梯度。無論采用平面八節點還是四節點梯度單元，計算精度均優于平面八節點常規單元，更優于平面四節點常規單元。為簡潔起見，本節只給出10×10 平面八節點梯度單元與常規單元的計算結果對比。本文采用的常規單元均為ABAQUS 自帶的全積分單元，八節點常規單元類型為CPS8，四節點常規單元類型為CPS4。

3.3.1 工況1 計算結果對比

工況1 下線性和指數型功能梯度材料的應力分布分別如圖7、8 所示，提取坐標y=0 上的應力值，與式（20）、式（23）的解析計算結果對比分別如圖9、10 所示。從圖中可以看出，梯度單元的計算結果與解析解吻合較好，而常規單元與解析解誤差較大。對于線性功能梯度材料，10×10 常規單元（圖7（b））的計算精度低于1×1 四節點或八節點梯度單元（圖5（a）或（b））的計算精度。而對于指數型功能梯度材料，10×10 常規單元（圖8（b））的計算精度低于1×1 八節點梯度單元（圖6（b））的計算精度；也低于3×3 四節點梯度單元（圖6（c））的計算精度。

圖7 線性功能梯度材料等梯度單元與常規單元應力云圖（工況1）Fig.7 Stress contour of graded element and conventional el?ement of linearly functional graded material (Case 1)

圖8 指數型功能梯度材料等梯度單元與常規單元應力云圖（工況1）Fig.8 Stress contour of graded element and conventional el?ement of exponentially functional graded material(Case 1)

圖9 坐標y=0 處線性功能梯度材料應力云圖（工況1）Fig.9 Stress value of linearly functional graded material at y=0 (Case 1)

圖10 坐標y=0 處指數型功能梯度材料應力云圖（工況1）Fig.10 Stress value of exponentially functional graded ma?terial at y=0 (Case 1)

3.3.2 工況2 計算結果對比

工況2 下線性和指數型功能梯度材料的應力分布分別如圖11～12 所示，提取坐標y=0 上的應力值，與式（21）、式（24）的解析計算結果對比分別如圖13～14 所示。從圖中可以看出，梯度單元的計算結果與解析解吻合較好，而常規單元與解析解誤差較大。

圖11 線性功能梯度材料等梯度單元與常規單元應力云圖（工況2）Fig.11 Stress contour of graded element and conventional element of linearly functional graded material(Case 2)

圖12 指數型功能梯度材料等梯度單元與常規單元應力云圖（工況2）Fig.12 Stress contour of graded element and conventional element of exponentially functional graded material(Case 2)

圖13 坐標y=0 處線性功能梯度材料應力云圖（工況2）Fig.13 Stress value of linearly functional graded material at y=0 (Case 2)

3.3.3 工況3 計算結果對比

圖14 坐標y=0 處指數型功能梯度材料應力云圖（工況2）Fig.14 Stress value of exponentially functional graded ma?terial at y=0 (Case 2)

工況3 下線性和指數型功能梯度材料的應力分布分別如圖15～16 所示，提取坐標y=0 上的應力值，與式（22）、式（25）的解析計算結果對比分別如圖17～18 所示。從圖中可以看出，梯度單元的計算結果與解析解吻合較好，而常規單元與解析解誤差較大。

圖15 線性功能梯度材料等梯度單元與常規單元應力云圖（工況3）Fig.15 Stress contour of graded element and conventional element of linearly functional graded material(Case 3)

圖16 指數型功能梯度材料等梯度單元與常規單元應力云圖（工況3）Fig.16 Stress contour of graded element and conventional element of exponentially functional graded material(Case 3)

圖17 坐標y=0 處線性功能梯度材料應力云圖（工況3）Fig.17 Stress vaule of linearly functional graded material at y=0 (Case 3)

圖18 坐標y=0 處指數型功能梯度材料應力云圖（工況3）Fig.18 Stress vaule of exponentially functional graded ma?terial at y=0 (Case 3)

根據工況1～3 的計算結果對比可以看出，常規單元由于采用分層法，每層單元的材料屬性不同，邊界處材料性能發生階躍，計算的應力場也在邊界處也產生階躍，與實際情況不符，特別是當載荷條件較復雜時，應力場嚴重失真，如圖11（b）與圖12（b）所示。而梯度單元將材料梯度定義在了單元內部，避免了單元邊界處的材料性能階躍，得到的應力場光滑連續。梯度單元可以采用較少數目的單元，即可得到比常規單元更為精確的計算結果，無論是在計算精度和效率上均優于常規單元。

4 結構分析算例

4.1 開孔結構拉伸應力分析算例1

4.1.1 模型說明

本節算例模型如圖19 所示［15］，為單位厚度的功能梯度材料帶孔平板結構件，端部受均布拉伸載荷q作用。該平板的基本組分材料為TiB 和Ti，TiB 的彈性模量和泊松比分別為375 GPa 和0.14，Ti 的彈性模量和泊松比分別為107 GPa 和0.34。該結構件材料性能分布關于x軸對稱、沿y軸方向呈指數型變化，如式（26～27）所示。

圖19 帶孔結構拉伸示意圖（單位：mm）Fig.19 Sketch of perforated structure stretching (Unit:mm)

式中W為試驗件寬度，其值為14 mm。

4.1.2 結果對比

采用八節點梯度單元和常規單元計算的應力結果分別如圖20～21 所示，圖中單元特征尺寸為該計算模型中最小單元面積的平方根，用于表征網格密度。圖中S為應力，S11表示x向應力分量。單元特征尺寸越大，網格密度越小，反之網格密度越大。

圖20（a～c）單元密度依次增大，應力結果收斂較快。當單元特征尺寸由 0.531 mm 變為0.266 mm 時，應力峰值變化量小于0.05%，在工程應用中可忽略，即八節點梯度單元采用單元特征尺寸為0.531 mm 的計算模型（圖20（b））即可得到收斂解。

圖20 8 節點梯度單元應力計算結果（算列1）Fig.20 Computed stress results of 8?node graded element(Example 1)

從圖21（a～d）可以看出，由于采用了分層單元法，單元邊界處材料性能發生階躍，常規單元計算的應力場呈現明顯的鋸齒狀，局部應力場失真。當采用常規單元時，需要更密的網格才能得到收斂解，計算精度和效率大大低于梯度單元。此外，當網格密度較小時，計算的最大應力位置位于左下方內圓弧上，如圖21（a）所示。而隨著網格密度增大，計算的最大應力位置轉移到由上方外圓弧上，如圖21（c，d）所示。可以看出，當網格密度不夠時，計算得到的最大應力位置和數值均與實際不符。

圖21 8 節點常規單元應力計算結果（算列1）Fig.21 Computed stress results of 8?node conventional ele?ment(Example 1)

4.2 懸臂梁彎曲應力分析算例2

4.2.1 模型說明

本節算例模型如圖22 所示，為單位厚度的功能梯度材料懸臂梁，端部受集中力作用而發生彎曲變形。該懸臂梁的基本組分材料同算例1，料性能分布關于x軸對稱、沿y軸方向呈指數型變化，如式（26～27）所示。本算例中W′為懸臂梁高度，其值為10 mm。

圖22 懸臂梁彎曲示意圖（單位：mm）Fig.22 Sketch of cantilever bending (Unit:mm)

4.2.2 結果對比

采用八節點梯度單元和常規單元計算的應力結果分別如圖23～24 所示。圖23 和圖24 中，（a～c）單元密度依次增大，采用的單元數目依次為2×6、4×12、10×30。隨著單元密度增大，圖23 應力結果收斂較快。當單元數目由4×12 變為10×30 時，應力峰值變化量小于1.1%，在工程應用中可忽略，即八節點梯度單元采用單元數目為4×12的計算模型（圖23（b））即可得到收斂解。而圖24隨著單元密度增大應力結果收斂較慢，當單元數目由4×12 變為10×30 時，應力峰值變化量大于12.7%，計算結果尚未收斂。由以上分析可以看出，當采用常規單元時，需要更密的網格才能得到收斂解，計算精度和效率大大低于梯度單元。

圖23 8 節點梯度單元應力計算結果（算列2）Fig.23 Computed stress results of 8?node graded element(Example 2)

圖24 8 節點常規單元應力計算結果（算列2）Fig.24 Computed stress results of 8?node conventional ele?ment(Example 2)

此外，從圖24（a～c）還可以看出，由于分層法單元邊界處材料性能階躍，當網格密度較小時，常規單元計算的應力場呈現鋸齒狀（單元數目2×6、4×12），直到網格密度增大到10×30，應力場的鋸齒狀才逐步消失。而圖23 采用梯度單元，當單元數目為2×6 時，所計算得到的應力場也比較光滑。將圖23 和圖24 中相同單元數目的計算模型對比可以看出，梯度單元的計算精度大大高于常規單元。

5 結 論

（1）將本文的平面四節點或八節點等參梯度單元用于功能梯度材料在位移載荷、應力載荷、彎曲載荷3 種工況下的計算，均可得到收斂解，但八節點單元比四節點單元需要更少的單元數/節點數即可得到收斂結果。

（2）將本文的梯度單元與常規單元對比，相同單元數目下，無論平面八節點還是四節點梯度單元，計算精度均優于平面八節點常規單元，更優于平面四節點常規單元。

（3）本文的梯度單元將材料梯度定義在單元內部，避免了單元邊界處的材料性能階躍，得到的應力場光滑連續，可以采用較少數目的單元，即可得到比常規單元更為精確的計算結果。

（4）將本文的梯度單元用于結構計算，需要較少數目的單元即可得到收斂解，計算精度和效率大大優于常規單元。