岳承宇 王立峰 趙永輝

(南京航空航天大學(xué)機械結(jié)構(gòu)力學(xué)及控制國家重點實驗室，南京 210016)

引言

顫振是嚴(yán)重的動氣動彈性不穩(wěn)定現(xiàn)象，當(dāng)飛行器飛行速度超過顫振臨界速度時，結(jié)構(gòu)和空氣動力相互作用，產(chǎn)生自激的振動，往往會引起災(zāi)難性的后果［1］.顫振主動抑制(AFS)技術(shù)是飛行器顫振研究的熱門領(lǐng)域，相比于顫振被動抑制，其優(yōu)點是適應(yīng)性強、效果好、較少增加結(jié)構(gòu)重量.基于主動氣動彈性機翼(AAW)的設(shè)計理念［2］，通常做法是在機翼上布置多個控制面，通過主動控制翼面的變形，改變作用在機翼上的氣動力，從而使系統(tǒng)趨于穩(wěn)定［3-5］.這種方法的不足之處在于能量轉(zhuǎn)換環(huán)節(jié)繁雜，附加部件多，可靠性不足［6］.

智能結(jié)構(gòu)用于顫振主動抑制是一種具有廣闊應(yīng)用前景的設(shè)計理念，含有功能材料的結(jié)構(gòu)除了具有承載能力，還有自我監(jiān)測與作動功能［7］.眾多功能材料中，壓電材料能提供電能與機械能的直接轉(zhuǎn)換，而且分布靈活，結(jié)構(gòu)簡單，響應(yīng)迅速，從而備受關(guān)注.Lazarus等分別使用鋁和石墨/環(huán)氧樹脂為基板材料，將壓電作動器分布于平直板的上下表面，對升力面的主動氣動彈性控制問題進行了理論和實驗研究，風(fēng)洞試驗結(jié)果表明顫振速度提高了11%［8］.在此基礎(chǔ)上，Han等進一步對壓電作動器進行了位置優(yōu)化并應(yīng)用魯棒控制理論以提高控制效果［9］.國內(nèi)由管德指導(dǎo)的小組進行了類似的工作，并且完成了風(fēng)洞試驗，成功利用壓電作動器實施了升力面的顫振主動抑制［10］.

大部分利用壓電材料進行顫振主動抑制的研究都是采用分布式壓電作動器和傳統(tǒng)的傳感器組成的反饋控制系統(tǒng)，而較少使用壓電傳感器.實際上，壓電傳感器基于正壓電效應(yīng)，具有靈敏度高，頻響范圍大的特點.另一個特點是其在空間上呈連續(xù)分布，可以避免傳統(tǒng)的傳感器由于放置位置不當(dāng)引起的個別模態(tài)觀測不到的問題［11］.

已有研究通常使用單輸入或較少輸入的控制系統(tǒng)，文獻［9］研究表明驅(qū)動力單一使得結(jié)構(gòu)缺少扭轉(zhuǎn)力的作用，從而顫振抑制效果受限.

本文使用多個壓電作動器與傳感器分布于整個機翼表面，作動器與傳感器對稱粘貼于上下兩側(cè)，并假設(shè)粘貼緊密.在已有的壓電材料中，壓電陶瓷(PZT)壓電應(yīng)變常數(shù)大，驅(qū)動力強，適合作為作動器，聚偏氟乙烯(PVDF)具有質(zhì)輕、柔韌、靈敏的特點，適合作為傳感器.故本文采用這兩種不同的材料分別作為作動器與傳感器.

復(fù)合材料層合板的有限元理論中，經(jīng)典層合板理論(CLPT)適合薄板，一階剪切變形板理論(FSDT)和高階剪切變形板理論(HSDT)適合薄板和中厚板［12］，在考慮了模型的準(zhǔn)確性和節(jié)省計算規(guī)模后，選擇使用一階剪切變形板理論來建立結(jié)構(gòu)的機電耦合模型，并取其前4階固有模態(tài)作為顫振的振型.本文對這種帶有多個壓電作動器與傳感器的板式智能結(jié)構(gòu)機翼進行了顫振主動抑制的研究.

1 有限元方程

