全動(dòng)舵系統(tǒng)柔性多體動(dòng)力學(xué)建模方法

楊執(zhí)鈞，張 忠，李 芮，高 博，郭 靜

（1.北京強(qiáng)度環(huán)境研究所 可靠性與環(huán)境工程技術(shù)重點(diǎn)實(shí)驗(yàn)室； 2.北京強(qiáng)度環(huán)境研究所：北京 100076）

0 引言

全動(dòng)舵系統(tǒng)是由舵機(jī)、傳動(dòng)機(jī)構(gòu)與氣動(dòng)力作用下的彈性舵面共同形成的耦合系統(tǒng)，存在著流、固、電、磁等多物理場的耦合作用，可能在一定的條件下引發(fā)顫振失穩(wěn)[1-3]；而一旦發(fā)生顫振失穩(wěn)，將會(huì)對(duì)航天飛行器造成非常嚴(yán)重的不良后果：結(jié)構(gòu)的振動(dòng)響應(yīng)會(huì)被不斷放大，甚至?xí)斐山Y(jié)構(gòu)破壞、飛行器姿態(tài)失控，乃至飛行失利。

目前國內(nèi)外的相關(guān)研究主要集中于氣動(dòng)力對(duì)于舵面顫振特性的影響方面，常采用計(jì)算流體動(dòng)力學(xué)/計(jì)算結(jié)構(gòu)動(dòng)力學(xué)（CFD/CSD）方法開展全動(dòng)舵面氣動(dòng)彈性響應(yīng)計(jì)算。針對(duì)舵軸間隙對(duì)顫振特性的影響[4-7]，學(xué)者們多基于虛擬質(zhì)量法把含間隙的非線性系統(tǒng)劃分為3 個(gè)線性子系統(tǒng)，采用虛擬質(zhì)量模態(tài)作為統(tǒng)一的坐標(biāo)去表示各子系統(tǒng)。黃程德等[8]發(fā)展了基于CFD 和虛擬質(zhì)量法耦合的間隙非線性氣動(dòng)彈性分析方法，發(fā)現(xiàn)全動(dòng)舵面極限環(huán)幅值隨著間隙角的增大而增大，且跨聲速極限環(huán)的臨界速度比線性顫振的跨聲速凹坑下降20%。但基于虛擬質(zhì)量法建模僅能考察舵軸間隙的影響，難以拓展到舵機(jī)等其他耦合系統(tǒng)。針對(duì)舵機(jī)對(duì)顫振特性影響的研究主要采用假設(shè)模態(tài)法，即將舵機(jī)系統(tǒng)簡化為具有復(fù)剛度的彈簧系統(tǒng)[9-10]。Shin 等[11]對(duì)帶有間隙環(huán)節(jié)和電動(dòng)舵機(jī)的舵面進(jìn)行了氣動(dòng)彈性分析，結(jié)果表明：考慮舵機(jī)動(dòng)剛度與不考慮時(shí)計(jì)算得到的顫振速度和顫振頻率均存在差別，而且非線性的影響會(huì)改變線性顫振特性，使得舵面出現(xiàn)極限環(huán)振蕩現(xiàn)象。丁偉濤等[12]基于典型舵面二階局部剛化模態(tài)和彈性模態(tài)建立了考慮舵機(jī)復(fù)剛度的舵系統(tǒng)動(dòng)力學(xué)模型，研究熱環(huán)境、電動(dòng)舵機(jī)設(shè)計(jì)參數(shù)以及指令信號(hào)幅值對(duì)顫振速度的影響，并提出在舵機(jī)電流環(huán)加入超前滯后環(huán)節(jié)的顫振抑制措施。但將舵機(jī)簡化為復(fù)剛度彈簧時(shí)，其機(jī)電耦合特性、傳動(dòng)特性等需進(jìn)行線性化處理，會(huì)降低仿真模型與真實(shí)模型的匹配程度。

