尹碩輝，余天堂，劉 鵬

等幾何有限元法是一種將CAD與CAE有機(jī)統(tǒng)一起來的新型有限元法［1］，該方法采用CAD中樣條基函數(shù)(如B樣條、T樣條、NURBS樣條等)作為形函數(shù)，從而能形成“精確”的計(jì)算網(wǎng)格。與常規(guī)有限單元法相比，等幾何有限元法具有高精度和高收斂率的優(yōu)點(diǎn)。此外，常規(guī)有限單元法近似通常僅具有C0連續(xù)性，不能有效地求解高階偏微分方程問題(如薄板殼等)。等幾何有限元法可以構(gòu)造任意高階連續(xù)的基函數(shù)，能有效地求解高階偏微分方程。等幾何有限元法已成功地求解流體力學(xué)［2］、板殼分析［3］、電磁場［4］、相場［5］、拓?fù)鋬?yōu)化［6］以及裂紋擴(kuò)展［7-8］等問題。

功能梯度材料(FGM)是一種材料性質(zhì)和功能呈梯度變化的新材料［9］，使結(jié)構(gòu)既具有良好的耐高溫和抗腐蝕性能又具有較高的強(qiáng)度和韌性，因此在航空航天、核能、生物醫(yī)學(xué)、機(jī)械、建筑等領(lǐng)域得到了廣泛應(yīng)用。功能梯度材料結(jié)構(gòu)在工作狀態(tài)下表現(xiàn)出復(fù)雜的力學(xué)行為。一些學(xué)者基于不同的板理論(如經(jīng)典板理論［10］、一階剪切變形板理論［11］和高階板理論［12］等)和數(shù)值方法(如有限元法［11］、無網(wǎng)格法［13-14］和光滑有限元法［15-16］等)分析了功能梯度板的力學(xué)行為，得到了一些有益的結(jié)論。

本文針對金屬-陶瓷功能梯度板，考慮板材料性能沿厚度方向呈梯度連續(xù)變化，基于一階剪切變形板理論，推導(dǎo)了非均勻有理B樣條(NURBS)等幾何有限元基本公式;采用NURBS等幾何有限元分析了梯度指數(shù)、邊界條件及長厚比對功能梯度板自振頻率的影響。

1 基于一階剪切變形理論的功能梯度板的基本方程

圖1所示為金屬-陶瓷功能梯度非均勻材料板，其材料性質(zhì)沿厚度方向呈現(xiàn)梯度變化。設(shè)從板的上表面到下表面，金屬的體積分?jǐn)?shù)從100%過渡到0%，陶瓷的體積分?jǐn)?shù)則由0%到100%。

基于一階剪切板理論［11］，板內(nèi)任意一點(diǎn)的位移為:

其中:u0、v0和 w0分別代表板中面某點(diǎn)的位移;βx、βy代表中面法線變形后在xz平面和yz平面上的轉(zhuǎn)角。

圖1 金屬-陶瓷功能梯度板示意圖Fig.1 The graph of metal/ceramic FGM plate

平面內(nèi)的應(yīng)變?yōu)?

其中:εm和εb分別為中面面內(nèi)應(yīng)變場和板的曲率場，

此外，板的橫向剪應(yīng)變場為:

根據(jù)胡克定律，在正軸下板內(nèi)的應(yīng)力-應(yīng)變關(guān)系可表示為:

其中:σ=［σxσyτxy］T，τ=［τxzτyz］T分別為平面內(nèi)應(yīng)力及橫向應(yīng)力。彈性矩陣為:

其中α為剪切修正因子(本文取5/6)。

彈性模量E、泊松比ν和物質(zhì)密度ρ等物性參數(shù)是坐標(biāo) z的函數(shù)［14］:

其中:n為梯度指數(shù);下標(biāo)c和m分別為陶瓷材料和金屬材料;t為板厚。

2 等幾何有限元法

2.1 節(jié)點(diǎn)矢量和B樣條基函數(shù)

