張國(guó)斌，張青斌，豐志偉，陳青全，楊濤

國(guó)防科技大學(xué) 空天科學(xué)學(xué)院，長(zhǎng)沙 410073

空中加油是指加油機(jī)和受油機(jī)在飛行過(guò)程中以管道形成對(duì)接，然后將燃油輸送到受油機(jī)的過(guò)程。通過(guò)空中加油操作，可以大幅提高受油機(jī)的飛行半徑和留空時(shí)間，從而提升受油機(jī)的使用效能[1-2]。空中加油主要有2種方法：軟管錐套式空中加油和伸縮管式空中加油，這2種方法也分別被稱作軟式空中加油和硬式空中加油。

如圖1[3]所示，在軟式空中加油過(guò)程中，加油機(jī)保持一定的速度水平飛行，從加油吊艙中釋放輸油軟管，軟管末端裝有傘狀錐套[1]，在氣動(dòng)力的作用下，錐套穩(wěn)定于加油機(jī)后側(cè)的某個(gè)區(qū)域。在受油機(jī)機(jī)頭側(cè)方位置安裝有受油管，駕駛員操縱受油機(jī)使受油管插入錐套完成對(duì)接。對(duì)接完成后，燃油閥門(mén)打開(kāi)，燃油通過(guò)加油機(jī)傳送至受油機(jī)。相比于硬式空中加油，軟式空中加油設(shè)備輕便小巧，輸油過(guò)程中允許加油機(jī)與受油機(jī)間存在一定的相對(duì)運(yùn)動(dòng)，安全性得到了一定的提高。然而，由于完成對(duì)接后的相對(duì)運(yùn)動(dòng)控制難度大、對(duì)接點(diǎn)結(jié)構(gòu)薄弱，軟式空中加油目前仍是一項(xiàng)危險(xiǎn)性高、成功率低的空中作業(yè)項(xiàng)目，對(duì)飛行員的技術(shù)水平和飛機(jī)性能均有較高的要求。

圖1 軟式空中加油[3]Fig.1 Diagram of hose-drogue aerial refueling[3]

在軟式空中加油過(guò)程中，位于受油機(jī)側(cè)方的受油管需要保持錐套相對(duì)受油機(jī)靜止，往往要承受較大的約束力，而其結(jié)構(gòu)尺寸又小，在輸油過(guò)程中容易斷裂或從錐套中脫離。長(zhǎng)期的空中加油實(shí)踐和仿真均表明：加油機(jī)和受油機(jī)間的相對(duì)運(yùn)動(dòng)、氣流與輸油軟管間的流固耦合作用以及其他干擾因素會(huì)誘發(fā)輸油軟管產(chǎn)生甩鞭現(xiàn)象(Hose Whipping Phenomenon，HWP)[4]。甩鞭現(xiàn)象表現(xiàn)為輸油軟管的劇烈甩動(dòng)，同時(shí)受油管和錐套的連接處還會(huì)產(chǎn)生較大的瞬時(shí)載荷。由于甩鞭現(xiàn)象的出現(xiàn)，空中加油過(guò)程中受油管極有可能因?yàn)槭艿捷^大載荷而出現(xiàn)斷裂，進(jìn)而導(dǎo)致加油任務(wù)失敗，斷裂的油管甚至還有可能飛入受油機(jī)進(jìn)氣道釀成慘重的空中安全事故。因此，建立軟式空中加油對(duì)接狀態(tài)下的動(dòng)力學(xué)模型，并采用數(shù)值仿真手段獲取受油管受到的對(duì)接約束力具有重要的工程應(yīng)用價(jià)值。