為了能夠較為簡單且有效地耦合舵機(jī)模型、舵軸間隙等因素的影響，本文采用剛?cè)狁詈隙囿w動(dòng)力學(xué)方法建立典型舵面動(dòng)力學(xué)模型，并采用地面顫振試驗(yàn)的集中式氣動(dòng)力模擬方法開展氣動(dòng)彈性耦合仿真，分析舵機(jī)模型、舵軸間隙等因素對(duì)動(dòng)力學(xué)響應(yīng)的影響。

1 全動(dòng)舵系統(tǒng)建模方法

1.1 剛?cè)狁詈隙囿w動(dòng)力學(xué)模型

多體動(dòng)力學(xué)的研究對(duì)象是由多個(gè)具有運(yùn)動(dòng)學(xué)約束、存在大范圍相對(duì)運(yùn)動(dòng)的剛性或柔性部件組成的復(fù)雜系統(tǒng)，主要研究內(nèi)容是這類系統(tǒng)的動(dòng)力學(xué)建模、計(jì)算和控制。對(duì)于構(gòu)件在運(yùn)動(dòng)過程中同時(shí)存在彈性變形的問題，基于Craig-Bampton 方法[13]，采用固連在結(jié)構(gòu)體上的體坐標(biāo)系運(yùn)動(dòng)來描述結(jié)構(gòu)體空間大范圍運(yùn)動(dòng)，相對(duì)于體坐標(biāo)系來描述柔性變形引起的運(yùn)動(dòng)，從而建立能夠考慮彈性變形與非線性大范圍運(yùn)動(dòng)相互耦合影響的柔性多體動(dòng)力學(xué)模型。剛?cè)狁詈隙囿w動(dòng)力學(xué)模型為基于模態(tài)截?cái)喾ǐ@得的降階有限元模型，假設(shè)有限元模型為

其中：M為有限元模型質(zhì)量矩陣；K為有限元模型剛度矩陣；üd和ud分別為加速度矢量和位移矢量；F為載荷矢量。選取模型受外載荷點(diǎn)和存在位移約束點(diǎn)為邊界點(diǎn)，將式中M、K的各行列重新整理，得到剛?cè)狁詈蟿?dòng)力學(xué)模型

其中，下標(biāo)b 代表邊界點(diǎn)，i 代表內(nèi)部節(jié)點(diǎn)。通過Craig-Bampton 方法[13]，保持邊界點(diǎn)6 自由度剛體位移，其余節(jié)點(diǎn)位移采用模態(tài)化處理，則剛?cè)狁詈蟿?dòng)力學(xué)模型的位移矢量轉(zhuǎn)化為

其中：η為邊界點(diǎn)固定時(shí)的模態(tài)自由度位移；H為Craig-Bampton 轉(zhuǎn)化矩陣，

其中：Ξ為邊界點(diǎn)約束模態(tài)，由有限元模型剛體模態(tài)給定；Φ為固定邊界點(diǎn)模態(tài)陣型，由有限元固定邊界點(diǎn)后的彈性模態(tài)給定。

假設(shè)動(dòng)力學(xué)分析頻率最高為fmax，在H矩陣中對(duì)Φ陣型展開模態(tài)截?cái)?，僅保留模態(tài)頻率低于fmax的模態(tài)，則剛?cè)狁詈蟿?dòng)力學(xué)模型的位移矢量轉(zhuǎn)化為

其中： η*為降階模態(tài)自由度位移；H*為降階Craig-Bampton 轉(zhuǎn)化矩陣。最終降階子結(jié)構(gòu)剛?cè)狁詈蟿?dòng)力學(xué)模型為

可獲得降階質(zhì)量和降階剛度的具體數(shù)值。

令

可將子結(jié)構(gòu)動(dòng)力學(xué)方程轉(zhuǎn)化為狀態(tài)空間方程的形式

其中：

1.2 舵軸間隙碰撞模型

建立多個(gè)子結(jié)構(gòu)模型狀態(tài)方程，以頂桿子結(jié)構(gòu)X和舵軸撥片子結(jié)構(gòu)Y為例，它們相互耦合時(shí)的位移及受力關(guān)系為

其中：下標(biāo)u、ü代表由位移、加速度引起的內(nèi)力；g為非線性映射關(guān)系；σ和 σ˙為子結(jié)構(gòu)模型的相對(duì)位移和相對(duì)速度。

