基于三階剪切變形理論的壓電功能梯度板靜力學等幾何分析

劉 濤， 李朝東， 汪 超， 蔣雅芬， 劉慶運

(1.上海大學 機電工程與自動化學院，上海 200072； 2.安徽工業大學 機械工程學院，馬鞍山 243002； 3.安徽工業大學 創新教育學院，馬鞍山 243002)

功能梯度材料(Functionally Graded Material, FGM)[1]是指沿著某方向連續性改變材料結構或組成的材料，如在厚度方向，材質從金屬過渡到陶瓷材料。因其獨特的力學優勢，目前被廣泛應用于航空、核能、生物醫學、機械等領域[2-4]。同樣作為智能材料，壓電材料獨有的正逆壓電效應[5]，使其在結構的健康監測[6]和微機電系統[7]領域中發揮著重要作用。因此，將壓電材料與FGM結合以獲得性能更加優越的智能結構和智能材料，引起了國內外學者的廣泛關注。

壓電智能結構的迅速發展使其結構形式變得更加多樣化、復雜化。因此，結構的自由振動[8-9]、彎曲[10-12]、屈曲[13]、振動[14]、形狀控制等問題的研究愈發受到重視[15]。壓電智能結構分析的建模方法主要有解析法建模、數值方法建模和實驗方法建模。壓電功能梯度板(Piezoelectric Functionally Graded Plate, PFGP)是最基本的壓電功能梯度材料結構形式之一，就其結構分析的數值分析方法看，采用常規的三維彈性理論框架，傳統的數值方法會耗費大量計算網格。若要提高計算效率，則需要根據材料屬性定制對應的形函數[16]改進數值方法(如梯度有限先法)。然而，這種改進方式仍然擺脫不了維度上帶來的網格自由度增量。為解決上述問題，國內外很多學者采用等效單層理論[17]框架下的數值方法去研究PFGP的力學行為。He等與Liew等分別在經典板理論(Classical Plate Theory, CPT)和一階剪切變形理論(First-Order Shear Deformation Theory, FSDT)的框架下，利用有限單元法(Finite Element Method, FEM)研究了PFGP在機-電-熱載荷下的靜力學、動力學特性和主動振動控制問題。此外，為解決FEM中因單元畸變降低計算精度[18-19]的問題，文獻[12-14]采用無網格法(Mesh-Free Method)研究了PFGP的屈曲、靜態響應、形狀控制以及振動等問題。Loja等[20]建立了高階剪切變形理論(Higher-Order Shear Deformation Theory, HSDT)框架下的B樣條有限元模型，文中有效分析了PFGP的靜態響應與自由振動問題。結合HSDT和馮·卡爾曼方程(Von Kármán Euqations)，Fakhari等[21-22]又采用有限單元法研究了熱環境下PFGP的非線性自由、強迫振動等問題。上述板理論中，CPT僅對薄板有效而不具有通用性；FSDT因階次過低而可能會導致分析精度稍弱；而HSDT構造過程過于復雜而不利于數值分析的離散化過程。相對而言，三階剪切變形理論(Third-Order Shear Deformation Theory, TSDT)因利用C1連續單元提高求解的精確度且無需引入剪切因子具有明顯的優勢。目前，已有一批學者將TSDT的數值方法運用到功能梯度梁、板及殼體的力學分析中[23-25]。在PFGP的研究中，Selim等[26]構建了TSDT框架下的無網格法，文中研究了兩種不同結構的PFGPs的主動振動控制問題。

但就上述的數值方法(有限單元法或無網格法)而言，若研究對象的幾何體稍作復雜，則會出現計算分析模型和幾何模型不一致的現象[27]。為解決此類問題，Hughes提出一種新的力學數值算法——等幾何分析(Isogeometric Analysis, IGA)[28]。其基本思想是將計算機輔助設計(Computer Aided Design, CAD)的非均勻有理B樣條(Non-Uniform Rational B-Splines, NURBS)作為計算機輔助工程(Computer Aided Engineering, CAE)的形函數，方便地實現任意連續階次的平滑性。與傳統的數值方法相比，IGA這種統一的數學方程，即省去了CAD與CAE數據之間的交互，又避免了如有限元使用拉格朗日多項式差值函數帶來的逼近誤差。可以說，理論上IGA形成了CAD與CAE的無縫連接。鑒于IGA方法自身的諸多優點以及其天然地滿足TSDT所需的C1連續性要求[29]的特性，使得部分學者將TSDT與IGA方法結合對FGM結構的靜力學[30]、動力學[31-32]、結構與材料參數優化[33-34]等問題進行了研究，驗證了該方法的精確性。