對(duì)輸油軟管的建模多采用離散的質(zhì)點(diǎn)模型或基于球鉸連接的剛性桿模型。胡孟權(quán)等[5]將輸油軟管質(zhì)點(diǎn)離散化，研究了大氣紊流對(duì)油管錐套組合體運(yùn)動(dòng)的影響。Ro和Kamman[6]將輸油軟管離散為若干采用球鉸連接的剛性桿，研究了加油系統(tǒng)的動(dòng)力學(xué)特性。王海濤等[7]提出了一種長(zhǎng)度可變的多級(jí)串聯(lián)“球-桿”動(dòng)力學(xué)模型，基于此模型研究了空中加油過(guò)程中的控制策略。在串聯(lián)“球-桿”動(dòng)力學(xué)模型的基礎(chǔ)上，王海濤等[8]通過(guò)在相鄰兩剛體間引入彎曲恢復(fù)力，建立了軟管錐套組合體的動(dòng)力學(xué)模型，依靠該模型對(duì)甩鞭現(xiàn)象進(jìn)行了初步分析。采用這類集中參數(shù)模型可以較為便捷地建立起輸油軟管的動(dòng)力學(xué)模型，并能較好地反映軟管在加油過(guò)程中的動(dòng)力學(xué)特性，但是難以描述油管上連續(xù)分布的應(yīng)力和應(yīng)變，與實(shí)際情況的差異性較大。在輸油軟管的振動(dòng)控制研究中，Liu等[9-10]利用偏微分方程來(lái)描述加油過(guò)程中輸油軟管的動(dòng)力學(xué)特性。這種方法的控制方程在小變形假設(shè)下推導(dǎo)，僅考慮輸油軟管的橫向振動(dòng)，針對(duì)輸油軟管的振動(dòng)控制問(wèn)題取得了較好的效果。采用非線性有限元建立輸油軟管動(dòng)力學(xué)模型也是目前常用的一種建模方法。盡管理論推導(dǎo)較為復(fù)雜，數(shù)值仿真的計(jì)算量也稍大，但是非線性有限元基于連續(xù)介質(zhì)力學(xué)對(duì)輸油軟管進(jìn)行描述，理論基礎(chǔ)更加堅(jiān)實(shí)，能夠更加真實(shí)地反映輸油軟管本身的物理特性。Zhu和Meguid[11]開(kāi)發(fā)了一種三節(jié)點(diǎn)梁?jiǎn)卧糜谲浭娇罩屑佑偷膭?dòng)力學(xué)仿真，并通過(guò)試驗(yàn)驗(yàn)證了模型的有效性。劉釩[12]基于絕對(duì)節(jié)點(diǎn)坐標(biāo)法(Absolute Nodal Coordinate Formulation，ANCF)[13-14]建立了輸油軟管動(dòng)力學(xué)模型，對(duì)軟管錐套組合體在給定流場(chǎng)的穩(wěn)定拖曳姿態(tài)和釋放過(guò)程進(jìn)行了研究。在輸油過(guò)程中，燃油相對(duì)油管的流動(dòng)會(huì)形成較大的慣性力，但目前考慮燃油流動(dòng)的動(dòng)力學(xué)模型還較少。通過(guò)采用ALE(Arbitrary Lagrange-Euler)方式描述管內(nèi)流體單元，并結(jié)合ANCF方法，洪迪峰[15]建立了一種可描述結(jié)構(gòu)大位移和大轉(zhuǎn)動(dòng)的輸流管道的動(dòng)力學(xué)模型，該模型與ANCF結(jié)構(gòu)單元的兼容性極好，易于進(jìn)行系統(tǒng)動(dòng)力學(xué)的組裝?？罩屑佑瓦^(guò)程的動(dòng)力學(xué)建模涉及的學(xué)科領(lǐng)域廣、建模難度大，傳統(tǒng)的集中參數(shù)模型難以適應(yīng)幾何非線性條下柔性體高精度建模的要求。非線性有限元技術(shù)在近幾十年來(lái)實(shí)現(xiàn)了飛速的發(fā)展，結(jié)合計(jì)算多體動(dòng)力學(xué)建模技術(shù)可建立剛-柔-液耦合的高保真空中加油動(dòng)力學(xué)模型。

飛行速度、輸油軟管長(zhǎng)度、輸油軟管力學(xué)特性參數(shù)和其他不確定因素均可對(duì)空中加油過(guò)程中受油管和錐套間的約束力造成影響。這類不確定因素在某一區(qū)間內(nèi)服從一定的概率分布，如要獲取所有情況下對(duì)接約束力大小，需要仿真數(shù)目極大的工況，從而消耗大量的計(jì)算資源和計(jì)算時(shí)間。多項(xiàng)式混沌展開(kāi)是一種基于譜分析的方法，可以在輸入?yún)?shù)滿足一定概率分布的情況下快速給出輸出參數(shù)的期望值、方差和概率分布，而不需要大量的抽樣計(jì)算，能夠極大地提高分析計(jì)算的效率[16]。

本文首先基于絕對(duì)節(jié)點(diǎn)坐標(biāo)法建立輸油軟管的非線性有限元模型，利用ALE描述方法在輸油軟管單元控制方程上耦合燃油的動(dòng)力學(xué)方程，結(jié)合多體系統(tǒng)建模方法建立空中加油系統(tǒng)剛-柔-液耦合的多體動(dòng)力學(xué)模型，給出受油管指定截面所受約束力的計(jì)算方法；其次，引入基于多項(xiàng)式混沌展開(kāi)的不確定性分析方法；再次，通過(guò)動(dòng)力學(xué)仿真說(shuō)明甩鞭現(xiàn)象對(duì)于對(duì)接約束力的影響，對(duì)比分析出現(xiàn)甩鞭現(xiàn)象時(shí)輸油軟管上的受力分布特點(diǎn)；最后，在不確定性條件下，基于多項(xiàng)式混沌方法給出在不同巡航高度下，加油機(jī)與受油機(jī)有無(wú)相對(duì)運(yùn)動(dòng)時(shí)，對(duì)接約束力的期望范圍。

1 動(dòng)力學(xué)建模

1.1 坐標(biāo)系定義與基本假設(shè)