1.1 位移場和電勢場假設(shè)

圖1為帶有壓電層的復(fù)合材料層合板示意圖，壓電鋪層沿z軸極化，板一共有N層，上下兩層分別為壓電作動器層與傳感器層，基板采用正交各向異性復(fù)合材料.根據(jù)一階剪切變形板理論的假設(shè)，板上任意一點的位移為

其中:u0、v0為板中面內(nèi)的位移，w0為橫向位移，θx、θy分別表示中面法線繞y、-x軸的轉(zhuǎn)角.

圖1 壓電層合板示意圖Fig.1 Laminated piezoelectric plate configuration

將式(1)寫為矩陣的形式為

其中:Δ為位移場向量，d為廣義位移向量;Δ、d和矩陣Z的具體表達(dá)見附錄A.

根據(jù)線彈性體的位移-應(yīng)變關(guān)系，由位移場(1)導(dǎo)出的應(yīng)變場表達(dá)式為

其中:εre為去掉了z方向正應(yīng)變的各應(yīng)變分量組成的向量為中面的面內(nèi)正應(yīng)變和切應(yīng)變組成的向量，κ為中面的曲率和扭率組成的向量，γ0為橫向切應(yīng)變組成的向量.式中矩陣的具體表達(dá)式見附錄B.

壓電層與基板接觸的那一面接地，在壓電片的厚度很薄的情況下，可以假設(shè)電勢φk(x，y，z，t)沿厚度的分布為線性的［12］(k=1，N)，寫出壓電作動器層和傳感器層的電勢場分布函數(shù)為

1.2 本構(gòu)方程

壓電材料的正壓電效應(yīng)和逆壓電效應(yīng)可以由如下的本構(gòu)方程表達(dá)［13］

其中:Ck、ek、?k分別為材料的彈性常數(shù)矩陣、壓電應(yīng)力常數(shù)矩陣、介電常數(shù)矩陣，Ek、Dk為電場強度矢量和電位移矢量.

一階剪切變形板理論假設(shè)板沿z方向無正應(yīng)力，故壓電材料的本構(gòu)方程從一般完整的本構(gòu)方程縮減為如下關(guān)系

其中:各矩陣的具體表達(dá)式見附錄C.

基板部分為正交各向異性材料，其本構(gòu)方程為

1.3 有限元插值

整個結(jié)構(gòu)采用四節(jié)點四邊形等參單元進行離散化，單元第i個節(jié)點的廣義位移向量參照式(2)寫為di(i=1，…，4)，而整個單元的節(jié)點廣義位移向量寫為

采用雙線性插值函數(shù)對節(jié)點位移進行插值，廣義位移向量可以由節(jié)點位移表示為

利用(9)式，(3)式中的廣義應(yīng)變可以寫為

其中:應(yīng)變矩陣Bε、Bκ、Bγ的具體表達(dá)式見附錄D.

在解決剪切鎖死的問題上，參考文獻［15-16］的方法，另行假設(shè)橫向剪切應(yīng)變場替代原有的應(yīng)變場γ0，使用單元采樣點處的橫向切應(yīng)變值對整個單元的橫向切應(yīng)變進行獨立插值，最后得到的替代橫向切應(yīng)變表示為

將(10)、(11)式一并寫為

電場強度為電勢的負(fù)梯度，即

將(4)式、(13)式代入(14)式可得電場強度的離散表達(dá)為

1.4 結(jié)構(gòu)動力學(xué)方程

單元的總勢能由彈性體應(yīng)變能和壓電層的電勢能組成，表示為

其中:Ωk為第k層的體積.

單元的動能為

其中:ρk為第k層的質(zhì)量密度.

單元的外力和外加電荷面密度做功為

其中:S為單元的表面積，fs為單元上表面的面力向量，為施加在作動器上的外加電荷面密度.

結(jié)構(gòu)的機電耦合動力學(xué)方程可以由Hamilton原理推出，Hamilton原理表示為

將(3)、(6)、(7)、(12)、(15)式代入(16)式，(2)、(9)式代入(17)式，(2)、(9)、(13)式代入(18)式，再將所得結(jié)果代入Hamilton方程(19)，得到單元的動力學(xué)方程為

