劉文一，李玉龍，吳訓(xùn)濤

（中國(guó)人民解放軍91550部隊(duì)41分隊(duì)，大連，116023）

0 引 言

飛行器水下垂直向上運(yùn)動(dòng)時(shí)，由于空泡的出現(xiàn)，會(huì)經(jīng)受復(fù)雜的水下載荷，在飛行器表面布置壓強(qiáng)傳感器來測(cè)量表面壓強(qiáng)，由于受傳感器安裝位置和采樣頻率的限制，測(cè)得的壓強(qiáng)數(shù)據(jù)在時(shí)間和空間上是離散的。為了對(duì)飛行器進(jìn)行動(dòng)力學(xué)響應(yīng)分析，必須將離散的數(shù)據(jù)連續(xù)化，才能獲得飛行器的載荷分布情況[1]。

獲得飛行器表面壓強(qiáng)場(chǎng)后，將表面壓強(qiáng)場(chǎng)加載在飛行器表面上進(jìn)動(dòng)力學(xué)響應(yīng)分析，以獲得飛行器對(duì)表面壓強(qiáng)的動(dòng)響應(yīng)，從而辨識(shí)得到飛行器的內(nèi)力，即飛行器各截面隨時(shí)間和空間變化的軸向力、剪切力和彎矩[2]。對(duì)于長(zhǎng)徑比較大的飛行器，剪切力和彎矩是其水下載荷的主要部分，這兩個(gè)值直接影響飛行能否正常出水和出水姿態(tài)，因此辨識(shí)得到的飛行器水下載荷具有重要的工程應(yīng)用價(jià)值[3]。

1 表面壓強(qiáng)

對(duì)于飛行器水下載荷分析而言，主要關(guān)心飛行器截面分布法向力。實(shí)驗(yàn)時(shí)在飛行器表面布置若干個(gè)壓強(qiáng)測(cè)點(diǎn)，然后對(duì)測(cè)得的表面壓強(qiáng)進(jìn)行插值求得飛行器連續(xù)表面壓強(qiáng)場(chǎng)。飛行器在出水過程中法向力分布會(huì)出現(xiàn)一個(gè)峰值，而且峰值并不在一個(gè)固定的部位，瞬時(shí)出現(xiàn)法向力峰值的主要原因是由于空泡產(chǎn)生了傾斜。對(duì)于同一截面，當(dāng)迎水面處于空泡末端高壓作用時(shí)，與其對(duì)應(yīng)的背水面已經(jīng)進(jìn)入空泡低壓區(qū)，形成較大的迎、背水壓差，因此產(chǎn)生了較大的瞬時(shí)截面分布法向力峰值。同時(shí)，在出水過程中，隨著速度增加，空泡數(shù)隨之降低，空泡末端渦旋回射產(chǎn)生的壓強(qiáng)峰逐漸后移，由此引起的法向力分布峰值也逐漸后移[4]。

飛行器在水下任一點(diǎn)的表面壓強(qiáng)可表示為

式中 CP為壓強(qiáng)系數(shù)；ρ為水的密度，kg/m2；V∞為飛行器運(yùn)動(dòng)速度，m/s；P0為環(huán)境壓強(qiáng)，Pa。

由于飛行器在水下會(huì)產(chǎn)生空化現(xiàn)象，因此載荷計(jì)算關(guān)心的是CP為負(fù)值的部分，引入符號(hào)ζ：

2 連續(xù)壓強(qiáng)場(chǎng)

飛行器體表面壓強(qiáng)測(cè)量數(shù)據(jù)測(cè)得了每個(gè)測(cè)點(diǎn)上壓強(qiáng)隨時(shí)間的變化歷程，但對(duì)于兩個(gè)測(cè)點(diǎn)之間的壓強(qiáng)分布，則必須對(duì)壓強(qiáng)測(cè)量數(shù)據(jù)進(jìn)行連續(xù)化處理，利用波形插值法可以將離散的壓強(qiáng)數(shù)據(jù)插值得到隨時(shí)間變化的壓強(qiáng)場(chǎng)，這是因?yàn)樵诔鏊^程中，空泡沿飛行器表面是按照一定速度推進(jìn)潰滅的[5]。

進(jìn)行波形插值時(shí)，為確定某點(diǎn)的壓強(qiáng)載荷，首先根據(jù)其位置確定出水時(shí)刻，移動(dòng)時(shí)間軸，將出水時(shí)刻的主波峰對(duì)齊，然后進(jìn)行線性插值得到峰值附近的壓強(qiáng)曲線；最后還原時(shí)間軸，得到該點(diǎn)的壓強(qiáng)載荷時(shí)間歷程曲線[6]，如圖1所示。