彈性接觸問題采用罰函數(shù)模擬[14]。假設(shè)頂桿與舵軸的初始間隙為δ，頂桿與舵軸撥片發(fā)生碰撞的判斷式為

當(dāng)頂桿與舵軸撥片的相對(duì)位移滿足該判斷條件，則存在碰撞彈性力Fs與碰撞阻尼力Fd：

其中最大碰撞彈性力約束σmax和最大碰撞阻尼力約束 σ ˙max為人工定義參數(shù)，碰撞彈性力作用于頂桿與舵軸雙模型上，且方向與碰撞點(diǎn)相反；而碰撞阻尼力僅作用于發(fā)生碰撞時(shí)相對(duì)速度較大的模型中。根據(jù)頂桿-舵軸撥片間隙碰撞模型，可由頂桿-舵軸撥片的相對(duì)位移σ和相對(duì)速度 σ˙，獲得碰撞彈性力Fs與碰撞阻尼力Fd。

1.3 伺服舵機(jī)機(jī)電耦合模型

采用傳遞函數(shù)形式建立舵機(jī)機(jī)電模型，舵機(jī)輸入電壓信號(hào)V為伺服系統(tǒng)給定信號(hào)，負(fù)載轉(zhuǎn)速ωB（減速齒輪的主動(dòng)齒輪轉(zhuǎn)速）與負(fù)載功率W的關(guān)系式為

其中，扭矩常數(shù)kt、反電動(dòng)勢(shì)常數(shù)kv和電阻R均為給定參數(shù)。

舵機(jī)系統(tǒng)傳動(dòng)模型采用傳遞函數(shù)形式建模，其中包含減速齒輪、滾珠絲杠和頂桿裝置。包含摩擦影響的減速齒輪模型輸出為

其中：ωB、ωF分別為主動(dòng)、從動(dòng)齒輪轉(zhuǎn)速；rB、rF分別為主動(dòng)、從動(dòng)齒輪半徑。從動(dòng)齒輪的輸出功率WF為

其中Wloss為摩擦損耗功率，

其中η為摩擦損耗系數(shù)，由人工指定。

滾珠絲杠模型輸出為

其中：v為絲杠的平動(dòng)速度；L為絲杠長度。假設(shè)滾珠絲杠為理想模型，則絲杠的輸出功率與從動(dòng)齒輪輸出功率WF相同。

同樣假設(shè)頂桿裝置為理想模型且與絲杠連接處不存在間隙摩擦，則頂桿裝置模型輸出速度為v，輸出功率與絲杠輸出功率WF相同，此時(shí)頂桿的輸出力為

舵機(jī)一般采用電流環(huán)、速度環(huán)、位置環(huán)三閉環(huán)控制，如圖1 所示，電流環(huán)居內(nèi)，速度環(huán)居中，位置環(huán)為最外環(huán)。電流環(huán)輸入為舵機(jī)內(nèi)電流，速度環(huán)輸入為舵機(jī)輸出速度，位置環(huán)輸入為舵軸偏轉(zhuǎn)角度，根據(jù)PID 控制算法調(diào)整相應(yīng)控制參數(shù)，可令舵機(jī)伺服特性滿足需求。

圖1 伺服系統(tǒng)電機(jī)控制程序示意Fig.1 Diagram of the motor control program for servo system

1.4 集中氣動(dòng)力模型

基于偶極子網(wǎng)格法和最小狀態(tài)法有理函數(shù)近似擬合獲得[15]廣義氣動(dòng)力模型的一般表達(dá)式

其中：s為拉普拉斯域變量；A0代表氣動(dòng)剛度；A1代表氣動(dòng)阻尼；A2代表非定常氣動(dòng)力表觀質(zhì)量；D、R和E為氣動(dòng)力響應(yīng)系數(shù)矩陣；I為單位矩陣。根據(jù)1.1 節(jié)多體動(dòng)力學(xué)模型中邊界點(diǎn)選取獲得邊界點(diǎn)前n階模態(tài)陣型為Ψ，其維度為n×N，其中N為舵面所受集中氣動(dòng)力個(gè)數(shù)。則舵面響應(yīng)位置上所受集中氣動(dòng)力為