節(jié)點(diǎn)矢量 U={ξ1，ξ2，…，ξn+p+1}由一系列參數(shù)坐標(biāo)點(diǎn)組成，其中 ξi為第 i個節(jié)點(diǎn)，ξi≤ξi+1，i=1，2，…，n+p+1，p為基函數(shù)的階次，n為形成B樣條曲線所需基函數(shù)數(shù)目(與控制節(jié)點(diǎn)數(shù)目相等)。當(dāng)節(jié)點(diǎn)矢量中所有相鄰節(jié)點(diǎn)間距相等(即ξi+1-ξi為常數(shù))，稱之為均勻節(jié)點(diǎn)矢量，反之為不均勻節(jié)點(diǎn)矢量。當(dāng)節(jié)點(diǎn)矢量兩端點(diǎn)(ξ1，ξn+p+1)重復(fù)度為p+1時，稱之為開式節(jié)點(diǎn)矢量。等幾何有限元法通常使用開式節(jié)點(diǎn)矢量，不為零的節(jié)點(diǎn)區(qū)間定義為單元［ξi，ξi+1］。B樣條基函數(shù)在節(jié)點(diǎn)區(qū)間［ξi，ξi+1］內(nèi)具有無窮階可導(dǎo)，在節(jié)點(diǎn) ξi處具有Cp-k連續(xù)，k為節(jié)點(diǎn)重復(fù)出現(xiàn)的次數(shù)，稱為重復(fù)度。給定節(jié)點(diǎn)矢量 U，B 樣條基函數(shù) Ni，p(ξ)定義為［17］:

當(dāng)B樣條基函數(shù)階次p=0和p=1時，等幾何有限元法即為常規(guī)有限元法;當(dāng)p≥2時，等幾何有限元法不同于常規(guī)有限元法。

2.2 NURBS曲面

B樣條曲線定義如下［17］:

其中:Pi=(xiyizi)為第i個控制節(jié)點(diǎn)坐標(biāo)，類似于常規(guī)有限元法的結(jié)點(diǎn)坐標(biāo);Ni，p(ξ)為定義在開式節(jié)點(diǎn)矢量上的p階基函數(shù)。

給定n×m個控制節(jié)點(diǎn)Pij及ξi和ηj方向的節(jié)點(diǎn)矢量 U={ξ1，ξ2，…，ξn+p+1}，V={η1，η2，…，ηm+q+1}。采用張量形式，二維B樣條曲面定義為:

其中:Ni，p(ξ)、Mj，q(η)分別為 ξi和 ηj方向節(jié)點(diǎn)矢量 U和V上定義的B樣條基函數(shù)。

在B樣條曲面中對每個控制節(jié)點(diǎn)Pij乘以一個相應(yīng)的權(quán)因子 ωij得到 NURBS樣條曲面。權(quán)函數(shù)定義為［17］:

NURBS曲面定義為:

2.3 等幾何有限元位移模式

等幾何有限元法采用NURBS基函數(shù)作為結(jié)點(diǎn)形函數(shù)，將控制節(jié)點(diǎn)下標(biāo)ij重新編號為A，其位移逼近可表示為:

其中:uA=［uAvAwAβxAβyA］T為控制節(jié)點(diǎn)A廣義位移矢量;NC表示控制節(jié)點(diǎn)數(shù)目;RA為控制節(jié)點(diǎn)A處形函數(shù)矩陣，是主對角線元素為RA的5階對角方陣。

3 功能梯度板自由振動等幾何有限元控制方程

首先給出基于一階剪切板理論等幾何有限元單元中應(yīng)變場和曲率場的計(jì)算。由式(4)～(6)、(20)得單元中面面內(nèi)應(yīng)變場、曲率場以及橫向剪應(yīng)變場為:

自由振動問題的特征方程為

其中:K和M分別為整體勁度矩陣和質(zhì)量矩陣，ω為自振頻率。

單元對K的貢獻(xiàn)可表示為:

單元對M的貢獻(xiàn)可表示為:

4 數(shù)值算例

4.1 四邊簡支Al/Al2 O3方板