其中:矩陣Me為單元質(zhì)量矩陣，矩陣為單元結(jié)構(gòu)剛度矩陣，矩陣、表示壓電作動器和傳感器的外加電載荷和單元節(jié)點電勢差之間的對應(yīng)關(guān)系，具有剛度屬性，矩陣、表示作動器層、傳感器層與基板之間的靜態(tài)機電耦合效應(yīng).這些單元特性矩陣由結(jié)構(gòu)的尺寸、壓電層的分布以及材料的物理特性決定;計算時先沿z方向積分，再采用2×2的Gauss積分在自然坐標(biāo)系下對單元表面進行積分［14］.為單元的等效外力向量表示施加在作動器層的電荷面密度產(chǎn)生的作用力.各矩陣的具體表達(dá)式見附錄F.

實際操作中，作動器層的電勢分布已知，通過控制電勢分布，作動器對結(jié)構(gòu)的作用可以寫為等效外力的形式，式(20)可以寫為

傳感器層的單元節(jié)點電勢差可由式(22)得到

對于沒有壓電層的單元，其運動方程可以由(23)式簡化得到.經(jīng)過組集得到結(jié)構(gòu)的整體運動方程和傳感器電壓輸出方程為

其中:X為結(jié)構(gòu)整體位移自由度向量，u為作動電壓輸入向量，y為傳感電壓輸出向量.

2 顫振特性分析

圖2為機翼模型示意圖，圖中的網(wǎng)格為結(jié)構(gòu)有限元網(wǎng)格劃分，機翼為懸臂式的復(fù)合材料層合板結(jié)構(gòu)，采用T300/976石墨/環(huán)氧樹脂材料，厚度為3 mm，共4層，單層厚度為0.75 mm，鋪層角度為［90/0/0/90］.機翼上分布15個壓電區(qū)，片狀的壓電作動器和傳感器對稱粘貼于機翼上下表面，作動器材料采用PZT-5A，厚度為0.4 mm，傳感器材料采用PVDF，厚度為0.1 mm.每個壓電片的作動電壓為獨立變量，且在單個壓電片上相同，傳感器輸出每個壓電片中部的電壓，計算時輸出電壓由單元節(jié)點電勢差插值得到.

各材料的特性參數(shù)見附錄G.

圖2 機翼模型示意圖Fig.2 Wing model configuration

對結(jié)構(gòu)進行動力特性分析，即求解廣義特征值問題

其中:Φ為固有振型矩陣，Λ為廣義特征值組成的對角陣.

利用(25)式中的K、M矩陣對結(jié)構(gòu)進行動力特性分析時，假設(shè)了電壓輸入u為0，即壓電作動器上分布的電勢差為0，這實際要求作動器閉路.計算機翼的前4階固有模態(tài)如圖3所示.

圖3 機翼固有頻率與固有振型Fig.3 Wing natural frequencies and mode shapes

采用偶極子網(wǎng)格法(DLM)計算頻域亞音速非定常空氣動力［17］，空氣動力網(wǎng)格如圖4所示，空氣動力網(wǎng)格點與結(jié)構(gòu)有限元網(wǎng)格點的連接采用無限板樣條(IPS)方法［18］.(25)式中的F1寫為與結(jié)構(gòu)有關(guān)的頻域氣動力形式

其中:qd為動壓，Qgg為廣義空氣動力矩陣，是減縮頻率和馬赫數(shù)M∞的函數(shù).

圖4 升力面上的空氣動力網(wǎng)格Fig.4 Aerodynamic grids on the lifting surface

引入模態(tài)坐標(biāo)變換

其中:Φ只取結(jié)構(gòu)的前4階振型，q為模態(tài)坐標(biāo)向量.

得到模態(tài)坐標(biāo)下的氣動彈性方程和傳感器電壓輸出方程為

采用V-g法對未控系統(tǒng)進行顫振分析［18］，取海平面處大氣密度ρa=1.225kg/m3，并假設(shè)空氣為不可壓縮流(M∞=0)，得到顫振V-g和V-f圖如圖5所示.從圖中可以看出機翼的顫振速度為51.74 m/s，顫振頻率為19.90 Hz.

