鄧林嘯 呂竹勇 王天舒 呂敬?

（1.北京航空航天大學(xué)航空科學(xué)與工程學(xué)院，北京 100191）（2.北京空天技術(shù)研究所，北京 100074）（3.清華大學(xué)航空航天學(xué)院 北京 100084）

引言

在飛行器中，除了部分導(dǎo)彈和少數(shù)運載火箭使用固體燃料以外，其余大部分的飛行器包括航天器和航空器都使用液體燃料.而隨著飛行器航程和任務(wù)復(fù)雜度的增加，液體燃料在飛行器整體質(zhì)量上的占比越來越高，則自由液面所導(dǎo)致的液體晃動對飛行器的影響也越來越不可忽視[1].從上世紀(jì)60年代起，相關(guān)的研究人員就開始研究儲箱內(nèi)的液體燃料對飛行器的姿態(tài)控制、飛行性能和穩(wěn)定性的影響[2～4].

當(dāng)前的研究中，對于飛行器整體耦合動力學(xué)建模過程中的晃動力計算通常采用CFD、SPH[5]、運動脈動球模型[6]、等效力學(xué)模型[7]等幾種方式.但是受限于計算資源和內(nèi)存，常用的CFD、SPH方法等在對充液儲箱內(nèi)的液體流場進行計算時用時太長，不利于整體耦合動力學(xué)方程的計算，所以通?？梢岳煤唵蔚牡刃ЯW(xué)模型來代替復(fù)雜的晃動流場計算.岳寶增等[4]利用廣義準(zhǔn)坐標(biāo)下的拉格朗日方法得到航天器剛體部分運動和液體燃料晃動的耦合動力學(xué)方程，采用復(fù)合控制方法對航天器姿態(tài)進行控制.張詩琪等[5]使用SPH方法計算晃動力，然后引入剛體約束方程得到固液耦合的動力學(xué)方程，分析了在周期激勵情況下的固液耦合特性.

對于飛行器姿態(tài)機動過程中流固耦合問題，目前多數(shù)學(xué)者都聚焦于真空微重環(huán)境下對航天器的研究，鮮少有對大氣中飛行的航空器的分析，故本文研究大氣飛行環(huán)境下液體晃動對飛行器控制系統(tǒng)的影響.

本文針對穩(wěn)態(tài)運動中的充液飛行器，此時液體晃動為小幅晃動，故而利用等效力學(xué)模型[8]描述液體晃動，將液體晃動問題轉(zhuǎn)化為多剛體運動動力學(xué)問題；在利用虛功率原理推導(dǎo)出多貯箱充液系統(tǒng)動力學(xué)方程后，結(jié)合線性小擾動理論，建立起充液飛行器線性小擾動方程組，并推導(dǎo)出飛行器縱向模態(tài)傳遞函數(shù).最后分析了橫放圓柱貯箱尺寸和充液比對飛行器縱向典型傳函系數(shù)的影響.

1 充液系統(tǒng)的增量方程

充液貯箱在小幅晃動情況下時，可以根據(jù)等效原則將液體晃動等效為等效擺力學(xué)模型，參考文獻(xiàn)[8]中液體晃動等效力學(xué)模型推導(dǎo)方法，利用虛功率原理推導(dǎo)出攜帶nf個貯箱充液飛行器系統(tǒng)動力學(xué)方程和第i階液體晃動方程如下：

上述中，Rc為飛行器在慣性系下位移，ωc為飛行器繞機體系轉(zhuǎn)動角速度，q為單擺相對飛行器機體系擺角；m、S、J為系統(tǒng)質(zhì)量、系統(tǒng)相對飛行器機體系原點靜矩、轉(zhuǎn)動慣量矩陣，Spi、Jpi為第i階單擺相對單擺系原點靜矩、轉(zhuǎn)動慣量；g為重力加速度；ωpi為第i階單擺相對飛行器機體系角速度；F、T為飛行器所受氣動力和氣動力矩；Qt、Qr分別為液體晃動導(dǎo)致的平動耦合力、轉(zhuǎn)動耦合力；rci、為第i階單擺相對單擺懸掛點矢量；Hpi為第i階單擺擺動面法線方向；SpHi、JpHi、Jrsi為與晃動擺相關(guān)項具體表達(dá)式可參考文獻(xiàn)[8].本文所用各坐標(biāo)系示意圖及其相互轉(zhuǎn)換關(guān)系如圖1所示，各坐標(biāo)系定義如下：

圖1 各坐標(biāo)系關(guān)系示意圖Fig.1 The diagram of the relationship of each coordinate system