2 全動(dòng)舵系統(tǒng)模型

本章驗(yàn)證采用剛?cè)狁詈辖惦A模型建立典型舵面動(dòng)力學(xué)模型的正確性，首先開展剛?cè)狁詈辖惦A模型的靜載荷激勵(lì)下響應(yīng)與商用有限元軟件計(jì)算結(jié)果對(duì)比。

典型舵面的幾何尺寸如圖2 所示，據(jù)此建立的幾何模型如圖3 所示，其有限元模型如圖4 所示，有限元單元個(gè)數(shù)為1375，節(jié)點(diǎn)個(gè)數(shù)為1460。

圖2 典型舵面模型幾何尺寸Fig.2 Geometric sizes of the typical rudder surface

圖3 典型舵面幾何模型Fig.3 Geometric model of the typical rudder surface

圖4 典型舵面有限元模型Fig.4 Finite element model of the typical rudder surface

根據(jù)有限元模型，固定節(jié)點(diǎn)后求解典型舵面模態(tài)振型，保留其前6 階模態(tài)，構(gòu)造典型舵面的剛?cè)狁詈辖惦A模型。固定該模型舵軸處的自由度，給定舵面翼尖后緣靜載荷激勵(lì)，獲得典型舵面翼尖后緣處節(jié)點(diǎn)的時(shí)域響應(yīng)，并與商用有限元軟件計(jì)算結(jié)果進(jìn)行對(duì)比。為了提升計(jì)算效率并降低數(shù)值求解發(fā)散的概率，經(jīng)不同算法對(duì)比后，本文基于剛?cè)狁詈辖惦A模型的動(dòng)力學(xué)數(shù)值算法均采用MatLab 中的變步長ode23t 算法。對(duì)比結(jié)果如圖5 所示，基于Craig-Bampton 方法建立的典型舵面剛?cè)狁詈辖惦A模型與基于商用軟件建立的有限元模型之動(dòng)力學(xué)響應(yīng)間的偏差較小，證明了采用基于Craig-Bampton 方法建立典型舵面降階模型的正確性。

圖5 無阻尼典型舵面靜載荷下時(shí)域響應(yīng)的計(jì)算結(jié)果對(duì)比Fig.5 Comparison of calculation results for the time domain response of a typical rudder surface under static load with no damping

典型舵面邊界點(diǎn)取6 個(gè)，以翼根前緣為坐標(biāo)原點(diǎn)，6 個(gè)節(jié)點(diǎn)的坐標(biāo)分別為(0, 0, 0)(550, 0, 0)(200,250, 0)(550, 250, 0)(280, -80, 0)(280, -40, 20)，其中：節(jié)點(diǎn)1～4 為激振器激勵(lì)節(jié)點(diǎn)，連接方式采用多體動(dòng)力學(xué)連接模型；節(jié)點(diǎn)5 為舵軸撥片節(jié)點(diǎn)，與頂桿連接，開展無間隙仿真時(shí)連接方式采用線性多體動(dòng)力學(xué)連接模型，開展間隙仿真時(shí)采用頂桿-舵軸撥片間隙碰撞模型；節(jié)點(diǎn)6 為舵軸根部節(jié)點(diǎn)，與支撐結(jié)構(gòu)相連，連接方式采用多體動(dòng)力學(xué)連接模型。

3 數(shù)值算例

3.1 全動(dòng)舵面顫振結(jié)果分析

建立典型舵面計(jì)及阻尼的剛?cè)狁詈辖惦A模型，其阻尼由如下比例給定，

其中，α取值為0，β取值為0.001。計(jì)及阻尼后，典型舵面基于偶極子網(wǎng)格法計(jì)算獲得的顫振速度為1256 m/s，基于本文建立的顫振模型在不同來流速度下的時(shí)域響應(yīng)如圖6 所示，仿真獲得顫振速度為1270 m/s，兩者間的偏差僅為1.11%，證明基于剛?cè)狁詈隙囿w動(dòng)力學(xué)建立的顫振模型精度符合要求。

