龔學兵，胡大慶，王永平，鄭 佩

(中國航天科技集團有限公司四院四十一所，西安 710025)

0 引言

隨著火箭總體設計技術的快速發(fā)展，火箭的長細比越來越大，相對彎曲剛度越來越小，箭體的氣動載荷分布更加復雜，控制模型的非線性特征更強，火箭氣動彈性控制問題越來越突出。氣動彈性控制技術是未來火箭總體設計的關鍵技術之一，該項技術能有效提升火箭控制精度、降低結構的消極質量、提升彈道規(guī)劃的能力，有利于分析總體參數(shù)拉偏范圍對產(chǎn)品性能的影響以及飛行試驗的危險因素來源。

針對導彈、飛機的氣動伺服彈性問題，美國航空航天學會作了一系列有關氣動彈性現(xiàn)象、氣動彈性分析和氣動彈性試驗等方面的專題綜述，并對氣動伺服彈性研究的未來發(fā)展趨勢提供指導[1]。目前，針對氣動彈性問題，主要從理論分析和數(shù)值仿真計算、試驗分析等三方面開展相關研究[2]。

理論分析方面，通過在飛行器動力學方程的基礎上引入結構動力學模型，并將氣動力通過有理函數(shù)(振型函數(shù))進行廣義非定常氣動力進行擬合[3-7]，從而實現(xiàn)了通過理論的數(shù)學模型研究柔性變形、橫向彎曲振動、旋轉鎖定等氣動彈性問題的機理，并研究彈性變形對箭體姿態(tài)、控制系統(tǒng)設計的影響；而針對箭體姿態(tài)運動與彈性變形的耦合作用，部分作者[8-10]通過引入模態(tài)化方程組、大運動彈性體線性化處理、非均勻梁彎曲運動方程組等數(shù)學模型實現(xiàn)擺動與彈體振動耦合分析、慣量矩陣的時變特征分析、以及彈性導彈的側向運動分析。雖然上述方法能夠從理論角度給出彈性運動的精確結果，但非線性特征更加突出，導致控制系統(tǒng)設計存在困難。

數(shù)值仿真計算方面，采用計算流體力學程序計算非線性非定常氣動力是當前跨聲速氣動伺服彈性系統(tǒng)建模研究的主要手段。氣動伺服彈性的非線性特征導致流固耦合計算陷入開環(huán)系統(tǒng)模型階數(shù)過高的問題，不利于控制率的設計。目前,相關學者[11-13]開展了相關模型降階、提升計算效率、增強算法的魯棒性和提供求解精度的相關學術研究。氣動仿真[14-17]通過數(shù)字求解，能將復雜的模型通過分段線性化插值的手段，將非線性問題轉為線性化問題，有利于準確的分析彈性變形對箭體的升力系數(shù)、阻力系數(shù)和俯仰力矩系數(shù)的影響，以及在不同馬赫數(shù)和攻角下的氣動力變化規(guī)律，同時進行結構強度校核與系統(tǒng)顫振速度分析，但上述方法的定量結果難以通過線性化的數(shù)學模型進行表征。

綜上所述，雖然國、內外研究學者在氣動彈性控制問題領域開展了大量的理論與仿真分析工作，但上述研究未從總體設計角度考慮，難以在線性化處理后的氣動伺服彈性動力學模型的基礎上，為控制系統(tǒng)設計提供氣動參數(shù)、氣動彈性系數(shù)的準確拉偏范圍。工程上常用的氣動伺服彈性分析方法是利用線性化的傳遞函數(shù)模型，分析彈箭開環(huán)系統(tǒng)的幅頻、相頻特性，開展控制系統(tǒng)的頻域PID控制系統(tǒng)設計，通過添加陷幅濾波器抑制彈箭開環(huán)系統(tǒng)與飛行控制系統(tǒng)的不利耦合作用。而彈性箭體的氣動參數(shù)具有非線性特征，不利于控制系統(tǒng)設計。本文通過理論分析部分提供控制系統(tǒng)所需的線性化彈性動力系數(shù)，并通過數(shù)值仿真技術為控制系統(tǒng)參數(shù)拉偏提供準確的變化范圍，結合伯德圖的穩(wěn)定性分析結果，為氣動彈性控制的總體設計提供指導。

1 未變形箭體的氣動計算