由于飛行器發(fā)射時(shí)有艇速存在，使飛行器迎水面

圖1 飛行器軸向表面壓強(qiáng)插值原理Fig.1 Principle of Axial Surface Pressure Interpolation for Aircraft P1～P4—實(shí)測(cè)壓強(qiáng)數(shù)據(jù)；PK—插值得到的壓強(qiáng)數(shù)據(jù)；xk，tk—測(cè)點(diǎn)位置和出水時(shí)刻

根據(jù)水動(dòng)壓強(qiáng)測(cè)壓點(diǎn)相對(duì)位置信息與飛行器結(jié)構(gòu)外形，建立柱坐標(biāo)系，利用柱坐標(biāo)描述測(cè)壓點(diǎn)的位置。將迎水面和背水面測(cè)壓點(diǎn)坐標(biāo)存成一個(gè)文件，作為坐標(biāo)輸入。對(duì)迎水面和背水面測(cè)壓點(diǎn)對(duì)應(yīng)的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行預(yù)處理，與彈道合并作為數(shù)據(jù)文件。水動(dòng)壓強(qiáng)載荷輸入數(shù)據(jù)預(yù)處理流程如圖2所示，從而到飛行器表面隨時(shí)間和空間變化的連續(xù)壓強(qiáng)場(chǎng)，如圖3所示。和背水面存在壓強(qiáng)差，因此同一截面迎水面和背水面壓強(qiáng)的分布和大小不同，對(duì)于周向壓強(qiáng)場(chǎng)，則采用余弦理論插值得到。

圖2 水動(dòng)壓強(qiáng)載荷數(shù)據(jù)預(yù)處理流程Fig.2 Preprocessing Flow of Hydrodynamic Pressure Load Date

圖3 飛行器表面時(shí)間-空間連續(xù)壓強(qiáng)場(chǎng)Fig.3 Time-space Continuous Surface Pressure Field of Aircraft

3 飛行器動(dòng)響應(yīng)分析

3.1 飛行器質(zhì)量分布

載荷計(jì)算時(shí)的一個(gè)基本的定義是對(duì)廣義動(dòng)能起作用的任何項(xiàng)必須在質(zhì)量矩陣中創(chuàng)建一個(gè)系數(shù)。另一個(gè)定義是任何廣義力 F與加速度項(xiàng)成比例，就產(chǎn)生了質(zhì)量矩陣M，即：

其中，加速度矢量u˙˙的每一個(gè)分量代表一個(gè)廣義自由度[7]。

載荷計(jì)算時(shí)，將飛行器簡(jiǎn)化為六自由度變質(zhì)量-變剛度梁模型，梁彎曲振動(dòng)的單元質(zhì)量矩陣 me可按如下公式計(jì)算：

式中 A為梁截面積；l為單元長(zhǎng)度。

對(duì)于質(zhì)量矩陣，在實(shí)際工程計(jì)算中，根據(jù)情況也可以把各節(jié)點(diǎn)分得的集中質(zhì)量直接填入總質(zhì)量對(duì)應(yīng)位置[8]。集中轉(zhuǎn)動(dòng)慣量按下式計(jì)算：

式中 me為該站上的集中質(zhì)量；R為飛行器半徑；h為處于梁端點(diǎn)時(shí)與相鄰點(diǎn)間距之半；處于中間點(diǎn)時(shí)取它相鄰的前后兩點(diǎn)距離之半[9]。

對(duì)于集中質(zhì)量的彈性梁，其廣義質(zhì)量為

式中 Mj為第j個(gè)模態(tài)的廣義質(zhì)量；mi為第i站的集中質(zhì)量。

3.2 飛行器剛度分布

飛行器簡(jiǎn)化成梁時(shí)，平衡方程為

式中eωP為功等效力；u為位移矩陣；P為壓力矩陣；Q為剪切力矩陣；eK為梁?jiǎn)卧膭偠染仃?，即?/p>

計(jì)算導(dǎo)數(shù)并積分后可得到梁的剛度矩陣為

式中 E為材料的彈性模量；G為材料的剪切模量。

3.3 飛行器有限元模型

通過飛行器的質(zhì)量和剛度分布特性可建立飛行器有限元模型。由于飛行器長(zhǎng)徑比較大，而且沿飛行器軸向各部段的剛度和質(zhì)量質(zhì)心是變化的，因此計(jì)算時(shí)需建立飛行器的變質(zhì)量-變剛度梁模型。梁?jiǎn)卧牡刃Ы孛娣e為0.5 m2，轉(zhuǎn)動(dòng)慣量IZZ為1 kg·m2。離散單元為 2節(jié)點(diǎn)六自由度無質(zhì)量歐拉梁?jiǎn)卧|(zhì)量按集中質(zhì)量處理，建立的飛行器變質(zhì)量-變剛度有限元模型如圖4所示。