最近，少部分學者將IGA方法引入到PFGP的分析中。如，Phung-Van等[35]分析了PFGP在機-電-熱負載下的非線性瞬態響應；Nguyen等[36]研究了壓電功能梯度多孔板的靜動態響應；Nguyen-Quang[37]則針對均勻、V型、X型以及O型四種體積分數分布形式，分析了集成壓電層的功能梯度碳納米管增強復合材料板的動態響應。值得注意的是，在PFGP靜力學問題中，構建三階剪切變形理論的等幾何分析方法還未被涉及，因此，本文擬用TSDT下的IGA方法，研究并探討PFGP的自由振動、彎曲以及變形控制等靜力學問題。

1 基于三階剪切變形理論的壓電功能梯度板模型

如圖1所示，PFGP長為a，寬為b，共分為上、中、下三層。中間層厚度為hf，是由金屬和陶瓷組成的功能梯度材料，在厚度方向上其材質具有光滑連續的變化特性。中間層上、下表面均被一厚度為hp的壓電層完全覆蓋，忽略中間層與壓電層之間的膠黏層。板的總厚度ht=hf+2hp。

圖1 壓電功能梯度板

1.1 功能梯度層的體積分數

假設金屬材料的體積分數含量在厚度方向上為冪函數形式分布

(1)

Vc(z)=1-Vm(z)

(2)

則功能梯度板內任意一點的材料的物性參數為

P(z)=(Pm-Pc)Vm(z)+Pc

(3)

式中:下標c和m分別為非金屬和金屬;z為厚度方向上的坐標(z∈[-hf/2,hf/2]);n為功能梯度指數n∈[0,∞];Pc和Pm分別為非金屬和金屬的相關材料屬性，如：密度、彈性模量或者泊松比等。

1.2 三階剪切變形理論

由Reddy[38]提出的TSDT，可得到PFGP內部任意一點的位移為

u(x,y,z)=u0(x,y)+zβx+c1z3(βx+w0,x)

v(x,y,z)=v0(x,y)+zβy+c1z3(βy+w0,y)

w(x,y,z)=w0(x,y)

(4)

由幾何方程可得幾何非零應變為

(5)

其中，面內應變與橫向剪切應變分別為

ε={εxxεyyγxy}T=ε0+zκ1+z3κ2

(6)

γ={γxzγyz}T=εs+z2κs

(7)

式中:ε0={u0,xv0,yu0,y+v0,x}T,

κ1={βx,xβy,yβx,y+βy,x}T，

κ2=c1{βx,x+w0,xxβy,y+w0,yyβx,y+βy,x+2w0,xy}T，

εs={βx+w0,xβy+w0,y}T

κs=3c1{βx+w0,xβy+w0,y}T

由廣義胡克定律，應力的本構關系可以表示為

(8)

彈性常量為

Q66=G12,Q44=G23,Q55=G13

(9)

式中，E，v分別為彈性模量與泊松比。

實際應用中，壓電材料常處于第二類邊界條件，即力學夾持和電學短路[39]。根據第二類壓電方程[40]，PFGP的本構方程可表示為

(10)

(11)

式中，壓電應力常數e31、e32可以通過式(12)獲得。

(12)

式中，d31，d32為壓電應變常數。

僅考慮z方向上的電場強度，對應的電場強度矩陣表示為

(13)

式中，φ為壓電層厚度方向上的電勢差。

根據上述分析，PFGP的應變能為

(14)

因此，彈性剛度矩陣為

(15)

式中，

(16)

(17)

1.3 弱形式

利用哈密頓變分原理[41]，板的弱形式表示為

(18)

式中，t1和t2分別為時間軸的某兩個時刻，L為總能量

(19)

2 壓電功能梯度板的等幾何分析

2.1 NURBS基函數

