基于等幾何有限元法的功能梯度板自由振動分析

尹碩輝，余天堂，劉 鵬

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

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

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

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

圖1所示為金屬-陶瓷功能梯度非均勻材料板，其材料性質沿厚度方向呈現梯度變化。設從板的上表面到下表面，金屬的體積分數從100%過渡到0%，陶瓷的體積分數則由0%到100%。

基于一階剪切板理論［11］，板內任意一點的位移為:

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

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

平面內的應變為:

其中:εm和εb分別為中面面內應變場和板的曲率場，

此外，板的橫向剪應變場為:

根據胡克定律，在正軸下板內的應力-應變關系可表示為:

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

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

彈性模量E、泊松比ν和物質密度ρ等物性參數是坐標 z的函數［14］:

其中:n為梯度指數;下標c和m分別為陶瓷材料和金屬材料;t為板厚。

2 等幾何有限元法

2.1 節點矢量和B樣條基函數

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

當B樣條基函數階次p=0和p=1時，等幾何有限元法即為常規有限元法;當p≥2時，等幾何有限元法不同于常規有限元法。

2.2 NURBS曲面

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

其中:Pi=(xiyizi)為第i個控制節點坐標，類似于常規有限元法的結點坐標;Ni，p(ξ)為定義在開式節點矢量上的p階基函數。

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

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

在B樣條曲面中對每個控制節點Pij乘以一個相應的權因子 ωij得到 NURBS樣條曲面。權函數定義為［17］:

NURBS曲面定義為:

2.3 等幾何有限元位移模式

等幾何有限元法采用NURBS基函數作為結點形函數，將控制節點下標ij重新編號為A，其位移逼近可表示為:

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

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

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

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

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

單元對K的貢獻可表示為:

單元對M的貢獻可表示為:

4 數值算例

4.1 四邊簡支Al/Al2 O3方板

圖2為長厚比L/t=5的四邊簡支Al/Al2O3方板。板的材料參數為: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。梯度指數 n =0、0.5、1、4、10，分別計算板的自振頻率。采用四種計算網格(8×8個控制節點、14×14個控制節點、20×20個控制節點、24×24個控制節點)，每種計算網格分別采用非線性參數化網格［1］和線性參數化網格［1］。8×8個控制節點的非線性參數化網格和線性參數化網格如圖3。表1和表2分別給出了不同階次非線性參數化網格和線性參數化網格求解的第一階歸一化自振頻率(由表1和表2可知，基于非線性參數化網格和線性參數化網格都能獲得收斂的解答，但線性參數化網格求解精度高和收斂速度快，與文獻［18］常規材料所得結果一致。本文后續算例均采用線性參數化網格，控制節點數為20×20，基函數階次為p=q=3。

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

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

圖3 8×8個控制節點的兩種網格(p=q=3)Fig.3 Two meshesof 8 ×8 control points(p=q=3)

表1 非線性參數化網格第一階歸一化頻率Tab.1 The normalized first frequency with nonlinear parameterization mesh

表2 線性參數化網格第一階歸一化頻率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;材料參數為: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分別給出了梯度指數為n=0和n=2 000時采用不同方法獲得的Ti-6Al-4V/Aluminum oxide功能梯度板前10階自然頻率。由表3和表4可知本文計算結果與 Bishop 等［21］、Zhao等［14］以及He等［10］計算結果相符。需要說明的是這兩個表中有些頻率值相同，但其振型并不同，由于篇幅限制，此處沒給出振型圖。

圖4 歸一化第一階頻率與梯度指數的關系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 梯度指數、邊界條件及長厚比對自振頻率的影響

考慮 A l/Al2O3，Al/ZrO2，Ti-6Al-4V/Aluminum oxide和SUS304/Si3N4四種材料的方板，分析梯度指數、邊界條件及長厚比對自振頻率的影響。Al、Al2O3、Ti-6Al-4V和Aluminum oxide的材料參數同前兩個算例;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 不同梯度指數時Al/Al2O3方板的歸一化頻率Tab.5 The normalized frequencies versus gradient index for square Al/Al2O3 plate

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

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

圖5 第一階歸一化頻率與梯度指數的關系Fig.5 The normalized first frequency versus gradient index

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

5 結論

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

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

［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，Fadaee 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.