圖2為長厚比L/t=5的四邊簡支Al/Al2O3方板。板的材料參數(shù)為:Aluminum(Al)的彈性模量Em=70 GPa，泊松比 νm=0.3，密度 ρm=2 707 kg·m-3;Alumina(Al2O3)的彈性模量 Ec=380 GPa，泊松比 νc=0.3，密度 ρc=3 800 kg·m-3。梯度指數(shù) n =0、0.5、1、4、10，分別計(jì)算板的自振頻率。采用四種計(jì)算網(wǎng)格(8×8個控制節(jié)點(diǎn)、14×14個控制節(jié)點(diǎn)、20×20個控制節(jié)點(diǎn)、24×24個控制節(jié)點(diǎn))，每種計(jì)算網(wǎng)格分別采用非線性參數(shù)化網(wǎng)格［1］和線性參數(shù)化網(wǎng)格［1］。8×8個控制節(jié)點(diǎn)的非線性參數(shù)化網(wǎng)格和線性參數(shù)化網(wǎng)格如圖3。表1和表2分別給出了不同階次非線性參數(shù)化網(wǎng)格和線性參數(shù)化網(wǎng)格求解的第一階歸一化自振頻率(由表1和表2可知，基于非線性參數(shù)化網(wǎng)格和線性參數(shù)化網(wǎng)格都能獲得收斂的解答，但線性參數(shù)化網(wǎng)格求解精度高和收斂速度快，與文獻(xiàn)［18］常規(guī)材料所得結(jié)果一致。本文后續(xù)算例均采用線性參數(shù)化網(wǎng)格，控制節(jié)點(diǎn)數(shù)為20×20，基函數(shù)階次為p=q=3。

圖2 四邊簡支的方形板Fig.2 Simply supported square plate

針對三種不同的板長厚比，分析了梯度指數(shù)對第一階 自振頻 率的影 響，并 與采 用 ES-DSG3［16］、MITC4［16］、HSDT［19］、kp-Ritz［14］和解析法［20］的結(jié)果做比較，見圖4。歸一化頻率β=ω/ω*，ω為第一階自振頻率的數(shù)值解，ω*第一階自振頻率的解析解［20］。由圖4可知，本文方法求解精度優(yōu)于其它數(shù)值方法。L/t=20時，各種數(shù)值解在n=0.5時的結(jié)果與精確解的偏差明顯大于其他梯度指數(shù)時所對應(yīng)的偏差，其原因尚未有合理的解釋。

圖3 8×8個控制節(jié)點(diǎn)的兩種網(wǎng)格(p=q=3)Fig.3 Two meshesof 8 ×8 control points(p=q=3)

表1 非線性參數(shù)化網(wǎng)格第一階歸一化頻率Tab.1 The normalized first frequency with nonlinear parameterization mesh

表2 線性參數(shù)化網(wǎng)格第一階歸一化頻率Tab.2 The normalized first frequency with linear parameterization mesh

4.2 四邊簡支Ti－6Al－4V/Aluminum oxide方板

四邊簡支的Ti-6Al-4V/Aluminum oxide方板的幾何尺寸為 L=0.4 m，t=0.005 m;材料參數(shù)為:Ti-6Al-4V的彈性模量 Em=105.7 GPa，泊松比 νm=0.298 1，密度 ρm=4 429 kg·m-3;Aluminum oxide 的彈性模量 Ec=320.2 GPa，泊松比 νc=0.26，密度 ρc=3 750 kg·m-3。表3和表4分別給出了梯度指數(shù)為n=0和n=2 000時采用不同方法獲得的Ti-6Al-4V/Aluminum oxide功能梯度板前10階自然頻率。由表3和表4可知本文計(jì)算結(jié)果與 Bishop 等［21］、Zhao等［14］以及He等［10］計(jì)算結(jié)果相符。需要說明的是這兩個表中有些頻率值相同，但其振型并不同，由于篇幅限制，此處沒給出振型圖。

圖4 歸一化第一階頻率與梯度指數(shù)的關(guān)系Fig.4 The normalized first frequency versus gradient index

表3 n=0前10階自然頻率(Hz)Tab.3 The former 10 natural frequencies with n=0(Hz)

表4 n=2 000前10階自然頻率(Hz)Tab.4 The former 10 natural frequencies with n=2 000(Hz)

4.3 梯度指數(shù)、邊界條件及長厚比對自振頻率的影響

考慮 A l/Al2O3，Al/ZrO2，Ti-6Al-4V/Aluminum oxide和SUS304/Si3N4四種材料的方板，分析梯度指數(shù)、邊界條件及長厚比對自振頻率的影響。Al、Al2O3、Ti-6Al-4V和Aluminum oxide的材料參數(shù)同前兩個算例;Zirconia(ZrO2)的彈性模量Ec=151 GPa，泊松比νc=0.3，密度 ρc=3 000 kg·m-3;SUS304 的彈性模量Em=207.78 GPa，泊松比 νm=0.317 7，密度 ρm=8 166 kg·m-3;Si3N4的彈性模量 Ec=322.27 GPa，泊松比 νc=0.24，密度 ρc=2 370 kg·m-3。歸一化頻率為 ω*=

表5 不同梯度指數(shù)時Al/Al2O3方板的歸一化頻率Tab.5 The normalized frequencies versus gradient index for square Al/Al2O3 plate

表5給出了長厚比L/t=10的Al/Al2O3板在四邊簡支(SSSS)，底邊固支三邊自由(CFFF)以及四邊固支(CCCC)條件下前四階歸一化頻率，該結(jié)果與zhao［14］計(jì)算結(jié)果相符。由表5可知，頻率隨著梯度指數(shù)的增加而減小。這與事實(shí)相符，因?yàn)樘荻戎笖?shù)增加，板中陶瓷成分減少，板的整體剛度減小。對于Al/ZrO2，Ti-6Al-4V/Aluminum oxide和SUS304/Si3N4板可以得到類似的結(jié)果，如圖5。從圖5可知，隨著梯度指數(shù)增加，自振頻率減小;Al/Al2O3和Ti-6Al-4V/Aluminum oxide自振頻率相接近;Al/ZrO2自振頻率最大;SUS304/Si3N4自振頻率最小。當(dāng)梯度指數(shù)在0-2之間變化時，自振頻率降低最快;梯度指數(shù)大于5以后自振頻率變化不大。另外，由表5和圖5可知，板的約束越大，其頻率越大。