如圖2所示，在軟式空中加油過(guò)程中，加油機(jī)始終保持勻速平飛狀態(tài)，因此，將固定于加油機(jī)上的笛卡爾坐標(biāo)系oxyz認(rèn)定為慣性系，其中o點(diǎn)為輸油軟管在加油機(jī)上的釋放點(diǎn)，ox沿加油機(jī)水平飛行方向，oy軸指向重力反方向，oz軸按右手法則確定。在空中加油多體系統(tǒng)中，設(shè)輸油軟管oA采用球鉸分別與加油機(jī)和錐套連接于o點(diǎn)和A點(diǎn)。將錐套視為剛體，質(zhì)心為o1，體坐標(biāo)系為o1x1y1z1，在對(duì)接的初始狀態(tài)，o1x1y1z1與oxyz各坐標(biāo)軸平行。彎曲和拉壓是輸油軟管在加油過(guò)程中的主要變形模式，采用能夠描述彎曲和拉壓變形的絕對(duì)節(jié)點(diǎn)坐標(biāo)法梯度縮減梁?jiǎn)卧?Gradient Deficient ANCF Beam)對(duì)輸油軟管進(jìn)行建模[17]，這種單元采用了歐拉-伯努利梁假設(shè)(梁截面始終垂直于中心軸線且保持為平面)。當(dāng)油管中包含燃油時(shí)，根據(jù)彈性力學(xué)中對(duì)圓管受均布?jí)毫?wèn)題的推導(dǎo)[18]，油壓所產(chǎn)生的應(yīng)力分量平行于梁截面，本文認(rèn)為油壓不對(duì)梁?jiǎn)卧膹澢屠瓑鹤冃萎a(chǎn)生影響。在輸油過(guò)程中，將管內(nèi)燃油視為一維不可壓縮的理想流體，并且流速處處相等。

圖2 空中加油系統(tǒng)坐標(biāo)系定義Fig.2 Definition of coordinate frame of air refueling system

1.2 基于絕對(duì)節(jié)點(diǎn)坐標(biāo)法的輸油軟管建模

絕對(duì)節(jié)點(diǎn)坐標(biāo)法是由美國(guó)學(xué)者Shabana等[13]提出的一種柔性多體系統(tǒng)建模方法。絕對(duì)節(jié)點(diǎn)坐標(biāo)法采用微分幾何的方法描述連續(xù)體變形，引入的人為假設(shè)少，相比于傳統(tǒng)方法具有更高的插值精度。此外，絕對(duì)節(jié)點(diǎn)坐標(biāo)法具有質(zhì)量矩陣為常數(shù)、數(shù)值穩(wěn)定性好、易于與其他多體系統(tǒng)組合等特點(diǎn)，目前已經(jīng)成為多柔體系統(tǒng)最為有效的建模方法之一[19-21]。

1.2.1 軟管結(jié)構(gòu)動(dòng)力學(xué)模型

采用絕對(duì)節(jié)點(diǎn)坐標(biāo)法梯度縮減梁?jiǎn)卧⑤斢蛙浌芙Y(jié)構(gòu)的非線性有限元模型，該單元的參考構(gòu)型坐標(biāo)空間是一維的：只采用沿梁中心軸線的弧長(zhǎng)s作為標(biāo)定材料點(diǎn)的物質(zhì)坐標(biāo)，單元原長(zhǎng)為L(zhǎng)，如圖3(a)所示。單元包含2個(gè)節(jié)點(diǎn)，節(jié)點(diǎn)坐標(biāo)包含一個(gè)位置向量ri和梯度向量ri,s(i=1,2)，如圖3(b)所示。絕對(duì)節(jié)點(diǎn)坐標(biāo)法中所有矢量均在慣性系下表示。

圖3 絕對(duì)節(jié)點(diǎn)坐標(biāo)法縮減梁?jiǎn)卧狥ig.3 Gradient deficient ANCF beam

基于有限元技術(shù)，單元內(nèi)任意物質(zhì)點(diǎn)s在當(dāng)前構(gòu)形的空間位置向量用形函數(shù)矩陣和節(jié)點(diǎn)坐標(biāo)向量表示形式為

r(s)=N(s)qe

(1)

式中：qe為單元節(jié)點(diǎn)坐標(biāo)向量；N(s)為形函數(shù)矩陣；qe的具體形式為

(2)

N(s)具體為

(3)

其中:I為3×3單位陣;4個(gè)形函數(shù)分別為

(4)

式中：ξ=s/L。

采用虛功原理或第2類Lagrange方程均可以推導(dǎo)出單元的動(dòng)力學(xué)方程：

(5)

(6)

式(6)在參考構(gòu)型中積分，其中：ρh和Ah分別為輸油軟管在變形前的密度和截面積。

單元的廣義內(nèi)力由單元變形能V對(duì)單元節(jié)點(diǎn)坐標(biāo)求導(dǎo)得到，即

(7)

由于采用了歐拉-伯努利梁假設(shè)，并考慮到輸油軟管截面形狀為圓環(huán)，單元的變形能為

(8)

式(8)同樣在參考構(gòu)型下積分，其中：E為輸油軟管材料的楊氏模量；Ih為輸油軟管截面慣性矩；ε為軸向應(yīng)變；κ為曲率。通過(guò)微分幾何可知

(9)

(10)