節點矢量k(ξ)由一組單調不遞減的實數序列組成[42]。定義為k(ξ)={ξ1=0,…,ξi,…,ξn+p+1=1}，其中ξi稱為節點，n為基函數個數，p為樣條階次。Ni,p(ξ)表示第i個p次B樣條基函數

(20)

(21)

對于二維平面問題，NURBS基函數由兩個方向的一維B樣條基函數張量積構成

(22)

式中:wi,j為權重;Ni,p(ξ)和Nj,q(η)分別為ξ方向p階次和η方向q階次的B樣條基函數，Nj,q(η)形式和Ni,p(ξ)類似。

2.2 機械位移場離散化

利用NURBS基函數，PFGP的近似機械位移場u(ξ,η)可表示為

(23)

其中，dI=[u0Iv0IβxIβyIwI]T為未知量，是控制點I的位移矢量。RI(ξ,η)為式(21)所定義的形函數。將式(23)代入式(6)和(7)，得到

(24)

其中，

2.3 電勢場離散化

在近似電勢場時，可分別將上、下壓電層沿著厚度方向離散為多個子層單元。假設電勢在每個子層厚度方向上的變化為線性[43]，則壓電層的近似電勢場可以描述為

(26)

φi=[φi-1φi] (i=1,2,…,nsub)

(27)

式中，nsub為壓電子層的總數量。

這里，假定每個壓電子層單元同一高度的電勢相同，則由式(13)可得每個子層單元的電場為

(28)

2.4 控制方程

將式(10)、(24)、(26)、(28)代入式(18)，可得壓電功能梯度板的運動控制方程為

(29)

式中:Muu為質量矩陣；Kuu為PFGP的剛度矩陣，Kφφ為壓電材料的自適應剛度矩陣；Kuφ與Kφu為機電耦合剛度矩陣；f為機械載荷；Q為壓電層表面的電荷面密度等效節點電量。其值分別為

(30)

(31)

(32)

(33)

N={N1N2N3}T

(34)

(35a)

(35b)

(35c)

(36)

(37)

(38)

(39)

(40)

因此，Kuφ可以重寫為

(41)

將式(29)的第二行，代入第一行消去電勢φ，可以得到

(42a)

M=Muu

(42b)

(42c)

最后，靜力學、自由振動的問題可描述為

Kd=F

(43)

(44)

式中:ωn為第n階固有頻率;K和M分別為整體剛度矩陣、質量矩陣。

2.5 壓電功能梯度板變形控制

PFGP的變形控制主要分為開環控制和閉環控制。開環控制時，上下壓電層都作為驅動器時，當外部輸入壓電作用于壓電層，電壓將轉化為驅動力。通過調整輸入電壓的大小改變結構的位移，達到控制結構變形的目的。而在閉環控制中，一個壓電層作為驅動器，另一個壓電層作為傳感器。通過反饋控制原理，結構在機械載荷作用下產生機械變形，傳感器層感知結構的變形而產生相應的輸出電壓。該感應電壓通過控制器控制增益放大后，作為驅動器的輸入電壓，驅動器由于逆壓電效應產生相應的反驅動力抑制結構的變形，從而達到控制結構變形的目的。

因此，閉環控制中，式(29)的第二行可以改寫為

[Kφu]ad-[Kφφ]aφa=Qa

(45)

[Kφu]sd-[Kφφ]sφs=Qs

(46)

式中,下標a，s分別為驅動器層和傳感器層。

根據位移反饋控制規律，驅動器的輸入電壓φa表示為

φa=Gdφs

(47)

式中，Gd為位移反饋控制增益系數。

靜態變形中，對于傳感器層，外部電荷Qs=0。忽略逆壓電效應，則由式(46)可得，結構變形傳感器層所產生的電勢為

(48)

將式(47)、(48)代入式(45)可得

(49)

最后，將上式代入(43)，整理可得

(50)

(51)

3 數值算例

為方便起見，在本文后續章節中，將板的機械邊界條件簡寫為：簡支(S)，固支(C)，自由(F)。此外，基函數階次均取為p=q=3。

3.1 收斂性分析