圖6給出了梯度指數(shù)和長厚比對四邊簡支Al/ZrO2板和四邊固支SUS304/Si3N4板第一階歸一化頻率的影響。由圖6可知，相同梯度指數(shù)，自振頻率隨長厚比增加而增大，當(dāng)長厚比L/t=20后自振頻率變化不大。

圖5 第一階歸一化頻率與梯度指數(shù)的關(guān)系Fig.5 The normalized first frequency versus gradient index

圖6 第一階歸一化頻率與長厚比的關(guān)系Fig.6 The normalized first frequency versus length-to-thickness ratio

5 結(jié)論

基于一階剪切變形理論，建立了分析了功能梯度材料板自由振動問題的NURBS等幾何有限元格式，分析了四種典型功能梯度材料板自振頻率。算例分析表明等幾何有限元法求解精度高和收斂速度快的優(yōu)點(diǎn)。本文推廣了等幾何有限元應(yīng)用范圍。

本文僅討論了基于一階剪切變形理論的等幾何有限元法。下一步將等幾何有限元用于功能梯度板的幾何非線性分析;對于薄板，有效消除剪切自鎖的方法也有待進(jìn)一步研究。

［1］ Hughes T J R，Cottrell J A，Bazilevs Y.Isogeometric analysis:CAD，finite elements，NURBS，exact geometry and mesh refinement［J］. Computer Methods in Applied Mechanics and Engineering，2005，194(39-41):4135-4195.

