基于分布水動壓力的飛行器水下載荷計算方法研究

劉文一，李玉龍，吳訓濤

（中國人民解放軍91550部隊41分隊，大連，116023）

0 引 言

飛行器水下垂直向上運動時，由于空泡的出現，會經受復雜的水下載荷，在飛行器表面布置壓強傳感器來測量表面壓強，由于受傳感器安裝位置和采樣頻率的限制，測得的壓強數據在時間和空間上是離散的。為了對飛行器進行動力學響應分析，必須將離散的數據連續化，才能獲得飛行器的載荷分布情況[1]。

獲得飛行器表面壓強場后，將表面壓強場加載在飛行器表面上進動力學響應分析，以獲得飛行器對表面壓強的動響應，從而辨識得到飛行器的內力，即飛行器各截面隨時間和空間變化的軸向力、剪切力和彎矩[2]。對于長徑比較大的飛行器，剪切力和彎矩是其水下載荷的主要部分，這兩個值直接影響飛行能否正常出水和出水姿態，因此辨識得到的飛行器水下載荷具有重要的工程應用價值[3]。

1 表面壓強

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

飛行器在水下任一點的表面壓強可表示為

式中 CP為壓強系數；ρ為水的密度，kg/m2；V∞為飛行器運動速度，m/s；P0為環境壓強，Pa。

由于飛行器在水下會產生空化現象，因此載荷計算關心的是CP為負值的部分，引入符號ζ：

2 連續壓強場

飛行器體表面壓強測量數據測得了每個測點上壓強隨時間的變化歷程，但對于兩個測點之間的壓強分布，則必須對壓強測量數據進行連續化處理，利用波形插值法可以將離散的壓強數據插值得到隨時間變化的壓強場，這是因為在出水過程中，空泡沿飛行器表面是按照一定速度推進潰滅的[5]。

進行波形插值時，為確定某點的壓強載荷，首先根據其位置確定出水時刻，移動時間軸，將出水時刻的主波峰對齊，然后進行線性插值得到峰值附近的壓強曲線；最后還原時間軸，得到該點的壓強載荷時間歷程曲線[6]，如圖1所示。

由于飛行器發射時有艇速存在，使飛行器迎水面

圖1 飛行器軸向表面壓強插值原理Fig.1 Principle of Axial Surface Pressure Interpolation for Aircraft P1～P4—實測壓強數據；PK—插值得到的壓強數據；xk，tk—測點位置和出水時刻

根據水動壓強測壓點相對位置信息與飛行器結構外形，建立柱坐標系，利用柱坐標描述測壓點的位置。將迎水面和背水面測壓點坐標存成一個文件，作為坐標輸入。對迎水面和背水面測壓點對應的實驗數據進行預處理，與彈道合并作為數據文件。水動壓強載荷輸入數據預處理流程如圖2所示，從而到飛行器表面隨時間和空間變化的連續壓強場，如圖3所示。和背水面存在壓強差，因此同一截面迎水面和背水面壓強的分布和大小不同，對于周向壓強場，則采用余弦理論插值得到。

圖2 水動壓強載荷數據預處理流程Fig.2 Preprocessing Flow of Hydrodynamic Pressure Load Date

圖3 飛行器表面時間-空間連續壓強場Fig.3 Time-space Continuous Surface Pressure Field of Aircraft

3 飛行器動響應分析

3.1 飛行器質量分布

載荷計算時的一個基本的定義是對廣義動能起作用的任何項必須在質量矩陣中創建一個系數。另一個定義是任何廣義力 F與加速度項成比例，就產生了質量矩陣M，即：

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

載荷計算時，將飛行器簡化為六自由度變質量-變剛度梁模型，梁彎曲振動的單元質量矩陣 me可按如下公式計算：

式中 A為梁截面積；l為單元長度。

對于質量矩陣，在實際工程計算中，根據情況也可以把各節點分得的集中質量直接填入總質量對應位置[8]。集中轉動慣量按下式計算：

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

對于集中質量的彈性梁，其廣義質量為

式中 Mj為第j個模態的廣義質量；mi為第i站的集中質量。

3.2 飛行器剛度分布

飛行器簡化成梁時，平衡方程為

式中eωP為功等效力；u為位移矩陣；P為壓力矩陣；Q為剪切力矩陣；eK為梁單元的剛度矩陣，即：

計算導數并積分后可得到梁的剛度矩陣為

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

3.3 飛行器有限元模型

通過飛行器的質量和剛度分布特性可建立飛行器有限元模型。由于飛行器長徑比較大，而且沿飛行器軸向各部段的剛度和質量質心是變化的，因此計算時需建立飛行器的變質量-變剛度梁模型。梁單元的等效截面積為0.5 m2，轉動慣量IZZ為1 kg·m2。離散單元為 2節點六自由度無質量歐拉梁單元，質量按集中質量處理，建立的飛行器變質量-變剛度有限元模型如圖4所示。

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

3.4 飛行器動響應計算

作用在飛行器表面的力有縱向力和橫向力，根據縱向力分布系數和橫向力分布系數等參數，通過波形插值和余弦插值可以得到飛行器水下運動過程中作用在飛行器節點上的縱向外力場和橫向外力場，從而采用Newmark提出的逐步積分法進行結構分析。對飛行器整個水下運動過程利用動力學分析軟件進行瞬態動力學分析，從而計算得到飛行器各截面在水下運動過程中隨時間變化的力矩M和剪切力f。

對于一個n維自由度的振動系統，其動力學方程為

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

其特征方程為

其中，特征值所對應的特征向量{?r}是正交的，同時{?r}對剛度矩陣K及質量矩陣M也是正交的，即：

為將物理坐標表示的動力學方程（13）解耦，須將其轉換到模態坐標系。

通過以上變換，原來耦合的動力學方程式就成為解耦的以模態坐標表示的模態方程：

通過求解式（16）表示的 N個獨立的模態坐標下的動力學方程，可得到模態坐標下的各階模態坐標向量 {ηr( t)}，將其代入式（14）便可得系統在物理坐標系下的位移響應{x}，將其加載在有限元模型上進行動力學響應分析，進而可求得系統的各截面隨時間變化的彎矩和剪切力，結果如圖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 結 論

a）對于軸向離散的表面壓強，可采用波形插值法插值，可得到計算所需軸向的時間連續壓強場；對于周向壓強場，則利用余弦理論插值得到。

b）將迎水面和背水面測壓點坐標合并、表面壓強測點的實驗數據預處理后與彈道合并，作為坐標和載荷的輸入值；再進行編程計算，得到飛行器表面隨時間和空間變化的連續壓強場。

c）對于大質量變剛度的飛行器，可根據其質量-剛度分布特征，可建立變剛度-變質量的梁模型代替三維模型進行動力學響應分析。

d）將插值得到的隨時間和空間變化的表面壓強場，加載到飛行器有限元模型上進行動力學響應分析，可得到飛行器各截面隨時間變化的彎矩、剪切力，其結果對于發射條件選擇和飛行器結構設計具有重要的參考價值。