為驗證所提方法的收斂性與適用性，本節以一個尺寸為400 mm×400 mm(a×b)的壓電功能梯度板為例。板的中間層是由Ti-6Al-4V和Aluminum oxide組成的功能梯度材料，厚度hf=5 mm。由式(1)可知，當功能梯度指數n=0時，中間層是材料為Ti-6Al-4V的均質板。當n=∞時，為Aluminum oxide。兩塊沿z向極化，厚度為hp= 0.1 mm的PZT-G1195N壓電層完全覆蓋在中間層的上、下表面。該板的材料參數如表1所示。

如圖2所示，本文利用4種不同控制點數目的IGA計算網格驗證所提方法的收斂性與精確性。表2為對應4種計算網格下，邊界條件為SSSS的壓電功能梯度板的第1、2、4、6、8階固有頻率。由表2可知，本文所提方法計算結果與文獻[44]的解析解計算結果吻合程度較高。此外，在控制節點大于19×19時，本文方法對計算結果的影響較小，可以判斷，該網格自由度下，本文方法已經達到了收斂狀態。因此，本文后續算例均采用的控制節點數目為19×19。表3列舉了控制節點數目19×19下，邊界條件為CFFF和CCCC的PFGPs的前6階固有頻率。通過與文獻[8]的有限單元法計算結果相比，驗證了所提方法對于其他機械邊界條件壓電功能梯度板的適用性。同時，由表1、2可知，隨著功能梯度指數n的增大，固有頻率增大。

表1 材料參數

(a) 9×9

(b) 15×15

(c) 19×19

(d) 21×21

表2 不同網格下，四邊簡支(SSSS)PFGP的1、2、4、6、8階固有頻率

表2(續) Hz

表3 懸臂(CFFF)、四邊固支(CCCC)的PFGPs的前6階固有頻率

3.2 自由振動問題分析

本節，以一個中間層由Al和Al2O3組成，上下壓電層的材料為PZT-4的PFGP為例，研究電學邊界條件(短路和開路)對其固有頻率的影響(PZT-4的材料參數見文獻[44])。值得注意的是，在本節中，當功能梯度指數n=0時，中間層是材料為Al2O3的均質板。當n=∞時，為Al。

表4列出了功能梯度層的寬厚比(hf/a)、壓電層與功能梯度層厚度比(hp/hf)、功能梯度指數n以及電學邊界條件對四邊簡支(SSSS)的PFGP的1、2階固有頻率的影響。表5列舉了寬厚比為0.02(hf/a=0.02)、壓電層與功能梯度層厚度比為0.1(hp/hf=0.1)、機械邊界條件為SSSC和SCSC的PFGP的前6階固有頻率。

由表4、5可知：開路狀態下板的固有頻率大于短路狀態；隨著功能梯度指數n增加，板的金屬材料成分增加，剛度減小，固有頻率也減小；另外，當功能梯度層厚度與板長之比(hf/a)由0.05增加到0.1時，板的固有頻率也增加近一倍。圖3則為開路狀態下板的前6階固有振型。

表4 短路和開路條件下，四邊簡支(SSSS) PFGP的前兩階固有頻率

表5 機械邊界條件為SSSC和SCSC的PFGP在短路和開路狀態下前6階固有頻率

(a) 一階振型

(b) 二階振型

(c) 三階振型

(d) 四階振型

(e) 五階振型

(f) 六階振型

3.3 靜態彎曲與變形控制

3.3.1 壓電雙晶懸臂梁的彎曲分析

如圖4所示，尺寸為100 mm×5 mm×1 mm的壓電雙晶懸臂梁，其上、下層壓電材料采用極化方向相反的PVDF材料。材料參數為E1=E2=2.0 GPa，G12=1.0 GPa，v12=0，e31=e32=0.046 C/m2，k11=k22=k33=10.62×10-9F/m。當對上下壓電層施加外電場時，懸臂梁因逆壓電效應而產生彎曲變形。

圖4 壓電雙晶梁(mm)

在單位電壓(1V)作用下，懸臂梁上不同位置的撓度值如表6所示，通過對比可知，本文方法與文獻[45-46]的數值方法以及文獻[47]的解析解結果基本吻合。

表6 單位電壓下的壓電雙晶梁的撓度

如圖5所示，隨著輸入電壓的增大，壓電雙晶梁的中心線撓度也隨之增大。不同電壓作用下梁末端的撓度值如表7所示，對比可知，本文與文獻[45,47]的計算結果基本一致。