圖4 飛行器變質(zhì)量-變剛度有限元模型Fig.4 Mass-variable and Stiffness-variable Beam Element Model of Aircraft

3.4 飛行器動(dòng)響應(yīng)計(jì)算

作用在飛行器表面的力有縱向力和橫向力，根據(jù)縱向力分布系數(shù)和橫向力分布系數(shù)等參數(shù)，通過波形插值和余弦插值可以得到飛行器水下運(yùn)動(dòng)過程中作用在飛行器節(jié)點(diǎn)上的縱向外力場(chǎng)和橫向外力場(chǎng)，從而采用Newmark提出的逐步積分法進(jìn)行結(jié)構(gòu)分析。對(duì)飛行器整個(gè)水下運(yùn)動(dòng)過程利用動(dòng)力學(xué)分析軟件進(jìn)行瞬態(tài)動(dòng)力學(xué)分析，從而計(jì)算得到飛行器各截面在水下運(yùn)動(dòng)過程中隨時(shí)間變化的力矩M和剪切力f。

對(duì)于一個(gè)n維自由度的振動(dòng)系統(tǒng)，其動(dòng)力學(xué)方程為

式中 M，C和K分別為系統(tǒng)質(zhì)量、阻尼和剛度矩陣；{˙x}，{x˙}，{x}和{f}分別為加速度、速度、位移和激勵(lì)力向量[10]。

其特征方程為

其中，特征值所對(duì)應(yīng)的特征向量{?r}是正交的，同時(shí){?r}對(duì)剛度矩陣K及質(zhì)量矩陣M也是正交的，即：

為將物理坐標(biāo)表示的動(dòng)力學(xué)方程（13）解耦，須將其轉(zhuǎn)換到模態(tài)坐標(biāo)系。

通過以上變換，原來耦合的動(dòng)力學(xué)方程式就成為解耦的以模態(tài)坐標(biāo)表示的模態(tài)方程：

通過求解式（16）表示的 N個(gè)獨(dú)立的模態(tài)坐標(biāo)下的動(dòng)力學(xué)方程，可得到模態(tài)坐標(biāo)下的各階模態(tài)坐標(biāo)向量 {ηr( t)}，將其代入式（14）便可得系統(tǒng)在物理坐標(biāo)系下的位移響應(yīng){x}，將其加載在有限元模型上進(jìn)行動(dòng)力學(xué)響應(yīng)分析，進(jìn)而可求得系統(tǒng)的各截面隨時(shí)間變化的彎矩和剪切力，結(jié)果如圖5和圖6所示。

圖5 飛行器某截面處的彎矩曲線Fig.5 Bending Moment at a Section of the Aircraft

圖6 飛行器某截面處的剪切力曲線Fig.6 Shear Force at a Section of the Aircraft

4 結(jié) 論

a）對(duì)于軸向離散的表面壓強(qiáng)，可采用波形插值法插值，可得到計(jì)算所需軸向的時(shí)間連續(xù)壓強(qiáng)場(chǎng)；對(duì)于周向壓強(qiáng)場(chǎng)，則利用余弦理論插值得到。

b）將迎水面和背水面測(cè)壓點(diǎn)坐標(biāo)合并、表面壓強(qiáng)測(cè)點(diǎn)的實(shí)驗(yàn)數(shù)據(jù)預(yù)處理后與彈道合并，作為坐標(biāo)和載荷的輸入值；再進(jìn)行編程計(jì)算，得到飛行器表面隨時(shí)間和空間變化的連續(xù)壓強(qiáng)場(chǎng)。

c）對(duì)于大質(zhì)量變剛度的飛行器，可根據(jù)其質(zhì)量-剛度分布特征，可建立變剛度-變質(zhì)量的梁模型代替三維模型進(jìn)行動(dòng)力學(xué)響應(yīng)分析。

d）將插值得到的隨時(shí)間和空間變化的表面壓強(qiáng)場(chǎng)，加載到飛行器有限元模型上進(jìn)行動(dòng)力學(xué)響應(yīng)分析，可得到飛行器各截面隨時(shí)間變化的彎矩、剪切力，其結(jié)果對(duì)于發(fā)射條件選擇和飛行器結(jié)構(gòu)設(shè)計(jì)具有重要的參考價(jià)值。