超空泡射彈結構響應的流固耦合仿真研究

張勁生，張嘉鐘，王 聰，魏英杰

（哈爾濱工業大學 航天學院，哈爾濱 150001）

1 引 言

超空泡減阻是目前實現水下航行體超高速運動的新原理、新途徑。利用超空泡減阻效應，水下航行體可以達到非常高的速度。超空泡射彈發射時會受到擾動，導致射彈在前進的同時，還繞著其頭部來回轉動，使射彈的尾部與空泡壁面發生連續的碰撞和反彈（尾拍）[1]。極大的軸向沖擊載荷和劇烈的尾拍使得超空泡射彈的結構響應與傳統的航行體有很大的不同。過大的動應力會導致射彈結構失效，過大的彎曲變形可以使射彈運動時發生翻滾，導致運動失穩[2]。因此研究超空泡射彈的結構響應具有重要意義，并為射彈的結構優化設計提供前提條件。

由于射彈空間體積狹小，同時極強的瞬態沖擊載荷對傳感器等測量器件又有潛在的破壞作用，因此用試驗測試方法，研究射彈的結構響應的難度比較大[3]。運用理論分析和數值計算，文獻[4-6]對存在尾拍時，航行體的動力學特性進行了研究，將航行體視為剛體，并對尾拍力的形式進行了簡化處理；文獻[7]將勻速運動的航行體作為彈性體，研究其結構響應，但對尾拍載荷的處理上與文獻[4-5]是一樣的。已有的對超空泡射彈結構響應的研究，都先對尾拍力采用高度簡化的數學關系式，忽略彈體與水的耦合作用，將彈體視為剛體或者讓彈體作勻速運動。

本文將彈體作為彈性體，并利用彈體與流體間的流固耦合作用來處理尾拍力，研究了依靠慣性作自由飛行的射彈的結構響應。先在Logvinovich獨立膨脹原理基礎上對流體域給出了一種簡化計算模型，然后采用有限元法對射彈的結構響應進行仿真計算。并分析了射彈的長細比對其速度響應、動應力響應和彎曲變形的影響。

2 射彈模型與受力分析

射彈模型前部為平頭圓錐體，平頭直徑為0.002 m，其后部為圓柱體。采用體積為3.4×10-6m3，但長細比不同的系列模型。材料為鋼，彈性模量為1.93×1011Pa，泊松比0.3，密度7750 kg/m3，且只考慮彈性變形。

本文不考慮空泡內氣體阻力和彈體重力的影響；彈體載荷包括前端面上的流體阻力FR和尾部的尾拍沖擊力FI；在慣性作用下，彈體作減速自由飛行，同時發生尾拍。如圖1所示。

FR方向沿模型軸向，其大小為[8]：FR=ρAcDv2/2，其中ρ，cD分別為流體的密度和阻力系數，A，v分別為模型前端面積和前進速度。模型的加速度a和慣性力FT由下式確定：FT=－FR=－ma。

由于運動的相對性，彈體模型的前端面中心點o看作靜止不動[4]，在20 rad/s的初始擾動角速度的作用下，模型繞該中心在鉛垂面內做上下擺動，尾部與空泡壁面不斷發生碰撞。

3 流體域的計算模型

3.1 流體域的簡化計算模型

由于多相流CFD計算以及空泡流與彈體非線性振動間的耦合計算耗時巨大，本文使用了流體域的簡化計算模型，可以顯著降低計算的時間耗費。

尾拍過程使流體域的內邊界（空泡壁面）不斷發生變化，使該處的邊界條件很難設定。所以本文將整個流體域設定為多物質ALE算法，只需要對流體域的外邊界進行邊界條件的設定。將空泡的內部空間設為空白物質單元。如圖2（a）所示。

根據Logvinovich的獨立膨脹原理[9]，超空泡的每個橫截面在其所在的平面內由中心幾乎獨立地向外擴張，而在射彈飛行方向（縱向）上幾乎沒有速度，這已經為試驗所證實[10]。又考慮到尾拍的瞬時性，所以尾拍過程可以看作是運動著的彈體與靜止的水面之間的拍擊過程。也可以看作是，彈體模型做定軸轉動時其尾部與高速前進的水面間的撞擊。所以，通過圓筒形的運動流體來代替射彈周圍的水流、用圓筒的內表面來代替空泡壁面（見圖2中的a圖），就完全可以模擬模型的尾拍過程。

