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

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

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

令

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

其中：

1.2 舵軸間隙碰撞模型

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

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

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

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

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

1.3 伺服舵機機電耦合模型

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

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

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

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

其中Wloss為摩擦損耗功率，

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

滾珠絲杠模型輸出為

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

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

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

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

1.4 集中氣動力模型

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

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

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

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

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

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

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

3 數(shù)值算例

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

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

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

圖6 典型舵面在不同來流速度下的第二階廣義位移時域響應(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 計及舵軸間隙的全動舵面顫振結(jié)果分析

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

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

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

圖7 δ=0.38°時典型舵面在不同來流速度下的第二階廣義位移時域響應(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°時典型舵面在不同來流速度下的舵軸轉(zhuǎn)角時域響應(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°時，基于計及舵軸間隙顫振模型，典型舵面在不同來流速度下的時域響應(yīng)如圖9 和圖10 所示，可以發(fā)現(xiàn)：此時來流速度高于顫振速度后出現(xiàn)既不收斂也不發(fā)散的結(jié)構(gòu)振蕩，且并未隨仿真時間的延長而出現(xiàn)振蕩發(fā)散現(xiàn)象，推測產(chǎn)生該現(xiàn)象的原因是仿真時間較短。

圖9 δ=0.95°時典型舵面在不同來流速度下的第二階廣義位移時域響應(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°時典型舵面在不同來流速度下的舵軸轉(zhuǎn)角時域響應(yīng)Fig.10 Time domain response of rotating angle of a typical rudder surface under different free stream velocities for δ=0.95°

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

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

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

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

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

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

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

4 結(jié)束語

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

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

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

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