式中：rs表示r對(duì)s求一次偏導(dǎo)；rss表示r對(duì)s求2次偏導(dǎo)。由虛功原理，作用在材料點(diǎn)s上的外力F所做的虛功為

(11)

故力F所對(duì)應(yīng)的廣義力為

QF=N(s)TF

(12)

由重力g所產(chǎn)生的廣義力為

(13)

1.2.2 管內(nèi)燃油動(dòng)力學(xué)模型

加油過(guò)程中，輸油軟管中的燃油可認(rèn)為是一維不可壓縮的理想流體，將燃油相對(duì)油管的流速設(shè)為U。如圖4所示，在輸油軟管單元兩端的節(jié)點(diǎn)上，燃油的物質(zhì)坐標(biāo)分別為

圖4 管內(nèi)燃油模型Fig.4 Model for oil inside the hose

(14)

式中：U為燃油相對(duì)油管的流速；t為流動(dòng)時(shí)間。

軟管內(nèi)燃油的運(yùn)動(dòng)學(xué)描述同樣采用絕對(duì)節(jié)點(diǎn)坐標(biāo)法實(shí)現(xiàn)，并與軟管結(jié)構(gòu)單元共用形函數(shù)和節(jié)點(diǎn)坐標(biāo)。燃油物質(zhì)點(diǎn)的物質(zhì)坐標(biāo)p與軟管單元形函數(shù)參數(shù)s的關(guān)系為

s=p+Ut

(15)

式(15)表征了燃油在輸油軟管內(nèi)的流動(dòng)。

軟管內(nèi)燃油物質(zhì)點(diǎn)的空間位置可由形函數(shù)表達(dá)為

r(p)=N(s)qe

(16)

對(duì)式(16)求物質(zhì)導(dǎo)數(shù)，得到燃油物質(zhì)點(diǎn)的速度和加速度分別為

(17)

其中：

(18)

不考慮燃油的黏性，基于虛功原理，燃油的慣性力所做虛功和重力所做的虛功之和為0，即

(19)

式中：ρo為燃油密度；Ao為油管空腔截面積。經(jīng)過(guò)代數(shù)運(yùn)算后可得

(20)

其中：

(21)

結(jié)合軟管結(jié)構(gòu)的動(dòng)力學(xué)方程，可以得到考慮燃油流動(dòng)的輸油軟管單元?jiǎng)恿W(xué)微分方程為

(22)

1.3 流場(chǎng)模型與氣動(dòng)力估算

1.3.1 流場(chǎng)模型

根據(jù)1.1節(jié)中的假設(shè)，由于采用固連于加油機(jī)上的慣性坐標(biāo)系oxyz來(lái)度量空中加油過(guò)程中的所有物理量，故加油機(jī)后方流場(chǎng)可近似為直勻流疊加加油機(jī)尾渦流場(chǎng)?，F(xiàn)有的尾渦流場(chǎng)模型包括Rankine vortex模型、Adapted vortex模型、Lamb-Oseen vortex模型、Smooth blending vortex profile模型和Hallock-Burnham vortex模型等[22-23]。在空中加油的仿真計(jì)算中，通常采用Hallock-Burnham vortex模型計(jì)算尾渦產(chǎn)生的速度場(chǎng)[24-25]。

圖5 尾渦速度場(chǎng)模型Fig.5 Model for wake vortex velocity field

(23)

(24)

取加油機(jī)重量為185 t，平飛速度為105 m/s，飛機(jī)翼展為49.5 m，基于Hallock-Burnham模型計(jì)算由尾渦誘導(dǎo)出的流場(chǎng)，在尾渦渦心所在的鉛垂面上，尾流場(chǎng)如圖6所示。相比于加油機(jī)平飛速度105 m/s，尾渦所誘導(dǎo)出的最大速度小于12 m/s，二者的合成速度為105.7 m/s，可見(jiàn)尾渦對(duì)相對(duì)流動(dòng)速度的影響較小。

圖6 尾渦中心所在截面的尾渦流場(chǎng)鉛錘切面Fig.6 Vertical slice of wake flow fieldat wake center

1.3.2 氣動(dòng)力計(jì)算

輸油軟管微段dl受到的氣動(dòng)力可以表示為如下形式：

(25)

式中：Cd為輸油軟管受到的壓差阻力系數(shù)；Cf為輸油軟管的摩擦阻力系數(shù)；ρa(bǔ)為氣流密度；d為輸油軟管外徑；vn為輸油軟管微段相對(duì)氣流的法向速度；vt為輸油軟管相對(duì)氣流的切向速度。采用試驗(yàn)公式[27-28]計(jì)算油管受到的摩擦力系數(shù)Cf和壓差力系數(shù)Cd：

(26)

Cd=

(27)

式中：Ref和Rep為雷諾數(shù)，其計(jì)算式分別為

其中：ν為空氣黏度系數(shù)。

結(jié)合式(12)，輸油軟管單元受到氣動(dòng)力所對(duì)應(yīng)的廣義外力為

(28)