地面坐標(biāo)系Oxyz：坐標(biāo)系原點取在飛行器發(fā)射點上，x軸沿彈道面與地面交線指向目標(biāo)點處，z軸沿鉛垂面向下，y軸與其它兩軸垂直并構(gòu)成右手系.

機體坐標(biāo)系Ox1y1z1：原點在飛行器質(zhì)心，Ox1軸與飛行器縱軸重合，指向頭部為正，Oy1軸位于飛行器橫向?qū)ΨQ面內(nèi)與Ox1軸垂直，Oz1軸垂直于另外兩軸，方向按右手系確定.

彈道坐標(biāo)系Ox2y2z2：原點在飛行器瞬時質(zhì)心上，Ox2軸與導(dǎo)彈的速度矢量V重合，Oy2軸垂直于包含速度矢量V的鉛垂面，Oz2軸與其他兩軸垂直構(gòu)成右手坐標(biāo)系并位于包含速度矢量的鉛垂面內(nèi).

速度坐標(biāo)系Ox3y3z3：坐標(biāo)系原點取在飛行器質(zhì)心上，Ox3軸與飛行器質(zhì)心的速度矢量V重合，Oy3軸位于飛行器橫向?qū)ΨQ面內(nèi)與Ox3軸垂直，Oz3軸垂直于其它兩軸并構(gòu)成右手直角坐標(biāo)系.

為了簡化推導(dǎo)，引入下列假設(shè).

假設(shè)一：飛行器攜帶軸對稱貯箱，且貯箱幾何中心放置于飛行器質(zhì)心處，此時對于軸對稱貯箱，其前兩階等效擺面法線分別與飛行器機體系x、y共線；

假設(shè)二：飛行器推力沿機體系x軸方向，由于飛行器本體姿態(tài)角變化很小，可以認(rèn)為三個姿態(tài)通道是解耦的.

假設(shè)三：把速度看做時間的已知函數(shù)V(t)=V0(t)，且小擾動、未擾動運動的側(cè)向參數(shù)及縱向角速度足夠小.

將飛行器機體系原點放置于系統(tǒng)質(zhì)心，將平動方程投影到彈道坐標(biāo)系，轉(zhuǎn)動方程投影到飛行器機體系、液體晃動方程投影到單擺本體系，并引入相關(guān)力與力系數(shù)，得到系統(tǒng)平動方程、轉(zhuǎn)動方程、第i階液體晃動方程增量方程如式（4）～（6）所示：

其中，θ、ψV為彈道傾角和彈道偏角；δy、δz為升降舵偏角和方向舵偏角；α和β為飛行迎角和側(cè)滑角；γ、φ、?為滾轉(zhuǎn)角、偏航角、俯仰角；P、X、Y、Z為發(fā)動機推力、氣動力在速度坐標(biāo)系下分量；G、V為系統(tǒng)重力和飛行器速度；ωxωyωz為飛行器角速度在機體坐標(biāo)系分量；JxJyJz為系統(tǒng)繞飛行器機體系軸轉(zhuǎn)動慣量；HixHiyHiz為第i階晃動擺相關(guān)項Jrsi在飛行器機體系下投影；PV為使飛行器速度V產(chǎn)生單位變化時所需推力增量，其余上標(biāo)含義類似.

2 飛行器縱向傳遞函數(shù)推導(dǎo)

3 升降舵到俯仰角速度控制分析

3.1 開環(huán)性能分析

對于有翼式飛行器，若不計重力影響，即a33=0，同時略去下洗動力系數(shù)a'24，則不考慮液體晃動時，由式（24）略去液體晃動參數(shù)an可得飛行器升降舵到俯仰角開環(huán)傳遞函數(shù)，并表示為典型基本環(huán)節(jié)傳遞函數(shù)：

3.2 剛體飛行器算例

預(yù)設(shè)飛行器質(zhì)量為1000kg，繞機體軸y軸轉(zhuǎn)動慣量為1000kg·m2.未擾時平飛運動速度V=100m/s；飛行器相關(guān)動力系數(shù)取值為：a22=-0.071，a24=-3.84，a25=-8，a35=0.005，a33=0.帶入上述系數(shù)，未受液體晃動擾動前飛行器升降舵到俯仰角速度開環(huán)傳函如下所示

