龔學(xué)兵，胡大慶，王永平，鄭 佩

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

0 引言

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

針對(duì)導(dǎo)彈、飛機(jī)的氣動(dòng)伺服彈性問(wèn)題，美國(guó)航空航天學(xué)會(huì)作了一系列有關(guān)氣動(dòng)彈性現(xiàn)象、氣動(dòng)彈性分析和氣動(dòng)彈性試驗(yàn)等方面的專(zhuān)題綜述，并對(duì)氣動(dòng)伺服彈性研究的未來(lái)發(fā)展趨勢(shì)提供指導(dǎo)[1]。目前，針對(duì)氣動(dòng)彈性問(wèn)題，主要從理論分析和數(shù)值仿真計(jì)算、試驗(yàn)分析等三方面開(kāi)展相關(guān)研究[2]。

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

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

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

1 未變形箭體的氣動(dòng)計(jì)算

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

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

圖2 未變形箭體的ICEMCF氣動(dòng)網(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

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

氣動(dòng)參數(shù)計(jì)算結(jié)果如表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 變形箭體的理論氣動(dòng)參數(shù)計(jì)算

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

2.1 彈性動(dòng)力系數(shù)的理論計(jì)算公式

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

(1)

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

(2)

(3)

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

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

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

(1)對(duì)于未變形的箭體

(2)對(duì)于變形后的箭體

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

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

(4)

其中，

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

2.2 模態(tài)分析

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

圖4 結(jié)構(gòu)圓孔尺寸保留示意圖Fig.4 Schematic diagram of structural with round hole size reservation

模型簡(jiǎn)化后，為了保證質(zhì)量分布符合實(shí)體結(jié)構(gòu)的質(zhì)量分布特征，必須在模型中加入實(shí)體模型的其他附屬件的質(zhì)量，質(zhì)量分布如圖5所示。在質(zhì)量分布設(shè)置過(guò)程中，應(yīng)該根據(jù)質(zhì)量分布以及剛度影響等因素，設(shè)置相應(yīng)的點(diǎn)質(zhì)量和面質(zhì)量。

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

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

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

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

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

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

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

2.4 基于理論分析的氣動(dòng)升力對(duì)比分析

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

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

3 變形箭體的氣動(dòng)對(duì)比分析

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

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

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

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

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

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

對(duì)比表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ù)條件下的阻力會(huì)逐漸減小；

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

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

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

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

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

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

表7 變形火箭一階振型修正的法向力系數(shù)斜率的理論計(jì)算結(jié)果與Fluent計(jì)算結(jié)果對(duì)比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)定性分析

通過(guò)式(2)～式(4)，結(jié)合表7給出的參數(shù)，即可完成彈性動(dòng)力學(xué)系數(shù)的計(jì)算。

表8 一級(jí)飛行橫向運(yùn)動(dòng)方程式系數(shù)Table 8 Coefficient of lateral motion equation of first flight

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

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

4 結(jié)論

通過(guò)針對(duì)某型火箭開(kāi)展的相關(guān)氣動(dòng)彈性仿真計(jì)算及結(jié)果分析，得出以下結(jié)論,為箭體結(jié)構(gòu)設(shè)計(jì)及氣動(dòng)彈性控制系統(tǒng)設(shè)計(jì)提供依據(jù)。

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

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

(3)通過(guò)伯德圖可以看出本文大長(zhǎng)細(xì)比火箭氣動(dòng)彈性控制系統(tǒng)設(shè)計(jì)能夠?qū)崿F(xiàn)對(duì)箭體彈性振動(dòng)的主動(dòng)抑制，確保火箭飛行箭體結(jié)構(gòu)的穩(wěn)定性。