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

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

0 引言

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

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

為了能夠較為簡單且有效地耦合舵機模型、舵軸間隙等因素的影響，本文采用剛柔耦合多體動力學方法建立典型舵面動力學模型，并采用地面顫振試驗的集中式氣動力模擬方法開展氣動彈性耦合仿真，分析舵機模型、舵軸間隙等因素對動力學響應的影響。

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

1.1 剛柔耦合多體動力學模型

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

其中：M為有限元模型質量矩陣；K為有限元模型剛度矩陣；üd和ud分別為加速度矢量和位移矢量；F為載荷矢量。選取模型受外載荷點和存在位移約束點為邊界點，將式中M、K的各行列重新整理，得到剛柔耦合動力學模型

其中，下標b 代表邊界點，i 代表內部節(jié)點。通過Craig-Bampton 方法[13]，保持邊界點6 自由度剛體位移，其余節(jié)點位移采用模態(tài)化處理，則剛柔耦合動力學模型的位移矢量轉化為

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

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

假設動力學分析頻率最高為fmax，在H矩陣中對Φ陣型展開模態(tài)截斷，僅保留模態(tài)頻率低于fmax的模態(tài)，則剛柔耦合動力學模型的位移矢量轉化為

其中： η*為降階模態(tài)自由度位移；H*為降階Craig-Bampton 轉化矩陣。最終降階子結構剛柔耦合動力學模型為

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

令

可將子結構動力學方程轉化為狀態(tài)空間方程的形式

其中：

1.2 舵軸間隙碰撞模型

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

其中：下標u、ü代表由位移、加速度引起的內力；g為非線性映射關系；σ和 σ˙為子結構模型的相對位移和相對速度。

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

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

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

1.3 伺服舵機機電耦合模型

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

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

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

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

其中Wloss為摩擦損耗功率，

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

滾珠絲杠模型輸出為

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

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

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

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

1.4 集中氣動力模型

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

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

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

本章驗證采用剛柔耦合降階模型建立典型舵面動力學模型的正確性，首先開展剛柔耦合降階模型的靜載荷激勵下響應與商用有限元軟件計算結果對比。

典型舵面的幾何尺寸如圖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)，構造典型舵面的剛柔耦合降階模型。固定該模型舵軸處的自由度，給定舵面翼尖后緣靜載荷激勵，獲得典型舵面翼尖后緣處節(jié)點的時域響應，并與商用有限元軟件計算結果進行對比。為了提升計算效率并降低數(shù)值求解發(fā)散的概率，經不同算法對比后，本文基于剛柔耦合降階模型的動力學數(shù)值算法均采用MatLab 中的變步長ode23t 算法。對比結果如圖5 所示，基于Craig-Bampton 方法建立的典型舵面剛柔耦合降階模型與基于商用軟件建立的有限元模型之動力學響應間的偏差較小，證明了采用基于Craig-Bampton 方法建立典型舵面降階模型的正確性。

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

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

3 數(shù)值算例

3.1 全動舵面顫振結果分析

建立典型舵面計及阻尼的剛柔耦合降階模型，其阻尼由如下比例給定，

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

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

3.2 計及舵軸間隙的全動舵面顫振結果分析

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

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

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

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

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

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

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

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

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

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

圖11 全動舵系統(tǒ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)在不同來流速度下的舵軸轉角時域響應Fig.12 Time domain response of rotating angle of an allmovable rudder under different free stream velocities

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

4 結束語

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

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

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

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