［2］ Hsu M C，Akkerman I，Bazilevs Y.High-performance computing of wind turbine aerodynamics using isogeometric analysis［J］.Computers＆ Fluids，2011，49(1):93-100.

［3］ Benson D J，Bazilevs Y，Hsu M C，et al.Isogeometric shell analysis:The Reissner-Mindlin shell［J］.Computer Methods in Applied Mechanics and Engineering，2009，199(5-8):276-289.

［4］ Buffa A，Sangalli G，Vázquez R.Isogeometric analysis in electromagnetics:B-splines approximation ［J］.Computer Methods in Applied Mechanics and Engineering，2010，199(17):1143-1152.

［5］Gómez H，Calo V M，Bazilevs Y，et al.Isogeometric analysis of the Cahn-Hilliard phase-field model［J］.Computer Methods in Applied Mechanics and Engineering，2008，197(49-50):4333-4352.

［6］ Seo Y D，Kim H J，Youn S K.Isogeometric topology optimization using trimmed spline surfaces［J］.Computer Methods in Applied Mechanics and Engineering，2010，199(49-52):3270-3296.

［7］ VerhooselC V，Scott M A，Renéde Borst，et al.An isogeometric approach to cohesive zone modeling ［J］.International Journal for Numerical Methods in Engineering，2011，87(1-5):336-360.

［8］ Ghorashi S S， Valizadeh N， Mohammadi S. Extended isogeometric analysis for simulation of stationary and propagating cracks［J］.International Journal for Numerical Methods in Engineering，2011，89(9):1069-1101.

［9］ Koizumi M.FGM activities inJapan［J］.Composites，1997，28(1–2):1-4.

［10］ He X Q，Ng T Y，Sivashanker S，et al.Active control of FGM plates with integrated piezoelectric sensors and actuators［J］.International Journal of Solids and Structures，2001，38(9):1641-1655.

［11］ Sundararajan N，Prakash T，Ganapathi M.Nonlinear free fl exural vibrations of functionally graded rectangular and skew plates under thermal environments［J］.Finite Elements in Analysis and Design，2005，42(2):152-168.

［12］ Reddy J N.Analysis of functionally graded plates［J］.International Journal for Numerical Methods in Engineering，2000，47(1-3):663-684.

［13］ Lee Y Y，Zhao X，Liew K M.Thermoelastic analysis of functionally graded plates using the element-free kp-Ritz method［J］.Smart Materials and Structures，2009，18(3):035007.1-035007.15.

［14］ Zhao X，Lee Y Y，Liew K M.Free vibration analysis of functionally graded plates using the element-free kp-Ritz method［J］.Journal of Sound and Vibration，2009，319(3-5):918-939.

［15］ Nguyen-Xuan H，Tran L V，Thai H，et al.Analysis of functionally graded plates by an eきcient finite element method with node-based strain smoothing［J］. Thin-Walled Structures，2012，54:1-18.

［16］ Nguyen-Xuan H，Tran L V，Nguyen-Thoi T，et al.Analysis of functionally graded plates using an edge-based smoothed finite element method ［J］.Composite Structures，2011，93(11):3019-3039.

［17］Piegl L A，Tiller W.The NURBSBook［M］.Second ed.，New York，Springer-Verlag，1997.

［18］ Shojaee S，Valizadeh N，Izadpanah E，et al.Free vibration and buckling analysis of laminated composite plates using the NURBS-based isogeometric finite element method ［J］.Composite Structures，2012，94(5):1677-1693.

［19］ Matsunaga H.Free vibration and stability of functionally graded plates according to a 2D higher-order deformation theory［J］.Composite Structures，2008，82(4):499-512.

［20］ Hosseini-Hashemi Sh，F(xiàn)adaee M，Atashipour S R.A new exact analytical approach for free vibration of Reissner-Mindlin functionally graded rectangular plates ［J］.International Journal of Mechanical Sciences，2011，53(1):11-22.

［21］ Bishop R E D，Johnson D C.The mechanics of vibration［M］.Cambridge，United Kingdom，Cambridge University press，2011.