3 時域狀態(tài)空間方程

在設(shè)計顫振主動抑制的控制律時，需要用到時域的空氣動力，(30)式中由偶極子網(wǎng)格法計算出的廣義空氣動力矩陣Q，M∞)是在機翼做簡諧運動的條件下得到的，屬于頻域氣動力.為了得到時域氣動力，通常的做法是利用若干離散頻率點處的廣義空氣動力矩陣，將空氣動力矩陣Q在Laplace域用有理函數(shù)擬合表達(dá)，再通過引入空氣動力狀態(tài)變量，將整個方程轉(zhuǎn)換到時域，最后得到被控系統(tǒng)時域下的狀態(tài)空間方程.

本文采用Roger近似法擬合空氣動力矩陣［18］，矩陣Q，M∞)在Laplace域中近似表示為

其中:A0，…，為待定系數(shù)矩陣，γr-2為空氣動力滯后項=sbR/V，s為拉氏變量，bR為參考半弦長，V為來流速度.各待定系數(shù)矩陣根據(jù)已知頻率點處的空氣動力矩陣Q，M∞)求得.

引入空氣動力狀態(tài)變量

由(33)式可得空氣動力狀態(tài)變量自身滿足的時域微分方程為

圖5 顫振V-g和V-f圖Fig.5 FlutterV-g and V-f plots

利用式(32)將(30)式右邊第一項在Laplace域近似表示，引入空氣動力狀態(tài)變量后再作Laplace逆變換，最終(30)式可以寫為

其中:Acon、Bcon、Ccon分別為系統(tǒng)矩陣、輸入矩陣與輸出矩陣，其具體表達(dá)式見附錄H.

4 控制律設(shè)計與仿真分析

最優(yōu)控制理論中［19］，線性定常系統(tǒng)的輸出調(diào)節(jié)器問題表述為尋找最優(yōu)控制律

使得如下二次型性能指標(biāo)達(dá)到極小

其中:Q1、R1為對應(yīng)輸入和輸出的權(quán)系數(shù)矩陣，為了使問題簡化，均取為正定對角陣.權(quán)矩陣的取值由經(jīng)驗給定，合適的權(quán)矩陣使系統(tǒng)盡快鎮(zhèn)定的同時，輸入電壓應(yīng)保持在壓電材料可以承受的范圍.

增益反饋矩陣Fcon由下式求得

其中:P1為如下Riccati方程的解

針對本文的氣動彈性系統(tǒng)，狀態(tài)變量x中包含結(jié)構(gòu)的模態(tài)坐標(biāo)及其對時間的一階導(dǎo)數(shù)，另外還有引入的空氣動力狀態(tài)變量，狀態(tài)變量x實際上無法直接測得.LQG控制增加了最優(yōu)濾波器設(shè)計過程，通過系統(tǒng)的輸入u和輸出y，可得出狀態(tài)變量x的最優(yōu)估計值，表示為

根據(jù)Kalman濾波理論，若系統(tǒng)存在噪聲，對狀態(tài)變量和輸出變量測量過程中Gauss白噪聲的協(xié)方差矩陣分別為Q2和R2，濾波器增益矩陣H為

其中:P2為如下Riccati方程的解

LQG控制律表示為

包括了控制器的閉環(huán)系統(tǒng)結(jié)構(gòu)由圖6所示，系統(tǒng)的狀態(tài)方程為

圖6 LQG控制器結(jié)構(gòu)Fig.6 Block diagram of LQG controller

為了精確擬合氣動力，本文取4個空氣動力滯后項，分別為γ1=0.34，γ2=0.89，γ3=2.56，γ4=13.46，此時系統(tǒng)矩陣Acon為24階.對開環(huán)、閉環(huán)系統(tǒng)進行穩(wěn)定性分析，隨著來流速度的增加，系統(tǒng)特征根軌跡如圖7所示，圖中刪去了零頻率附近無意義模態(tài)的特征根.

圖7 隨著來流速度增加開閉環(huán)系統(tǒng)的根軌跡Fig.7 Root locus of open-loop and close-loop system with increasing incoming flow velocity