式中：第1個(gè)積分式的積分域?yàn)楫?dāng)前構(gòu)型；l為當(dāng)前單元長(zhǎng)度；第2個(gè)積分式的積分域?yàn)閰⒖紭?gòu)型。在動(dòng)力學(xué)程序中，對(duì)式(28)的積分可采用Gauss方法進(jìn)行數(shù)值計(jì)算[29]。

為簡(jiǎn)化計(jì)算，假設(shè)錐套受到的氣動(dòng)力通過(guò)錐套質(zhì)心，計(jì)算公式為

(29)

式中：v為氣流相對(duì)錐套質(zhì)心的速度；Sd為錐套的特征阻力面積；Cd為錐套的氣動(dòng)阻力系數(shù)。

1.4 多體系統(tǒng)組裝與對(duì)接約束力求解

1.4.1 多體系統(tǒng)組裝

利用絕對(duì)節(jié)點(diǎn)坐標(biāo)法建立描述輸油軟管和燃油動(dòng)力學(xué)的單元，對(duì)輸油軟管劃分單元再進(jìn)行組裝后可得到輸油軟管部分的動(dòng)力學(xué)方程為

(30)

式中：Mh為輸油軟管及內(nèi)部燃油所對(duì)應(yīng)的廣義質(zhì)量陣；qh為單元節(jié)點(diǎn)坐標(biāo)對(duì)應(yīng)的列向量；Qe為油管彈性內(nèi)力所對(duì)應(yīng)的廣義力列向量；Qa為油管所受到的氣動(dòng)力列向量；QU為燃油流動(dòng)所形成的附加廣義慣性力；Qg為油管及內(nèi)部燃油所對(duì)應(yīng)的重力廣義力列向量。當(dāng)輸油軟管的劃分單元數(shù)為ne時(shí)，廣義質(zhì)量陣Mh所對(duì)應(yīng)的階數(shù)為6(ne+1)。

(31)

式中：md為錐套的質(zhì)量；Jx、Jy和Jz分別為錐套在本體系上主慣量。錐套的動(dòng)力學(xué)方程為

(32)

(33)

式中:λ0、λ1、λ2和λ3分別為姿態(tài)四元數(shù)ε對(duì)應(yīng)的4個(gè)分量。

Φ(q)=012×1

對(duì)約束方程關(guān)于時(shí)間求一次導(dǎo)數(shù)，并利用式(33) 將四元數(shù)的導(dǎo)數(shù)替換為錐套角速度可以得到

繼續(xù)關(guān)于時(shí)間求導(dǎo)可得

(34)

結(jié)合式(30)、式(32)和式(34)，可得空中加油多體系統(tǒng)動(dòng)力學(xué)微分方程為

(35)

式中：

(36)

式(35)為指標(biāo)1微分代數(shù)方程，可以通過(guò)常微分求解器做積分計(jì)算。

1.4.2 對(duì)接約束力求解

(37)

使廣義約束力QΦ的分量對(duì)應(yīng)于相應(yīng)的廣義坐標(biāo)，則可得到相應(yīng)的約束力分量。

在對(duì)受油管進(jìn)行強(qiáng)度分析時(shí)，需要知道對(duì)接約束力作用于受油管截面的拉壓力、剪切力及彎矩。由式(37)，可以推算得到受油管作用于錐套質(zhì)心的約束力向量FΦ(投影在全局坐標(biāo)系oxyz)及對(duì)錐套的約束力矩向量MΦ(投影在錐套本體系o1x1y1z1)。根據(jù)牛頓第三定律，可知錐套施加給受油管的力系可等效為作用位置在錐套質(zhì)心的集中力-FΦ以及力矩-MΦ。在對(duì)接過(guò)程中，受油管與錐套連接為一體，令受油管的本體系與錐套本體系方向一致，如圖7(a)所示，則受油管截面B-B受到的彎矩為

(38)

圖7 受油管受力狀態(tài)Fig.7 Force on oil tube

(39)

2 不確定性分析方法

空中加油過(guò)程對(duì)接約束力的求解需要涉及眾多的輸入變量，包括：飛行速度、輸油軟管長(zhǎng)度、飛行高度和相對(duì)運(yùn)動(dòng)速度等。如果采用經(jīng)典的Monte-Carlo方法，則需要規(guī)模巨大的仿真任務(wù)，耗費(fèi)大量的計(jì)算資源和時(shí)間。高效的不確定性分析方法能夠用較少仿真次數(shù)得到多個(gè)輸入量在一定概率分布下輸出量的統(tǒng)計(jì)特征，從而大大提高工作效率。多項(xiàng)式混沌方法是目前較為成熟的一種不確定性分析方法，在實(shí)際工程問(wèn)題中得到了廣泛的應(yīng)用[31-33]。

2.1 多項(xiàng)式混沌展開(kāi)理論

多項(xiàng)式混沌展開(kāi)的基本思想是用含獨(dú)立隨機(jī)變量的正交多項(xiàng)式之和來(lái)近似表示隨機(jī)過(guò)程[34]。即將系統(tǒng)輸出量y和不確定輸入向量x的關(guān)系表示為

(40)