某火箭的箭體結構與氣動外形見圖1，長細比達到16，氣動彈性問題較為突出。在飛行過程中，各點的彈道和氣動力參數(shù)均處于不斷的變化之中，在氣動彈性參數(shù)計算時不可能對全彈道各點的參數(shù)都進行考察。因此，本文給出了火箭在飛行過程中攻角取0°、2°、4°和馬赫數(shù)分別為2、3、4、5狀態(tài)組合下的氣動參數(shù)。

利用ICEM CFD對上述理論化模型進行網(wǎng)格劃分。其中，尾部流場區(qū)域選取箭體參考長度的7.5倍左右，頭部區(qū)域選取箭體參考長度3.75倍左右，徑向外流場取箭體參考長度的3倍左右(圖2)。在Fluent中設置了壁面函數(shù)，在劃分網(wǎng)格之前根據(jù)相似型號的工程經(jīng)驗確定第一層網(wǎng)格的厚度值(商業(yè)上可以用Pointwise與CFD-online計算yplus值，從而確定第一層網(wǎng)格厚度值。k-ε模型的yplus值為30～300之間)。本文的首層網(wǎng)格的厚度值設置為1 mm，附面層的厚度為15 mm，網(wǎng)格的厚度增長率取1.1(圖3)。網(wǎng)格的無關性是通過對比多個設計師的仿真計算結果，確認網(wǎng)格數(shù)量和網(wǎng)格質量不影響氣動參數(shù)的計算精度。

圖2 未變形箭體的ICEMCF氣動網(wǎng)格劃分Fig.2 Aerodynamic grid generation of un- deformed rocket by ICEMCFD

圖3 尾翼位置的局部網(wǎng)格及邊界層網(wǎng)格示意圖Fig.3 Aerodynamic grid generation of un- deformed rocket by ICEMCFD

氣動計算的控制方程為雷諾平均Navier-Stokes(RANS)方程，空間離散格式為一階迎風格式，采用雙時間步推進進行求解。通過k-ε湍流模型和壁面函數(shù)法求解近壁面處邊界層，遠場邊界條件為壓力遠場，箭體壁面邊界條件為無滑移壁面。在開展箭體的氣動參數(shù)仿真計算時，單個狀態(tài)的計算需要迭代到5000步。針對高超聲速狀態(tài)，需要選取不同的Courant數(shù)，其默認值是5。在實際計算中，初始Courant設置為1，逐步增加為3、5、10，每次增加Courant參數(shù)后，在上一步計算結果的基礎上重新進行氣動仿真計算。每個Case的工況由馬赫數(shù)、大氣壓共同決定。

氣動參數(shù)計算結果如表1所示。

表1 未變形的軸向力系數(shù)Ca1、向力系數(shù)Cn1及俯仰力矩系數(shù)Mz1Table 1 Axial force coefficient(Ca1)，normal force coefficient(Cn1)and pitching moment coefficient(Mz1) of un-deformed rocket

2 變形箭體的理論氣動參數(shù)計算

變形箭體的理論計算為控制系統(tǒng)設計提供線性化的理論模型，便于控制系統(tǒng)通過經(jīng)典PID控制系統(tǒng)設計方法進行氣動彈性控制系統(tǒng)設計。

2.1 彈性動力系數(shù)的理論計算公式

假設火箭彈性振動用廣義坐標qi(t)表示，它是隨時間變化的量，由下列二階常微分方程確定：

(1)

式中ωi為火箭第i次振型的固有頻率；ξi為火箭第i次振型的阻尼系數(shù)；Qi為火箭對應的第i次振型的廣義力；Mi為火箭對應的第i次振型的廣義質量。

(2)

(3)

式中fy1(x)為作用于火箭的外力在箭體軸y1上的投影；m(x)為沿箭體縱軸的彈體質量分布。

頻率ωi和固有振型函數(shù)Wi(x)為殼體結構特性和沿縱軸質量分布的函數(shù)；Qi為影響彈體振動的廣義力。由式(1)可知，只要知道fy1(x,t)即可求出其表達式，fy1(x,t)包括控制力、氣動力、推力以及舵的慣性力等諸力沿箭體軸oy1的投影。

在未變形的氣動參數(shù)基礎上，進行變形箭體的氣動參數(shù)計算。通過振型函數(shù)、模態(tài)質量進行氣動參數(shù)的重新計算，主要依據(jù)如下：