其傳遞系數(shù)為KM=-0.1614：對于階躍升降舵輸入，響應(yīng)穩(wěn)定值為-0.1614；時間常數(shù)TM=0.5099，其固有頻率ωc=1/TM=1.961；相對阻尼系數(shù)ξM=0.0385；氣動力時間常數(shù)T1=12.1359.根據(jù)式（33）做出其零極點、根軌跡、伯德圖和階躍響應(yīng)曲線如圖2所示：如圖2所示，左上角為開環(huán)傳函零極點分布圖，可以看到開環(huán)極點實數(shù)小于零，所以開環(huán)系統(tǒng)穩(wěn)定.右上角為比例系數(shù)的根軌跡圖，可以看到負(fù)反饋閉環(huán)極點均在虛軸左側(cè)，故而其閉環(huán)穩(wěn)定.左下角為開環(huán)傳函的伯德圖，其相位裕量為79.4deg，增益裕量為15.8dB.左下角為其開環(huán)階躍響應(yīng).可以看到開環(huán)系統(tǒng)穩(wěn)定，但其階躍響應(yīng)超調(diào)量大，調(diào)節(jié)時間長，穩(wěn)定誤差大.

圖2 開環(huán)響應(yīng)曲線Fig.2 Open loop response curve

3.3 液體晃動對系統(tǒng)階躍響應(yīng)的影響

考慮圖3所示橫放圓柱貯箱，其長度為L，半徑為R，選取L=1m，R=0.5m，充液比50%的貯箱.貯箱中心放置于飛行器質(zhì)心處，軸向和飛行器機體軸重合；液體為航空煤油，此時液體重量為306.3kg；分析其對系統(tǒng)縱向傳遞函數(shù)的影響.式（27）和（32）分別為剛體飛行器和充液飛行器升降舵到俯仰角速度開環(huán)傳遞函數(shù).對比兩式，未充液系統(tǒng)（原始模型）在增加貯箱液體后，首先需要考慮飛行器質(zhì)量和慣性張量變化對動力系數(shù)的影響，即貯箱液體的重量和轉(zhuǎn)動慣量對傳函的影響（質(zhì)量慣量修正模型）；然后再考慮晃動參數(shù)ah對其的影響，即液體晃動對傳函的影響（液體晃動模型），其示意圖如圖4所示.

圖3 50%充液比貯箱內(nèi)液體的有限元模型Fig.3 Finite element model of liquid in tank with 50% liquid filling ratio

圖4 不同模型示意圖Fig.4 Schematic diagram of different models

未充液飛行器模型如式（33）所示，考慮橫放圓柱貯箱質(zhì)量及慣量，對動力系數(shù)a22，a24，a25，a34，a35，a33進行修正，得到其修正模型開環(huán)傳遞函數(shù)為：

最后考慮液體晃動相關(guān)參數(shù)ah對其的影響，得到其液體晃動模型開環(huán)傳遞函數(shù)為：

做出上述三種模型開環(huán)傳函零極點的圖像如圖5所示，可看到所有極點實部均小于零，即開環(huán)系統(tǒng)穩(wěn)定；修正模型相比于原始模型而言，在附加液體質(zhì)量和轉(zhuǎn)動慣量的影響下，開環(huán)傳函零極點均右移，其阻尼和自然頻率均減小，阻尼頻率變化不大，所以其階躍響應(yīng)波動增大，而周期基本不變；液體晃動模型相比于修正模型而言，其多出了兩對偶極子，零點沒有出現(xiàn)移動，主極點左移趨勢明顯，晃動頻率基本不變，但阻尼卻大幅降低，其階躍響應(yīng)振蕩加劇，動態(tài)響應(yīng)將進一步惡化.在后續(xù)分析中，忽略偶極子的影響，將液體晃動模型簡化為二階晃動模態(tài)進行分析，如式（35）所示.

圖5 不同模型零極點分布圖Fig.5 Distribution of poles and zeros of different models

不同模型的飛行器傳遞系數(shù)、時間常數(shù)、相對阻尼系數(shù)、氣動力時間常數(shù)如圖6所示，可以看到修正模型跟原始模型之間變化的主要是傳遞系數(shù)降低，即飛行器操縱性降低；氣動力時間常數(shù)增大，即零點出現(xiàn)大幅偏移.而晃動模型相對于修正模型，相對阻尼系數(shù)出現(xiàn)很大的變化，即其階躍響應(yīng)超調(diào)量將會增加，飛行器操縱時過載將會增加.不同模型的階躍響應(yīng)曲線如圖7所示，液體晃動對飛行器開環(huán)響應(yīng)影響很大，其超調(diào)量和調(diào)節(jié)時間均大幅增加.同時引起飛行器操縱性降低，過載增加.

圖6 不同模型系數(shù)變化示意圖Fig.6 Variation diagram of different model coefficients

圖7 開環(huán)階躍響應(yīng)Fig.7 Open loop step response

3.4 不同充液比對傳遞函數(shù)的影響

