基于等參梯度單元的功能梯度材料結(jié)構(gòu)分析

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

（中國空氣動力研究與發(fā)展中心設(shè)備設(shè)計及測試技術(shù)研究所，綿陽 621000）

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

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

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

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

1 功能梯度材料模型

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

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

在x=0 邊界處，材料性能與材料A 相同；在x=L邊界處，材料性能與材料B 相同；在0＜x＜L區(qū)域內(nèi)，根據(jù)組分含量百分比的變化程度不同，材料性能可為坐標(biāo)x的線性函數(shù)或指數(shù)函數(shù)，即

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

2 等參梯度單元

2.1 有限元基本理論方程

在有限元理論中，單元內(nèi)任意一點(diǎn)處的位移可以表示為

式中：Ni為形函數(shù)；uei為編號為i的節(jié)點(diǎn)處位移；n為該單元的節(jié)點(diǎn)數(shù)量。

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

式中：B為應(yīng)變矩陣；D為材料剛度矩陣。

根據(jù)虛功原理，可推導(dǎo)有限單元剛度矩陣方程如下［20］

式中：fe為單元節(jié)點(diǎn)載荷；ke為單元剛度矩陣，可由式（7）得到

式中Ωe表示單元的控制域。將以上單元水平內(nèi)的推導(dǎo)方程擴(kuò)展到整個結(jié)構(gòu)，即可得到總體的求解方程。

2.2 平面四邊形八節(jié)點(diǎn)與四節(jié)點(diǎn)等參單元

由于實(shí)際結(jié)構(gòu)形狀復(fù)雜，實(shí)際單元可以為任意形狀。通過坐標(biāo)變換，可以將總體坐標(biāo)系xy下任意形狀的四邊形單元變換成局部坐標(biāo)系ξη下的正方形單元。以平面八節(jié)點(diǎn)單元為例，如圖2所示。

圖2 坐標(biāo)變化示意圖Fig.2 Sketch of coordinate change

采用與式（3）中相同的形函數(shù)Ni作為坐標(biāo)變換插值函數(shù)，則單元幾何形狀和單元內(nèi)場函數(shù)采用了相同數(shù)目的節(jié)點(diǎn)參數(shù)及相同的插值函數(shù)進(jìn)行變換，即為等參變換，響應(yīng)的單元被稱為等參單元。因此，坐標(biāo)變換也可表示為

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

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

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

式（12）中的D為材料剛度矩陣，各項(xiàng)同性材料在平面應(yīng)力狀態(tài)下的材料剛度矩陣為

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

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

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

2.4 數(shù)值處理方法

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

式中：m為每個積分方向上的積分點(diǎn)個數(shù)；ξi、ηj為積分點(diǎn)坐標(biāo)；Wi、Wj為積分點(diǎn)處相應(yīng)的權(quán)系數(shù)；f(ξi，ηj)為式（19）函數(shù)在積分點(diǎn)坐標(biāo)ξi、ηj處的取值。

對于二維平面四邊形四節(jié)點(diǎn)和八節(jié)點(diǎn)單元，其內(nèi)部高斯積分點(diǎn)的分布如圖3 所示。四節(jié)點(diǎn)單元內(nèi)部共4 個高斯積分點(diǎn)，即每個積分方向2 個積分點(diǎn)；八節(jié)點(diǎn)單元內(nèi)部共9 個高斯積分點(diǎn)，即每個積分方向3 個積分點(diǎn)。積分點(diǎn)坐標(biāo)和權(quán)系數(shù)如表1 所示。

圖3 單元積分點(diǎn)Fig.3 Integration points of the element

表1 積分點(diǎn)坐標(biāo)與權(quán)系數(shù)Table 1 Integration point coordinates and weight coeffi?cient

求得單元剛度矩陣ke后，可利用ABAQUS 進(jìn)行二次開發(fā)功能，將等參梯度單元模型編寫為單元用戶子程序UEL，應(yīng)用于功能梯度材料結(jié)構(gòu)的數(shù)值模擬。單元用戶子程序UEL 的介紹及處理流程見參考文獻(xiàn)［18，21］。

2.5 后處理云圖顯示方法

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

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

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

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

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

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

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

3 單元驗(yàn)證

3.1 模型說明與解析解

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

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

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

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

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

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

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

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

3.2 收斂性

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

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

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

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

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

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

3.3 計算精度和效率

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

3.3.1 工況1 計算結(jié)果對比

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

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

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

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

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

3.3.2 工況2 計算結(jié)果對比

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

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

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

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

3.3.3 工況3 計算結(jié)果對比

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

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

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

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

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

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

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

4 結(jié)構(gòu)分析算例

4.1 開孔結(jié)構(gòu)拉伸應(yīng)力分析算例1

4.1.1 模型說明

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

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

式中W為試驗(yàn)件寬度，其值為14 mm。

4.1.2 結(jié)果對比

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

圖20（a～c）單元密度依次增大，應(yīng)力結(jié)果收斂較快。當(dāng)單元特征尺寸由 0.531 mm 變?yōu)?.266 mm 時，應(yīng)力峰值變化量小于0.05%，在工程應(yīng)用中可忽略，即八節(jié)點(diǎn)梯度單元采用單元特征尺寸為0.531 mm 的計算模型（圖20（b））即可得到收斂解。

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

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

圖21 8 節(jié)點(diǎn)常規(guī)單元應(yīng)力計算結(jié)果（算列1）Fig.21 Computed stress results of 8?node conventional ele?ment(Example 1)

4.2 懸臂梁彎曲應(yīng)力分析算例2

4.2.1 模型說明

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

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

4.2.2 結(jié)果對比

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

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

圖24 8 節(jié)點(diǎn)常規(guī)單元應(yīng)力計算結(jié)果（算列2）Fig.24 Computed stress results of 8?node conventional ele?ment(Example 2)

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

5 結(jié) 論

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

（2）將本文的梯度單元與常規(guī)單元對比，相同單元數(shù)目下，無論平面八節(jié)點(diǎn)還是四節(jié)點(diǎn)梯度單元，計算精度均優(yōu)于平面八節(jié)點(diǎn)常規(guī)單元，更優(yōu)于平面四節(jié)點(diǎn)常規(guī)單元。

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

（4）將本文的梯度單元用于結(jié)構(gòu)計算，需要較少數(shù)目的單元即可得到收斂解，計算精度和效率大大優(yōu)于常規(guī)單元。