式中：βi為常系數(shù)；φi(x)為多元多項(xiàng)式基函數(shù)；n為非常數(shù)多項(xiàng)式基函數(shù)的個(gè)數(shù)。

(41)

式中：〈·〉代表內(nèi)積運(yùn)算；Ω為積分域；γi=〈φi(x),φi(x)〉為標(biāo)準(zhǔn)化常量；δij為克羅內(nèi)克張量。當(dāng)給定隨機(jī)變量的概率密度函數(shù)ω(x)后，定義正交多項(xiàng)式基的第1項(xiàng)φ0(x)=1，再基于Gram-Schimidt正交化方法，利用式(41)所定義的內(nèi)積運(yùn)算，對(duì)線性無(wú)關(guān)的函數(shù)系

{φi|φi(x)=xi,i=0,1,2,…}

做正交化處理，即可得到在給定概率分布下的正交多項(xiàng)式基。照此方法可知：正態(tài)分布隨機(jī)變量對(duì)應(yīng)的正交多項(xiàng)式基為Hermite多項(xiàng)式；均勻分布隨機(jī)變量對(duì)應(yīng)的正交多項(xiàng)式基為L(zhǎng)egendre多項(xiàng)式；貝塔分布隨機(jī)變量對(duì)應(yīng)的正交多項(xiàng)式為Jacobi多項(xiàng)式。Legendre多項(xiàng)式的前4項(xiàng)如表1所示。

表1 前4項(xiàng)Legendre多項(xiàng)式Table 1 First four Legendre polynomials

若系統(tǒng)輸入含有多個(gè)獨(dú)立隨機(jī)變量，則可將單變量對(duì)應(yīng)的正交多項(xiàng)式按獨(dú)立組合的原則連乘再相加得到多維隨機(jī)變量的多項(xiàng)式混沌展開(kāi)[33]：

(42)

式中：q為獨(dú)立隨機(jī)變量個(gè)數(shù)；βi1,…,iq為多項(xiàng)式混沌展開(kāi)系數(shù)；nk(k=1,2,…,q)是變量xk的正交多項(xiàng)式最大階數(shù)，基函數(shù)可表示為

ψi1,…,iq(x)=φi1(x1)φi2(x2)…φiq(xq)

其中：φik(xk)是隨機(jī)變量xk的第ik次項(xiàng)。在不確定性分析問(wèn)題的研究中，一般對(duì)式(42)按照最大階次進(jìn)行截?cái)?。若選定的最大階次為p，則式(40) 中需要求的待定系數(shù)個(gè)數(shù)為

(43)

對(duì)系統(tǒng)輸出量完成多項(xiàng)式展開(kāi)后，由正交多項(xiàng)式的基本性質(zhì)，系統(tǒng)輸出量y的期望和方差分別為

E(y)=β0

(44)

(45)

2.2 多項(xiàng)式混沌展開(kāi)系數(shù)計(jì)算

多項(xiàng)式混沌展開(kāi)模型的系數(shù)確定是多項(xiàng)式混沌方法的核心步驟?，F(xiàn)有的系數(shù)求解方法可分為嵌入法(Intrusive Method)和非嵌入法(Non-intrusive Method)[35]。嵌入法在使用時(shí)需要修改動(dòng)力學(xué)模型，不適用于大型動(dòng)力學(xué)仿真程序，而非嵌入法則不需要修改動(dòng)力學(xué)模型，應(yīng)用范圍更廣。

本文討論非嵌入法中的隨機(jī)響應(yīng)面法[34]。對(duì)系統(tǒng)輸入隨機(jī)變量x構(gòu)成的多維空間，每次動(dòng)力學(xué)計(jì)算的輸入向量就對(duì)應(yīng)這個(gè)空間中的一個(gè)點(diǎn)，通常稱這個(gè)點(diǎn)為采用點(diǎn)或配置點(diǎn)。要求解多項(xiàng)式混沌展開(kāi)中的未知系數(shù)，需選取足夠多的配點(diǎn)。對(duì)q維p階多項(xiàng)混沌展開(kāi)，隨機(jī)響應(yīng)面法的取點(diǎn)原則是：對(duì)每個(gè)變量取p+1階多項(xiàng)式的根，然后將這些根進(jìn)行組合，最后得到(p+1)q個(gè)采樣點(diǎn)。將所有的采樣點(diǎn)輸入動(dòng)力學(xué)仿真程序計(jì)算將得到(p+1)q個(gè)系統(tǒng)輸出yi，參考式(40)，可以得到方程組：

(46)

記式(46)為Hβ=y，則可對(duì)多項(xiàng)混沌系數(shù)進(jìn)行最優(yōu)估計(jì)：

β=(HTH)-1HTy

(47)

3 不確定性條件下的仿真分析