未控系統(tǒng)在來流速度達(dá)到51.74m/s時開始出現(xiàn)正實部的特征根，即發(fā)生了顫振，其結(jié)果與上文V-g法顫振分析得到的結(jié)果相同.閉環(huán)系統(tǒng)的特征根全部在復(fù)平面的左半平面，表明受控系統(tǒng)是穩(wěn)定的.在流速較大時(V=100m/s)，系統(tǒng)仍保持穩(wěn)定，這個結(jié)果是在未考慮壓電作動器的最大承載電壓的情況下得到的，實際上在流速較大時，輕微的擾動也需要作動器提供較大的驅(qū)動力才能使系統(tǒng)鎮(zhèn)定，作動電壓有可能超出了壓電材料所能承受的極限，此種情況本文未作考慮.另一方面，本文使用了15個壓電作動器分布在整個機翼表面，對機翼前數(shù)階模態(tài)的聯(lián)合顫振有較全面的控制能力，在更多模態(tài)參與下的顫振抑制效果還有待研究.

對系統(tǒng)進行動響應(yīng)仿真計算，采用Runge-Kutta法求解系統(tǒng)的動響應(yīng)，設(shè)定來流速度V=54m/s，在t=0.5s時施加控制，結(jié)果如圖8所示.

圖中位移響應(yīng)為翼梢后緣點的橫向位移時間歷程，可以看到系統(tǒng)在初始擾動下發(fā)散運動，施加控制后大約在t=1s時系統(tǒng)趨于鎮(zhèn)定.圖中給出的控制電壓是翼根處的三塊壓電作動器的控制電壓，一般情況下所有作動器的控制電壓最大值不超過第一塊作動器的最大控制電壓，LQG控制律的算法是讓總的控制能量盡可能的小，第一塊壓電作動器在控制過程中起到了較大作用.

圖8 系統(tǒng)的時間歷程(V=54m/s)Fig.8 Time histories of the system at V=54m/s

5 結(jié)論

本文研究了一個帶有分布式壓電作動器與傳感器的復(fù)合材料板結(jié)構(gòu)機翼的升力面顫振的主動抑制問題.多對壓電作動器與傳感器對稱粘貼在機翼上下表面，作動器和傳感器采用不同的壓電材料，傳感器輸出由于結(jié)構(gòu)變形產(chǎn)生的電壓，作動器主動輸入控制電壓使作動器對結(jié)構(gòu)產(chǎn)生作用力.通過主動控制律的實施，升力面顫振被較好地抑制，證明了本文顫振主動抑制方法的有效性.

本文較詳細(xì)地推導(dǎo)了復(fù)合材料壓電層合板的有限元方程，建模方法不局限于本文中上下兩面粘貼有壓電片的板結(jié)構(gòu)模型，還可以推廣到壓電片嵌于復(fù)合材料板之中的結(jié)構(gòu)，壓電層數(shù)也可以是多層.

今后可研究作動器和傳感器的位置優(yōu)化以提高控制效率，研究考慮模型不確定性的魯棒控制或自適應(yīng)控制以接近工程實際.

1 Librescu L，Marzocca P.Advances in the linear/nonlinear control of aeroelastic structural systems.Acta Mechanica，2005，178(3-4):147～186

2 吳志剛，楊超.主動氣動彈性機翼的顫振主動抑制與陣風(fēng)減緩研究.機械強度，2003，25(1):32～35(Wu Z G，Yang C.Investigation on active flutter suppression and gust alleviation for an active aeroelastic wing.Journal of Mechanical Strength，2003，25(1):32～35(in Chinese))

3 Borglund D，Kuttenkeuler J.Active wing flutter suppression using a trailing edge flap.Journal of Fluids and Structures，2002，16(3):271～294

4 Zhao Y H.Flutter suppression of a high aspect-ratio wing with multiple control surfaces.Journal of Sound and Vibration，2009，324(3):490～513

5 楊仕福，王琪，趙永輝等.大展弦比機翼μ控制.動力學(xué)與控制學(xué)報，2011，9(1):79～85(Yang SF，Wang Q，Zhao Y H，et al.μControl of a high-aspect-ratio wing using multiple actuators.Journal of Dynamics and Control，2011，9(1):79～85(in Chinese))