選取長度L=1，半徑R=0.5m的橫放圓柱貯箱，貯箱充液比變化時，分析其開環(huán)響應(yīng).充液比取0%～90%時，修正模型和晃動模型飛行器傳遞系數(shù)、時間常數(shù)、相對阻尼系數(shù)、氣動力時間常數(shù)如圖8所示，其充液0%時為無充液飛行器模型相關(guān)參數(shù).隨著充液比的增加，因為液體附加質(zhì)量和慣性的影響，飛行器傳遞系數(shù)逐漸降低，操縱性降低；飛行器時間常數(shù)變化不大，即飛行器固有頻率基本不變；相對阻尼系數(shù)逐漸降低，液體晃動對相對阻尼系數(shù)影響最為劇烈，即液體晃動會導(dǎo)致階躍響應(yīng)超調(diào)量和調(diào)節(jié)時間增加，動態(tài)穩(wěn)定性降低；氣動力時間常數(shù)逐漸增加，主要是由附加液體質(zhì)量和慣量引起，即隨著充液比的增加，零點絕對值逐漸減小.

3.5 不同長度L對傳遞函數(shù)的影響

保持圓柱半徑R=0.5m，充液比為50%.L分別取1m，1.2m1.4m，1.6m，1.8m，分析其對系統(tǒng)控制的影響.加入橫放圓柱貯箱后，隨著貯箱長度的增加L的增加，相關(guān)參數(shù)的變化如圖9所示，表中給出了不同長度下的修正模型和晃動模型飛行器傳遞系數(shù)、時間常數(shù)、相對阻尼系數(shù)、氣動力時間常數(shù)的數(shù)值.可以看到隨著貯箱長度L的增加，液體附加質(zhì)量和慣量引起了傳遞系數(shù)降低，時間常數(shù)增加，相對阻尼比降低，氣動力時間常數(shù)增加；分別代表飛行器操縱性降低，固有頻率降低，相對穩(wěn)定性變差，零點絕對值降低.而晃動模型相對修正模型，仍然是出現(xiàn)了相對阻尼系數(shù)大幅降低，階躍響應(yīng)超調(diào)量增加，調(diào)節(jié)時間變長，相對穩(wěn)定性進一步變差.即對于階躍響應(yīng)而言，振蕩峰值將會隨著L的增加而加劇，阻尼振蕩頻率減小，振蕩周期將增加.

圖9 不同L下模型傳函系數(shù)變化圖Fig.9 Variation of model transfer coefficient with different L

3.6 不同半徑R對開環(huán)響應(yīng)的影響

保持橫放貯箱長度L=1m，充液比為50%，充液液體為航空煤油，在常重g=10的情況下分析.半徑R分別取0.1～0.5m，分析其對系統(tǒng)控制的影響.隨著貯箱半徑尺寸的變化，修正模型和液體晃動模型典型基本環(huán)節(jié)表示的傳函相關(guān)參數(shù)如圖10所示。如圖10所示.隨著R的增加，修正模型和晃動模型飛行器傳遞系數(shù)減小，操縱性降低，機動性能降低；飛行器時間常數(shù)輕微增加，飛行器固有頻率稍稍降低；氣動力時間常數(shù)增大，零點絕對值減??；相對阻尼系數(shù)主要還是由液體晃動引起，階躍響應(yīng)超調(diào)量和調(diào)節(jié)時間均大幅增加，相對穩(wěn)定性降低.

圖10 不同R下模型傳函系數(shù)變化圖Fig.10 Variation of model transfer coefficient with different R

4 小結(jié)

本文針對攜帶液體貯箱的飛行器，基于等效力學(xué)模型將穩(wěn)態(tài)工作下的液體小幅晃動問題轉(zhuǎn)化為多體系統(tǒng)動力學(xué)問題，通過虛功率原理推導(dǎo)了攜帶多貯箱的多階耦合充液系統(tǒng)動力學(xué)方程，實現(xiàn)了液體、固體的實時耦合；并利用線性小擾動理論給出了充液飛行器流固耦合的線性小擾動方程，給出了縱向模態(tài)升降舵偏角到飛行器姿態(tài)傳函.

最后針對攜帶橫放圓柱貯箱的飛行器，利用主導(dǎo)極點概念，分別從貯箱內(nèi)液體充液比和貯箱尺寸的角度分析了液體晃動對飛行器縱向典型傳遞函數(shù)系數(shù)的影響.液體晃動主要影響了典型傳函中的相對阻尼系數(shù)，舵偏角階躍響應(yīng)超調(diào)量增加，調(diào)節(jié)時間變長，相對穩(wěn)定性變差.