(1)對于未變形的箭體

(2)對于變形后的箭體

除氣動力外，上面提到的各力都作用在火箭箭體的某一位置，所以是集中力，由(1)和(2)可直接得到相應的廣義力。

將廣義力的具體表達式代入式(4)，略去各次振型之間的相互影響，即得到變形火箭在某次振型影響下的振動方程：

(4)

其中，

通過式(4)結合彈性振動的彈體小擾動方程，獲取彈性振動下舵偏角到彈體姿態(tài)角的傳遞函數(shù)，結合伯德圖的幅值裕度和相角裕度，即可完成彈性變形對火箭飛行穩(wěn)定性的影響。

2.2 模態(tài)分析

模態(tài)分析主要用于提取振型，對實體幾何模型優(yōu)化，要求對大的圓孔進行保留，確保變形斜率的準確性。圖4給出了實體結構模型的大尺寸的圓孔等特征保留的示意圖。

圖4 結構圓孔尺寸保留示意圖Fig.4 Schematic diagram of structural with round hole size reservation

模型簡化后，為了保證質量分布符合實體結構的質量分布特征，必須在模型中加入實體模型的其他附屬件的質量，質量分布如圖5所示。在質量分布設置過程中，應該根據(jù)質量分布以及剛度影響等因素，設置相應的點質量和面質量。

圖5 箭體質量分布示意圖Fig.5 Diagram of mass distribution of rocket

全箭模態(tài)仿真結果如表2所示，同一階彎曲模態(tài)包含兩個正交方向上的模態(tài)振型。由于Y、Z方向的質量分布存在差異，導致不同象限(I、III與II、IV)的頻率分布存在一定的差異性，該頻率計算結果與相同型號產(chǎn)品地面試驗結果修正后的頻率偏差在2 Hz范圍以內。通常振型斜率應由速率陀螺進行測量，但由于試驗條件限制，當前只能以理論分析為主。

表2 箭體模態(tài)仿真結果Table 2 Modal analysis results of rocket body Hz

2.3 振型函數(shù)的提取

根據(jù)上述結構模態(tài)的分析結果，采用APDL語言，提取彈體坐標位置以及相應的振型幅值，圖6為箭體振型隨著彈體位置變化的振型幅值函數(shù)圖。坐標原點選取彈體頂點作為參考點。

(a) Ma=2 (b) Ma=3

(c) Ma=4 (d) Ma=5圖6 振型函數(shù)Fig.6 Modal shape function

2.4 基于理論分析的氣動升力對比分析

表3 一階振型修正的法向力系數(shù)斜率Table 3 Slope of normal force revised by first bend mode

對比分析計算結果可知，由于馬赫數(shù)的增加，導致理論計算的彈體氣動升力斜率逐漸加大。這是由于馬赫數(shù)的增加，導致一階振型的幅值變大，振型引起的附加攻角增大，變形箭體的法向力系數(shù)斜率逐漸增大。

3 變形箭體的氣動對比分析

雖然上述理論分析可以用于氣動彈性控制的動力系數(shù)計算，但變形后的箭體氣動特性與未變形箭體的氣動特性存在一定的差異性，例如：彈體壓心前移、氣動系數(shù)不對稱等，而且線性的方程與非線性氣動彈性方程求解也存在一定的誤差，需要開展風洞試驗進行修正，而風洞試驗費用昂貴。本文通過Fluent仿真模擬試驗為理論計算結果提供拉偏范圍，用于降低氣動伺服彈性控制系統(tǒng)的設計風險，并給控制系統(tǒng)的彈道拉偏提供合理的變化范圍。

根據(jù)彈體變形的結果，獲得彈體一階變形的振型函數(shù)。忽略彈翼、舵面的變形。變形后的幾何模型如圖7所示。

圖7 變形后的箭體模型Fig.7 Geometric model of deformed rocket

根據(jù)變形后的彈體通過ICEM CFD軟件劃分變形火箭流場網(wǎng)格，如圖8所示。

圖8 變形箭體的氣動網(wǎng)格劃分結果Fig.8 Aerodynamic meshing result of the deformed rocket

計算獲得火箭在一階變形后的氣動阻力、升力、俯仰力矩系數(shù)(表4)，并與理論計算結果進行對比。

對比表1和表4可以看出：