6 任勇生，劉立厚，韓景龍等.飛行器非線性氣動彈性和顫振主動控制研究進展.力學(xué)季刊，2003，24(4):534～540(Ren Y S，Liu L H，Han JL，et al.Advance in nonlinear aeroelasticity study and active flutter suppression of aircraft.Chinese Quarterly of Mechanics，2003，24(4):534～540(in Chinese))

7 陳勇，熊克，王鑫偉等.飛行器智能結(jié)構(gòu)系統(tǒng)研究進展與關(guān)鍵問題.航空學(xué)報，2004，25(1):21～25(Chen Y，Xiong K，Wang X W，et al.Progress and challenges in aeronautical smart structure systems.Acta Aeronautica et Astronautica Sinica，2004，25(1):21～25(in Chinese))

8 Lazarus K B，Crawley E F，Lin C Y.Multivariable active lifting surface control using strain actuation:analytical and experimental results.Journal of Aircraft，1997，34(3):313～321

9 Han JH，Tani J，Qiu J.Active flutter suppression of a lifting surface using piezoelectric actuation and modern control theory.Journal of Sound and Vibration，2006，291(3):706～722

10 陳偉民，管德，李敏等.采用分布式壓電驅(qū)動器升力面的顫振主動抑制.力學(xué)學(xué)報，2002，34(5):756～763(Chen W M，Guan D，Li M，et al.Flutter suppression using distributed piezoelectric actuators.Acta Mechanica Sinica，2002，34(5):756～763(in Chinese))

11 Tzou H S，Tseng CI.Distributed piezoelectric sensor/actuator design for dynamic measurement/control of distributed parameter systems:a piezoelectric finite element approach.Journal of Sound and Vibration，1990，138(1):17～34

12 Samanta B，Ray M C，Bhattacharyya R.Finite element model for active control of intelligent structures.AIAA Journal，1996，34(9):1885～1893

13 Benjeddou A，DeüJ F，Letombe S.Free vibrations of simply-supported piezoelectric adaptive plates:an exact sandwich formulation.Thin-Walled Structures，2002，40(7-8):573～593

14 王勖成.有限單元法.北京:清華大學(xué)出版社，2003:111～112，145～147(Wang X C.Finite element method.Beijing:Tsinghua University Press，2003:111～112，145～147(in Chinese))

15 Hinton E，Huang H C.A family of quadrilateral Mindlin plate elements with substitute shear strain fields.Computers＆structures，1986，23(3):409～431

16 Nayak A K，Moy SSJ，Shenoi R A.Free vibration analysis of composite sandwich plates based on Reddy＇s higher-order theory.Composites Part B:Engineering，2002，33(7):505～519

17 Albano E，Rodden WP.A doublet-lattice method for calculating lift distributions on oscillating surfaces in subsonic flows.AIAA Journal，1969，7(2):279～285

18 趙永輝.氣動彈性力學(xué)與控制.北京:科學(xué)出版社，2007:246～250，278～280，310～315(Zhao Y H.Aeroelasticity and control.Beijing:Science Press，2007:246～250，278～280，310～315(in Chinese))

19 顧仲權(quán)，馬扣根，陳衛(wèi)東.振動主動控制.北京:國防工業(yè)出版社，1997:68～71(Gu Z Q，Ma K G，Chen W D.Active vibration control.Beijing:National Defense Industry Press，1997:68～71(in Chinese))

20 Ha SK，Keilers C，Chang FK.Finite element analysis of composite structures containing distributed piezoceramic sensors and actuators.AIAAJournal，1992，30(3):772～780

21 Xu K，Noor A K，Tang Y Y.Three-dimensional solutions for free vibrations of initially-stressed thermoelectroelastic multilayered plates.Computer Methods in Applied Mechanics and Engineering，1997，141(1):125～139

22 Shen S，Kuang Z B.An active control model of laminated piezothermoelastic plate.International Journal of Solids and Structures，1999，36(13):1925～1947

23 Mitchell JA，Reddy JN.A refined hybrid plate theory for composite laminates with piezoelectric laminae.International Journal of Solids and Structures，1995，32(16):2345～2367