3.1節(jié)首先對(duì)正常對(duì)接與出現(xiàn)甩鞭現(xiàn)象時(shí)的對(duì)接情況進(jìn)行對(duì)比分析。本文假設(shè)加油機(jī)飛行速度、輸油軟管長(zhǎng)度和對(duì)接相對(duì)運(yùn)動(dòng)速度大小都服從指定區(qū)間的均勻分布。3.2節(jié)分析在不同高度下，油管長(zhǎng)度和飛行速度為隨機(jī)變量時(shí)，受油管指定截面上的拉壓力、剪切力和彎矩的期望和標(biāo)準(zhǔn)差。3.3節(jié)分析不同高度下，受油機(jī)相對(duì)運(yùn)動(dòng)速度大小和輸油軟管長(zhǎng)度為隨機(jī)變量時(shí)，受油管指定截面的拉壓力、剪切力和彎矩的期望和標(biāo)準(zhǔn)差。

3.1 甩鞭現(xiàn)象對(duì)應(yīng)的約束力變化

本文針對(duì)某型加油機(jī)的空中加油過(guò)程進(jìn)行研究。在500 m的飛行高度下，環(huán)境參數(shù)為：空氣密度為1.17 kg/m3、空氣黏度系數(shù)為1.79×10-5、重力加速度為9.8 m/s2。取加油機(jī)飛行速度為155 m/s，油管長(zhǎng)度為24 m。其他參數(shù)如表2所示。油管劃分單元數(shù)為10，通過(guò)仿真分別得到了在正常對(duì)接和受油機(jī)相對(duì)速度為(1, 0, 0)m/s(投影在oxyz系)時(shí)受油管某一截面約束力隨時(shí)間變化曲線，如圖8所示。從圖8中可見(jiàn)，受油管截面的壓力、剪力和彎矩都出現(xiàn)了明顯的增大，特別是剪力和彎矩的增幅最為明顯，出現(xiàn)甩鞭現(xiàn)象后最大彎矩增加至正常對(duì)接的6.1倍；從曲線的波動(dòng)變化情況可以推測(cè)此時(shí)輸油軟管上出現(xiàn)了較高頻率的振動(dòng)。甩鞭現(xiàn)象對(duì)應(yīng)的錐套油管組合體的運(yùn)動(dòng)過(guò)程如圖9所示，圖9中預(yù)測(cè)的輸油軟管位形與實(shí)際甩鞭現(xiàn)象中拍攝的照片吻合度較好[8]。

圖8 受油管某考察截面的受力Fig.8 Force acting on the section of probe under investigation

圖9 甩鞭現(xiàn)象對(duì)應(yīng)的軟管位形Fig.9 Hose configuration corresponding to HWP

表2 系統(tǒng)參數(shù)Table 2 System parameters

輸油軟管出現(xiàn)甩鞭現(xiàn)象和正常對(duì)接時(shí)的彎矩分布云圖分別如圖10(a)和圖10(b)所示。從圖中可見(jiàn)，輸油軟管最大截面彎矩均出現(xiàn)在靠近錐套的部分，相比于正常對(duì)接狀態(tài)，出現(xiàn)甩鞭現(xiàn)象時(shí)輸油軟管彎矩明顯增大。輸油軟管截面的壓力如圖11所示(注意圖11中壓為正、拉為負(fù))，對(duì)比圖11(a) 和圖11(b)可知：發(fā)生甩鞭現(xiàn)象時(shí)，輸油軟管截面主要受壓，在輸油軟管曲率較大的部位會(huì)出現(xiàn)較大的拉力；而在正常對(duì)接狀態(tài)時(shí)，輸油軟管截面受拉，最大拉力出現(xiàn)在輸油軟管與加油機(jī)連接的部位；出現(xiàn)甩鞭現(xiàn)象時(shí)輸油軟管截面受到的拉壓力會(huì)比正常對(duì)接時(shí)的拉壓力偏大。

圖10 油管截面彎矩Fig.10 Bending moment on the hose section

圖11 油管截面壓力Fig.11 Pressure on the hose section

綜合3.1節(jié)的討論結(jié)果，甩鞭現(xiàn)象會(huì)導(dǎo)致作用在受油管截面的最大拉壓力、最大剪切力和最大彎矩出現(xiàn)顯著增大，約束力的變化頻率和變化幅度也會(huì)顯著的提高。此外，甩鞭現(xiàn)象還會(huì)使得輸油軟管受力狀態(tài)發(fā)生較大的改變，在輸油管軟管靠近錐套的部分受力增加尤為明顯。因此，需要重點(diǎn)考察加油機(jī)與受油機(jī)間相對(duì)運(yùn)動(dòng)而誘發(fā)的甩鞭現(xiàn)象。

3.2 輸油軟管長(zhǎng)度和飛行速度不確定

3.2節(jié)分析在給定高度下輸油軟管長(zhǎng)度和飛行速度為不確定量時(shí)，受油管指定截面的受力期望曲線和最大概率分布范圍。輸油軟管長(zhǎng)度為在區(qū)間[24, 27] m中均勻分布的隨機(jī)變量，加油機(jī)飛行速度同樣服從均勻分布，不同高度下，加油機(jī)飛行速度區(qū)間如表3所示。結(jié)合第2節(jié)介紹的基于多項(xiàng)式混沌的不確定性分析理論，取多項(xiàng)式混沌的最大展開(kāi)階數(shù)為兩階，利用2.2節(jié)中的隨機(jī)響應(yīng)面法確定多項(xiàng)式混沌的展開(kāi)系數(shù)。系統(tǒng)中的其他參數(shù)如表2所示。分析結(jié)果如圖12 所示。