圖6 典型舵面在不同來流速度下的第二階廣義位移時(shí)域響應(yīng)Fig.6 Time domain response of the second-order generalized displacement of a typical rudder surface under different free stream velocities

3.2 計(jì)及舵軸間隙的全動(dòng)舵面顫振結(jié)果分析

舵軸轉(zhuǎn)動(dòng)自由度釋放自由，舵軸撥片與撥叉采用碰撞模型，撥叉保持固定，即舵機(jī)模型不參與耦合，碰撞距離δ設(shè)定為舵軸轉(zhuǎn)動(dòng)角度（0.38°），采用傳感器檢測舵軸角度θ和舵面4 節(jié)點(diǎn)位移x～，則舵面廣義位移由各節(jié)點(diǎn)的修正位移x與模態(tài)逆矩陣相乘獲得，各節(jié)點(diǎn)的修正位移為

其中l(wèi)為節(jié)點(diǎn)到舵軸沿舵面的位移。

基于計(jì)及舵軸間隙的全動(dòng)舵面顫振模型，典型舵面在不同來流速度下的時(shí)域響應(yīng)如圖7 和圖8所示，可以發(fā)現(xiàn)：存在間隙時(shí)，在來流速度遠(yuǎn)低于顫振速度的工況下，全動(dòng)舵面會(huì)產(chǎn)生結(jié)構(gòu)極限環(huán)振蕩現(xiàn)象（圖7(a)、圖8(a)）；在來流速度高于顫振速度的工況下，最初在短時(shí)間內(nèi)出現(xiàn)既不收斂也不發(fā)散的結(jié)構(gòu)極限環(huán)振蕩現(xiàn)象，隨后振動(dòng)幅值逐步增大直至振蕩發(fā)散（圖7(b)、圖8(b)）。

圖7 δ=0.38°時(shí)典型舵面在不同來流速度下的第二階廣義位移時(shí)域響應(yīng)Fig.7 Time domain response of the second-order generalized displacement of a typical rudder surface under different free stream velocities for δ=0.38°

圖8 δ=0.38°時(shí)典型舵面在不同來流速度下的舵軸轉(zhuǎn)角時(shí)域響應(yīng)Fig.8 Time domain response of rotating angle of a typical rudder surface under different free stream velocities for δ=0.38°

采用上述相同參數(shù)，碰撞距離δ設(shè)定改為0.95°時(shí)，基于計(jì)及舵軸間隙顫振模型，典型舵面在不同來流速度下的時(shí)域響應(yīng)如圖9 和圖10 所示，可以發(fā)現(xiàn)：此時(shí)來流速度高于顫振速度后出現(xiàn)既不收斂也不發(fā)散的結(jié)構(gòu)振蕩，且并未隨仿真時(shí)間的延長而出現(xiàn)振蕩發(fā)散現(xiàn)象，推測產(chǎn)生該現(xiàn)象的原因是仿真時(shí)間較短。

圖9 δ=0.95°時(shí)典型舵面在不同來流速度下的第二階廣義位移時(shí)域響應(yīng)Fig.9 Time domain response of the second-order generalized displacement of a typical rudder surface under different free stream velocities for δ=0.95°

圖10 δ=0.95°時(shí)典型舵面在不同來流速度下的舵軸轉(zhuǎn)角時(shí)域響應(yīng)Fig.10 Time domain response of rotating angle of a typical rudder surface under different free stream velocities for δ=0.95°

綜上可知，存在舵軸間隙時(shí)，全動(dòng)舵面會(huì)在來流速度較低時(shí)產(chǎn)生極限環(huán)振蕩現(xiàn)象，但對(duì)于顫振現(xiàn)象的發(fā)生具有減緩作用，且隨著間隙的增大，全動(dòng)舵面產(chǎn)生顫振發(fā)散現(xiàn)象的時(shí)刻會(huì)推后。

3.3 全動(dòng)舵系統(tǒng)顫振結(jié)果分析