3.2 流體域尺寸與邊界條件

流體域長度0.74 m，內徑0.022 m，外徑0.2 m，其內徑由超空泡半徑的經驗公式來確定[11]?？梢运愕帽疚闹猩鋸椢膊康目张莅霃阶兓苄?，故采用0.022 m作為流體域內徑。

由于對稱性，流體域只取一半來計算。對稱面為對稱邊界條件，其他邊界設定見圖2（b）所示。邊界速度的變化區間取為300～1400 m/s，其大小隨時間變化，并按下式確定[4]：v（t）=v0/（1+ αv0t ），其中，v0，α分別為射彈的初始速度和系統常數。

4 數值方法

射彈模型采用拉格朗日算法，流體域采用ALE有限元算法，可以實現流固耦合的動態分析。ALE有限元矩陣方程[12]：

（1）質量守恒方程

式中：Mρ，Lρ和 Kρ分別為容量,轉換和散度矩陣。

（2）動量守恒方程

式中：M和L分別為廣義質量和傳遞矩陣，vr為對應于在參考構形描述下的速度；fint和fext分別為內力和外力向量。

本文采用LS-DYNA軟件進行仿真計算。彈體采用全積分S/R六面體單元，算法為拉格朗日算法。流體域采用中心單點積分的帶空白材料的單物質ALE單元，將空泡內部空間設為空白材料，使用歐拉算法。彈體與流體間的耦合采用歐拉—拉格朗日耦合算法。單元網格采用六面體結構化網格，固體單元與流體單元之間采用無侵蝕的罰函數耦合。為避免發生滲透現象對彈體尾部處的流體網格進行了局部加細。

5 計算結果與分析

射彈的長細比β分別取為：8.50、10.48、12.57、15.22、17.54和20。對彈體模型的速度響應、動應力響應和橫向彎曲變形進行了數值計算，并對結果進行了比較分析。

5.1 射彈的速度響應

對彈體中部一節點的速度響應進行快速傅立葉變換，得到其頻譜圖，將頻譜圖中幅值最大的頻率記為K。然后對不同長細比模型的速度響應的K值進行比較，可得到圖3所示的模型長細比β與頻率K的關系圖。

從圖中可以看出，隨著β的增加，K逐漸降低，即彈體振動的頻率變小了。β的增加使彈體變得更細長，抗彎剛度降低，使結構振動響應的周期也隨之變長，從而使K值變小。

5.2 射彈的動應力響應

對射彈模型的第一主應力σ1、第三主應力σ3、最大剪切應力σs和Von Mises等效應力σv-m進行計算。結果顯示，在每個單元內，這些應力的大小都近似于同步變化，并同時達到各自的最大值。對于彈體內的應力分布，由于彈體中部彎曲變形最顯著，導致各應力的最大值也都出現在彈體中部位置，如圖4所示。

該位置位于彈體中部的外表面上，通過應力分析可知，該處的應力狀態可近似為單向應力狀態，橫截面近似為應力主平面。當該處處于拉伸變形時，σ1＞0和σ3=0；當處于壓縮時，σ1=0和σ3＜0。仿真得到該位置的σ1隨時間的變化如圖5所示。從圖中可以看出，在彈體上下擺動的一個周期內，有近一半的時間σ1＞0，另一半時間σ1=0，這與彈體的動態彎曲變形過程是一致的。并且，在運動的初期，第一主應力σ1達到最大值，然后迅速減小。在該位置上，應力、σs和σv-m的變化情況也與之類似。

對于長細比不同的模型，隨著β的增加，彈體內σ1、、σs和σv-m的最大值的變化較復雜，忽高忽低。如圖6所示。這主要是由于，長細比的改變不僅僅使彈體的變形能力發生了變化，而且對尾拍載荷也產生了很大的影響，這些影響因素的共同作用使最大應力與長細比β的關系變得復雜、不明確。但從圖中可以看出，如果將β分段來看，在β不同的取值范圍內，長細比與最大應力的關系還是可以確定的。