表4 變形的軸向力系數(shù)Ca1、向力系數(shù)Cn1及俯仰力矩系數(shù)Mz1Table 4 Axial force coefficient(Ca1)，normal force coefficient(Cn1)and pitching moment coefficient(Mz1) of deformed rocket

(1)彈體在同一馬赫數(shù)條件下的阻力會逐漸減??；

(2)零度攻角下彈體的升力不為零，這是由于附加攻角導致彈體在理論0°攻角下仍然存在升力；

(3)零度攻角下彈體的俯仰力矩不為零，這是由于附加攻角導致彈體在理論0°攻角下仍然存在升力，進而產(chǎn)生俯仰力矩。

箭體載荷設計一般考慮極限條件下，本文選取4°攻角作為載荷設計依據(jù)，如表5所示。由表5可以看出，由于彈性變形，箭體的軸向力矩系數(shù)變小，而法相力矩系數(shù)和俯仰力矩系數(shù)變大。最大偏差分別為-15.59%、10.34%和16.91%。在箭體載荷設計時應考慮拉彎組合效應。同時在彈道拉偏過程中，應考慮箭體彈性變形引起的氣動參數(shù)變化。

此外，表6為變形前后的壓心系數(shù)，通過對比可以看出彈體壓心向前移動。隨著馬赫數(shù)增加，在4°攻角條件下，彈體壓心前移最大為13.4%，壓心的變化對控制系統(tǒng)的穩(wěn)定性設計應著重考慮。

表6 變形前、后的壓心系數(shù)Xp(α)Table 6 Pressure center coefficient(Xp(α))of undeformed rocket and deformed rocket

計算獲得彈體在一階變形后的氣動參數(shù)，并與未變形的氣動參數(shù)、氣彈理論計算的結果對比分析，如表7所示?？梢钥闯觯赡B(tài)疊加法獲得的彈體升力斜率隨著馬赫數(shù)增加而變小，理論計算結果與Fluent仿真計算結果在Ma=5條件下出現(xiàn)較大的偏差，達到14%。因此，在彈性動力系數(shù)計算時，應提供14%左右的拉偏量。

表7 變形火箭一階振型修正的法向力系數(shù)斜率的理論計算結果與Fluent計算結果對比Table 7 Slope of normal force coefficient revised by first bend mode of the deformed rocket with the theoretical and Fluent results

4 靜氣彈穩(wěn)定性分析

通過式(2)～式(4)，結合表7給出的參數(shù)，即可完成彈性動力學系數(shù)的計算。

表8 一級飛行橫向運動方程式系數(shù)Table 8 Coefficient of lateral motion equation of first flight

在彈性動力系數(shù)計算的基礎上，結合彈體的傳遞函數(shù)，即可通過伯德圖開展氣彈性穩(wěn)定裕度分析(圖9和表9)。由圖9可見，一階振型頻率對應的幅值曲線處于0 dB以下；由表9可見，一階彈性幅值裕度大于0 dB以上，滿足彈性運動幅值穩(wěn)定，系統(tǒng)相位滯后也在可以接受的范圍內。通過引入陷波濾波器，實現(xiàn)對箭體彈性振動的主動抑制。在拉偏下限狀態(tài)下，系統(tǒng)依然能夠滿足彈性幅值運動穩(wěn)定。

表9 伯德圖給出的相位、幅值裕度值Table 9 Amplitude and phase stability margin of bode diagram

4 結論

通過針對某型火箭開展的相關氣動彈性仿真計算及結果分析，得出以下結論,為箭體結構設計及氣動彈性控制系統(tǒng)設計提供依據(jù)。

(1)利用模態(tài)分析提供的振型函數(shù)，開展了箭體變形后的氣動特性分析。與未變形狀態(tài)相比，軸向力系數(shù)、法向力系數(shù)、俯仰力矩系數(shù)拉偏分別為-15.59%、10.34%和16.91%，此外箭體壓心最大前移13.4%。

(2)通過對比理論升力斜率與Fluent仿真計算的升力斜率，氣動彈性系數(shù)D1i、D2i、D3i的拉偏量級分別為20.4%、16.4%、15.2%。

(3)通過伯德圖可以看出本文大長細比火箭氣動彈性控制系統(tǒng)設計能夠實現(xiàn)對箭體彈性振動的主動抑制，確保火箭飛行箭體結構的穩(wěn)定性。