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

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

引言

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

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

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

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

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

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

上述中，Rc為飛行器在慣性系下位移，ωc為飛行器繞機(jī)體系轉(zhuǎn)動(dòng)角速度，q為單擺相對(duì)飛行器機(jī)體系擺角；m、S、J為系統(tǒng)質(zhì)量、系統(tǒng)相對(duì)飛行器機(jī)體系原點(diǎn)靜矩、轉(zhuǎn)動(dòng)慣量矩陣，Spi、Jpi為第i階單擺相對(duì)單擺系原點(diǎn)靜矩、轉(zhuǎn)動(dòng)慣量；g為重力加速度；ωpi為第i階單擺相對(duì)飛行器機(jī)體系角速度；F、T為飛行器所受氣動(dòng)力和氣動(dòng)力矩；Qt、Qr分別為液體晃動(dòng)導(dǎo)致的平動(dòng)耦合力、轉(zhuǎn)動(dòng)耦合力；rci、為第i階單擺相對(duì)單擺懸掛點(diǎn)矢量；Hpi為第i階單擺擺動(dòng)面法線方向；SpHi、JpHi、Jrsi為與晃動(dòng)擺相關(guān)項(xiàng)具體表達(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)系原點(diǎn)取在飛行器發(fā)射點(diǎn)上，x軸沿彈道面與地面交線指向目標(biāo)點(diǎn)處，z軸沿鉛垂面向下，y軸與其它兩軸垂直并構(gòu)成右手系.

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

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

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

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

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

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

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

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

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

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

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

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

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

3.2 剛體飛行器算例

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

其傳遞系數(shù)為KM=-0.1614：對(duì)于階躍升降舵輸入，響應(yīng)穩(wěn)定值為-0.1614；時(shí)間常數(shù)TM=0.5099，其固有頻率ωc=1/TM=1.961；相對(duì)阻尼系數(shù)ξM=0.0385；氣動(dòng)力時(shí)間常數(shù)T1=12.1359.根據(jù)式（33）做出其零極點(diǎn)、根軌跡、伯德圖和階躍響應(yīng)曲線如圖2所示：如圖2所示，左上角為開環(huán)傳函零極點(diǎn)分布圖，可以看到開環(huán)極點(diǎn)實(shí)數(shù)小于零，所以開環(huán)系統(tǒng)穩(wěn)定.右上角為比例系數(shù)的根軌跡圖，可以看到負(fù)反饋閉環(huán)極點(diǎ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é)時(shí)間長(zhǎng)，穩(wěn)定誤差大.

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

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

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

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

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

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

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

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

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

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

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

3.5 不同長(zhǎng)度L對(duì)傳遞函數(shù)的影響

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

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

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

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

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

4 小結(jié)

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

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