表3 不同高度下加油機(jī)平飛速度Table 3 Horizontal speed of the tanker at different altitudes

從圖12中的不確定性分析結(jié)果來(lái)看，不同高度下，受油管截面所受到的壓力、剪力和彎矩隨時(shí)間變化走勢(shì)基本相同，而且大小相差不大。這是因?yàn)樵诓煌叨认拢淖兊牧繛榭諝饷芏?、重力加速度及飛行速度，盡管在高空的飛行速度更高，但此時(shí)空氣密度和重力加速度都有減小，因此在不同高度下，油管截面的受力相差不大。油管截面受到的壓力上確界最大概率為6 880 N，剪力上確界最大概率為1 325 N，彎矩上確界最大概率為338 N·m。

圖12 受油管截面受力的不確定分析(無(wú)甩鞭現(xiàn)象)Fig.12 Uncertainty analysis of force on the probe section(without HWP)

3.3 輸油軟管長(zhǎng)度和對(duì)接速度不確定

從3.1節(jié)中的分析得知，甩鞭現(xiàn)象會(huì)使得受油管截面的受力大大增加。甩鞭現(xiàn)象的產(chǎn)生主要是因?yàn)槭苡蜋C(jī)相對(duì)加油機(jī)的運(yùn)動(dòng)導(dǎo)致了輸油軟管出現(xiàn)屈曲變形和高頻振動(dòng)。本節(jié)中以受油機(jī)相對(duì)加油的運(yùn)動(dòng)速度大小和油管長(zhǎng)度為系統(tǒng)輸入隨機(jī)變量，并假定相對(duì)運(yùn)動(dòng)速度在區(qū)間[0, 3] m/s內(nèi)均勻分布，油管長(zhǎng)度在區(qū)間[24, 27] m內(nèi)均勻分布。在500 m，6 000 m和9 000 m高度下按表3中的最大飛行速度對(duì)受油管截面受力進(jìn)行不確定性分析。系統(tǒng)其他參數(shù)按照表2輸入，油管截面的受到的約束力如圖13所示。對(duì)比圖13與圖12，出現(xiàn)甩鞭現(xiàn)象后，油管截面的受力大幅增加。截面受到的壓力上確界最大概率為10 770 N，剪力上確界最大概率為4 479 N，彎矩上確界最大概率為1 474 N·m。相比于不存在相對(duì)運(yùn)動(dòng)的對(duì)接過(guò)程，受油管在特定截面處的最大壓力、剪力和彎矩的期望值分別提高了50%、272%和772%，可見(jiàn)，加油機(jī)與受油機(jī)間的相對(duì)運(yùn)動(dòng)是相關(guān)部件強(qiáng)度設(shè)計(jì)需要考慮的重要因素。

圖13 油管截面受力的不確定分析(存在甩鞭現(xiàn)象)Fig.13 Uncertainty analysis of force on the probe section (with HWP)

4 結(jié) 論

1) 采用絕點(diǎn)節(jié)點(diǎn)坐標(biāo)法和ALE描述方式建立了輸油軟管的動(dòng)力學(xué)模型，采用剛體假設(shè)建立錐套模型，疊加直勻流和Hallock-Burnham尾渦模型建立簡(jiǎn)化流場(chǎng)模型，采用試驗(yàn)公式擬合油管受到的氣動(dòng)力，最終建立了空中加油系統(tǒng)剛-柔-液耦合的多體動(dòng)力學(xué)模型，該模型可用于計(jì)算空中加油對(duì)接約束力，并能夠預(yù)測(cè)甩鞭現(xiàn)象。相比于集中參數(shù)模型，本文采用非線性有限元技術(shù)建立的輸油軟管模型能更加真實(shí)地反映柔性油管的動(dòng)力學(xué)特性及燃油流動(dòng)所對(duì)系統(tǒng)動(dòng)力學(xué)的影響，并能得到輸油軟管上連續(xù)分布的受力狀態(tài)。

2) 仿真發(fā)現(xiàn)：相比于正常對(duì)接過(guò)程，出現(xiàn)甩鞭現(xiàn)象時(shí)，受油管截面受到的壓力、剪力和彎矩都將大幅增大；輸油軟管的受力狀態(tài)將會(huì)發(fā)生較大的改變，軸向變形由受拉轉(zhuǎn)變?yōu)橹饕軌?，截面彎矩將?huì)大幅增大，同時(shí)還會(huì)出現(xiàn)較高頻率的振動(dòng)。

3) 多項(xiàng)式混沌方法可用于輸入?yún)?shù)不確定情況下對(duì)接約束力的分析計(jì)算，能夠以較少的計(jì)算量給出對(duì)接約束力的期望值和標(biāo)準(zhǔn)差，估計(jì)出最大概率意義下的約束力的上界。