圖5 不同電壓作用下壓電雙晶梁中心線撓度

表7 不同電壓下的壓電雙晶梁末端撓度

3.3.2 壓電功能梯度方板的彎曲分析

本節對一個邊界條件為CFFF的PFGP進行分析。板的材料和尺寸均與3.1節算例相同。

圖6給出了板面受均布載荷q0=100 N/m2作用下，功能梯度指數n對板中心線撓度的影響。可以看出，隨著n的增大，懸臂板中心線的撓度值逐漸變小，該結果與文獻[8,48]的結果相吻合。圖7為PFGP的上、下壓電層僅受10 V的反向電壓時的撓度值變化，此時，壓電層均作為驅動器，即板的開環變形控制，可以看出本文方法與文獻[36,48]趨勢吻合。

圖6 均布載荷作用下n對懸臂板中心線撓度的影響

圖7 懸臂PFGP在10 V電壓作用下的中心線撓度

此外，為研究機電耦合情況下的撓度變化，對懸臂板同時施加q0=100 N/m2的均布機械載荷和不同的輸入電壓時板的中心線撓度變化如圖8所示，可以看出，中心線撓度變化的趨勢與文獻[36,48]相吻合。表8統計了對應的末端位移值，可以發現，隨著電壓的增大，板中心線末端撓度值逐漸減小。

(a) n=0

(b) n=0.5

(c) n=5

(d) n=∞

為探討機械邊界條件對PFGP形狀變化的影響，圖9分析了功能梯度指數n=2時，不同輸入電壓以及四種不同邊界條件下指定位置的撓度變化情況。可以看到，隨著電壓的增大，撓度呈線性變化。值得注意的是，對于懸臂PFGP，板末端向上翹起(末端撓度值為正)；但對于邊界條件為SFSF、SSSS、SCSC的PFGP，由于機械邊界條件的改變，板中心向下凹陷(中心點撓度值為負)。產生這種現象的原因如圖1所示，當上壓電層的輸入電壓與極化方向相反時，逆壓電效應使得其沿長度方向收縮。下壓電層厚度方向的輸入電壓與極化方向相同時，其沿長度方向伸長。機械邊界條件的改變使得板的彎曲方向發生了變化。

表8 均布載荷與不同電壓作用下懸臂PFGP中心線末端撓度

圖9 不同電壓作用下邊界條件對PFGP撓度的影響(n=2)

3.3.3 懸臂壓電功能梯度板的閉環變形控制分析

最后，分析了邊界條件為CFFF和SSSS的PFGP的閉環變形控制問題。算例中，板的尺寸和材料與3.1節相同。上壓電層作為驅動器，下壓電層作為傳感器。如圖10、11所示，在機械均布載荷q0=100 N/m2作用下，隨著增益Gd的增大，板中心線撓度逐漸減小。這是因為傳感器層因機械變形(正壓電效應)產生的輸出電壓經位移反饋控制增益Gd放大后作為驅動器層的輸入電壓。該輸入電壓產生相應的反驅動力抑制結構的變形。

4 結 論

本文為壓電功能梯度板的靜力學分析提供了一個新的高精度數值分析方法。該方法根據哈密頓變分原理以及相關壓電方程，構建了三階剪切變形理論下的功能梯度板等幾何有限元形式。為全面探討該數值方法的有效性，本文分析了自由振動的收斂性問題，并研究了功能梯度層的相關梯度指數n、寬厚比以及壓電層與功能梯度層的厚度比對壓電功能梯度板固有頻率的影響。此外，文中對機械載荷、電載荷、機-電載荷下壓電功能梯度板的靜態彎曲響應也進行了全面討論。最后，本文研究了壓電功能梯度板的閉環變形控制問題。

圖10 反饋增益Gd對懸臂PFGP中心線撓度的影響

圖11 反饋增益Gd對四邊簡支PFGP中心線撓度的影響

靜力特性是壓電功能梯度材料結構最基本的力學特性，但該材料的動力學響應亦是工程實際中經常面對的問題，因此，本文方法下的相關動力學研究是一個值得關注的課題。此外，在熱環境下壓電功能梯度材料的力學分析中，本文方法的適用性問題亦是一個值得關注的領域。