4.3 射彈的橫向彎曲變形

建立射彈模型的局部坐標系，坐標原點位于前端面中心點O，坐標系隨模型一起轉動，X軸始終由O指向彈體尾端中心點，Z軸與對稱面垂直。如圖所示。則彈體的橫向彎曲變形大小可用局部坐標系內Y方向的位移量Dy來衡量。

結果顯示，彈體橫向變形最大的位置在彈體的中部，如圖4所示的n6-n8截面處。以β=10.48的彈體為例，當t=0.006 4 s時，Dy達到最大值，在該時刻，彈體軸線的彎曲變形情況如圖7。對于其他長細比的模型，當出現最大彎曲時，彈體變形情況也是如此。

n6-n8截面處的Dy隨時間的變化如圖8所示。從圖中可見，彈體的橫向彎曲變形在運動開始不久就達到極值，然后很快衰減。這與彈體所受外載荷的變化情況相符合。

彈體內Dy的最大值與彈體長細比的關系如圖9所示。可見，隨著長細比增加，彈體變得更細長，從而使其彎曲變形有明顯增大的趨勢?？梢?，較細長的彈體飛行時的抖動更嚴重，這將使射彈的運動更加不穩定。

6 結 論

本文通過對射彈-流體的流固耦合仿真分析，得到如下結論：

（1）彈體的長細比變大將使其速度響應的頻率變小，彈體結構振動的頻率降低。

（2）彈體的第一主應力、第三主應力、最大剪切應力和Von Mises等效應力同時達到最大值，最大值出現在彈體中部。在運動的初期，各應力達到極值，然后迅速衰減。彈體長細比的變化對應力最大值的影響較復雜，應將長細比分段考慮。

（3）彈體橫向變形最大的位置在彈體的中部。射彈運動開始不久橫向變形達到最大值。隨著長細比增加，彎曲變形程度明顯增加。

[1]曹 偉,王 聰,魏英杰等．自然超空泡形態特性的射彈試驗研究[J].工程力學,2006,23(12):175-179.

[2]Hrubes J D.High-speed imaging of supercavitating underwater projectiles[J].Experiments in Fluids,2001,30:57-64.

[3]吳三靈,溫 波,于永強等.火炮動力學實驗[M].北京:國防工業出版社,2004.

[4]Rand R,Pratap R,Ramani D.Impact Dynamics of a supercavitating underwater projectile[C]//Proceedings of the 1997 ASME Design Engineering Technical Conferences.Sacramento:ASME,1997.

[5]Kulkarni S S,Pratap R.Studies on the dynamics of a supercavitating projectile[J].Applied Mathematical Modelling,2000,24:113-129.

[6]孟慶昌,張志宏,顧建農等.超空泡射彈尾拍分析與計算[J].爆炸與沖擊,2009,29(1):57-60.

[7]Ruzzene M,Soranna F.Impact dynamics of elastic supercavitating underwater vehicles[C]//9th AIAA/ISSMO Symposium on Multidisciplinary Analysis and Optimization.Atlanta:AIAA,2002.

[8]Kirschner I N,Kring D C,Stokes A W,et al.Control strategies for supercavitating vehicles[J].Journal of Vibration and Control,2002,8:219-242.

[9]Vasin A D.The Principle of independence of the cavity sections expansion(Logvinovich’s principle)as the basis for investigation on cavitation flows[C]//VKI Special Courses on Supercavitating Flows.Brussels:RTO-AVT and VKI,2001,RTO-EN-010(8).

[10]Semenenko V N.Artificial supercavitation:Physics and calculation[C]//VKI Special Courses on Supercavitating Flows.Brussels:RTO-AVT and VKI,2001,RTO-EN-010(11).

[11]Savchenko Y N.Experimental investigation of supercavitating motion of bodies[C]//VKI Special Courses on Supercavitating Flows.Brussels:RTO-AVT and VKI,2001,RTO-EN-010(4).

[12]李 政,金先龍,亓文果.流體-結構耦合問題的有限元并行計算研究[J].計算力學學報,2007,24(6):727-732.