隨后開展舵機(jī)系統(tǒng)對(duì)于全動(dòng)舵面顫振特性的研究，此時(shí)舵軸轉(zhuǎn)動(dòng)自由度釋放，舵軸撥片與撥叉采用多體動(dòng)力學(xué)連接模型，撥叉運(yùn)動(dòng)由舵機(jī)給定，舵機(jī)模型參與振動(dòng)耦合。同樣采用傳感器檢測舵軸角度θ和舵面4 節(jié)點(diǎn)位移x～，用于求解廣義位移（參見式(23)）。

基于全動(dòng)舵系統(tǒng)顫振模型的不同來流速度下第二階廣義位移時(shí)域響應(yīng)如圖11 和圖12 所示，分析發(fā)現(xiàn)：由于舵軸撥片未固定，采用伺服舵機(jī)驅(qū)動(dòng)導(dǎo)致某些時(shí)刻出現(xiàn)舵軸轉(zhuǎn)角偏離平衡位置較大距離的現(xiàn)象，并引發(fā)彈性振動(dòng)幅值急劇增加；但隨著舵軸轉(zhuǎn)角逐步至0°，彈性振動(dòng)的收斂和發(fā)散現(xiàn)象由來流速度決定。

圖11 全動(dòng)舵系統(tǒng)在不同來流速度下的第二階廣義位移時(shí)域響應(yīng)Fig.11 Time domain response of the second-order generalized displacement of an all-movable rudder under different free stream velocities

圖12 全動(dòng)舵系統(tǒng)在不同來流速度下的舵軸轉(zhuǎn)角時(shí)域響應(yīng)Fig.12 Time domain response of rotating angle of an allmovable rudder under different free stream velocities

仿真表明：在舵軸偏離平衡位置時(shí)，舵機(jī)伺服系統(tǒng)給定控制指令迫使舵軸偏轉(zhuǎn)角度逐步靠近0°，該控制效果由舵機(jī)伺服系統(tǒng)的控制參數(shù)給定，因此伺服控制系統(tǒng)中參數(shù)調(diào)整不合理時(shí)可能加劇顫振現(xiàn)象的發(fā)生，甚至降低顫振速度。

4 結(jié)束語

本文基于柔性多體動(dòng)力學(xué)方法建立全動(dòng)舵系統(tǒng)顫振模型，采用Craig-Bampton 方法建立典型舵面剛?cè)狁詈辖惦A模型，采用激振器式集中氣動(dòng)力模擬基于偶極子網(wǎng)格法獲得廣義氣動(dòng)力。數(shù)值仿真結(jié)果分析表明，該系統(tǒng)模型能夠快速與間隙、伺服系統(tǒng)等多種模型相互耦合，適用于研究各種因素對(duì)全動(dòng)舵面動(dòng)力學(xué)響應(yīng)的影響。本文主要研究舵軸間隙、舵機(jī)對(duì)全動(dòng)舵面顫振特性的影響，得到主要結(jié)論如下：

1）采用Craig-Bampton 方法建立全動(dòng)舵面動(dòng)力學(xué)模型與激振器式集中氣動(dòng)力模型耦合可獲得舵面顫振臨界速度，且與商用軟件預(yù)測對(duì)比的偏差小于2%；

2）全動(dòng)舵面的舵軸存在間隙時(shí)，發(fā)生極限環(huán)振蕩的臨界速度遠(yuǎn)遠(yuǎn)低于模型不存在間隙時(shí)的情況，且來流速度大于顫振速度后，短時(shí)間內(nèi)出現(xiàn)既不收斂也不發(fā)散的結(jié)構(gòu)極限環(huán)振蕩現(xiàn)象，且極限環(huán)振蕩持續(xù)時(shí)間隨著間隙增大而延長；

3）全動(dòng)舵系統(tǒng)的舵機(jī)影響將導(dǎo)致某些時(shí)刻出現(xiàn)舵軸轉(zhuǎn)角偏離平衡位置較大距離，由此引發(fā)彈性振動(dòng)幅值急劇增加，意味著伺服控制系統(tǒng)中參數(shù)調(diào)整不合理時(shí)將可能加劇結(jié)構(gòu)顫振，甚至降低顫振速度。