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

（哈爾濱工業(yè)大學(xué) 航天學(xué)院，哈爾濱 150001）

1 引 言

超空泡減阻是目前實(shí)現(xiàn)水下航行體超高速運(yùn)動(dòng)的新原理、新途徑。利用超空泡減阻效應(yīng)，水下航行體可以達(dá)到非常高的速度。超空泡射彈發(fā)射時(shí)會(huì)受到擾動(dòng)，導(dǎo)致射彈在前進(jìn)的同時(shí)，還繞著其頭部來(lái)回轉(zhuǎn)動(dòng)，使射彈的尾部與空泡壁面發(fā)生連續(xù)的碰撞和反彈（尾拍）[1]。極大的軸向沖擊載荷和劇烈的尾拍使得超空泡射彈的結(jié)構(gòu)響應(yīng)與傳統(tǒng)的航行體有很大的不同。過(guò)大的動(dòng)應(yīng)力會(huì)導(dǎo)致射彈結(jié)構(gòu)失效，過(guò)大的彎曲變形可以使射彈運(yùn)動(dòng)時(shí)發(fā)生翻滾，導(dǎo)致運(yùn)動(dòng)失穩(wěn)[2]。因此研究超空泡射彈的結(jié)構(gòu)響應(yīng)具有重要意義，并為射彈的結(jié)構(gòu)優(yōu)化設(shè)計(jì)提供前提條件。

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

本文將彈體作為彈性體，并利用彈體與流體間的流固耦合作用來(lái)處理尾拍力，研究了依靠慣性作自由飛行的射彈的結(jié)構(gòu)響應(yīng)。先在Logvinovich獨(dú)立膨脹原理基礎(chǔ)上對(duì)流體域給出了一種簡(jiǎn)化計(jì)算模型，然后采用有限元法對(duì)射彈的結(jié)構(gòu)響應(yīng)進(jìn)行仿真計(jì)算。并分析了射彈的長(zhǎng)細(xì)比對(duì)其速度響應(yīng)、動(dòng)應(yīng)力響應(yīng)和彎曲變形的影響。

2 射彈模型與受力分析

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

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

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

由于運(yùn)動(dòng)的相對(duì)性，彈體模型的前端面中心點(diǎn)o看作靜止不動(dòng)[4]，在20 rad/s的初始擾動(dòng)角速度的作用下，模型繞該中心在鉛垂面內(nèi)做上下擺動(dòng)，尾部與空泡壁面不斷發(fā)生碰撞。

3 流體域的計(jì)算模型

3.1 流體域的簡(jiǎn)化計(jì)算模型

由于多相流CFD計(jì)算以及空泡流與彈體非線性振動(dòng)間的耦合計(jì)算耗時(shí)巨大，本文使用了流體域的簡(jiǎn)化計(jì)算模型，可以顯著降低計(jì)算的時(shí)間耗費(fèi)。

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

根據(jù)Logvinovich的獨(dú)立膨脹原理[9]，超空泡的每個(gè)橫截面在其所在的平面內(nèi)由中心幾乎獨(dú)立地向外擴(kuò)張，而在射彈飛行方向（縱向）上幾乎沒(méi)有速度，這已經(jīng)為試驗(yàn)所證實(shí)[10]。又考慮到尾拍的瞬時(shí)性，所以尾拍過(guò)程可以看作是運(yùn)動(dòng)著的彈體與靜止的水面之間的拍擊過(guò)程。也可以看作是，彈體模型做定軸轉(zhuǎn)動(dòng)時(shí)其尾部與高速前進(jìn)的水面間的撞擊。所以，通過(guò)圓筒形的運(yùn)動(dòng)流體來(lái)代替射彈周?chē)乃?、用圓筒的內(nèi)表面來(lái)代替空泡壁面（見(jiàn)圖2中的a圖），就完全可以模擬模型的尾拍過(guò)程。

3.2 流體域尺寸與邊界條件

流體域長(zhǎng)度0.74 m，內(nèi)徑0.022 m，外徑0.2 m，其內(nèi)徑由超空泡半徑的經(jīng)驗(yàn)公式來(lái)確定[11]。可以算得本文中射彈尾部的空泡半徑變化很小，故采用0.022 m作為流體域內(nèi)徑。

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

4 數(shù)值方法

射彈模型采用拉格朗日算法，流體域采用ALE有限元算法，可以實(shí)現(xiàn)流固耦合的動(dòng)態(tài)分析。ALE有限元矩陣方程[12]：

（1）質(zhì)量守恒方程

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

（2）動(dòng)量守恒方程

式中：M和L分別為廣義質(zhì)量和傳遞矩陣，vr為對(duì)應(yīng)于在參考構(gòu)形描述下的速度；fint和fext分別為內(nèi)力和外力向量。

本文采用LS-DYNA軟件進(jìn)行仿真計(jì)算。彈體采用全積分S/R六面體單元，算法為拉格朗日算法。流體域采用中心單點(diǎn)積分的帶空白材料的單物質(zhì)ALE單元，將空泡內(nèi)部空間設(shè)為空白材料，使用歐拉算法。彈體與流體間的耦合采用歐拉—拉格朗日耦合算法。單元網(wǎng)格采用六面體結(jié)構(gòu)化網(wǎng)格，固體單元與流體單元之間采用無(wú)侵蝕的罰函數(shù)耦合。為避免發(fā)生滲透現(xiàn)象對(duì)彈體尾部處的流體網(wǎng)格進(jìn)行了局部加細(xì)。

5 計(jì)算結(jié)果與分析

射彈的長(zhǎng)細(xì)比β分別取為：8.50、10.48、12.57、15.22、17.54和20。對(duì)彈體模型的速度響應(yīng)、動(dòng)應(yīng)力響應(yīng)和橫向彎曲變形進(jìn)行了數(shù)值計(jì)算，并對(duì)結(jié)果進(jìn)行了比較分析。

5.1 射彈的速度響應(yīng)

對(duì)彈體中部一節(jié)點(diǎn)的速度響應(yīng)進(jìn)行快速傅立葉變換，得到其頻譜圖，將頻譜圖中幅值最大的頻率記為K。然后對(duì)不同長(zhǎng)細(xì)比模型的速度響應(yīng)的K值進(jìn)行比較，可得到圖3所示的模型長(zhǎng)細(xì)比β與頻率K的關(guān)系圖。

從圖中可以看出，隨著β的增加，K逐漸降低，即彈體振動(dòng)的頻率變小了。β的增加使彈體變得更細(xì)長(zhǎng)，抗彎剛度降低，使結(jié)構(gòu)振動(dòng)響應(yīng)的周期也隨之變長(zhǎng)，從而使K值變小。

5.2 射彈的動(dòng)應(yīng)力響應(yīng)

對(duì)射彈模型的第一主應(yīng)力σ1、第三主應(yīng)力σ3、最大剪切應(yīng)力σs和Von Mises等效應(yīng)力σv-m進(jìn)行計(jì)算。結(jié)果顯示，在每個(gè)單元內(nèi)，這些應(yīng)力的大小都近似于同步變化，并同時(shí)達(dá)到各自的最大值。對(duì)于彈體內(nèi)的應(yīng)力分布，由于彈體中部彎曲變形最顯著，導(dǎo)致各應(yīng)力的最大值也都出現(xiàn)在彈體中部位置，如圖4所示。

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

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

4.3 射彈的橫向彎曲變形

建立射彈模型的局部坐標(biāo)系，坐標(biāo)原點(diǎn)位于前端面中心點(diǎn)O，坐標(biāo)系隨模型一起轉(zhuǎn)動(dòng)，X軸始終由O指向彈體尾端中心點(diǎn)，Z軸與對(duì)稱面垂直。如圖所示。則彈體的橫向彎曲變形大小可用局部坐標(biāo)系內(nèi)Y方向的位移量Dy來(lái)衡量。

結(jié)果顯示，彈體橫向變形最大的位置在彈體的中部，如圖4所示的n6-n8截面處。以β=10.48的彈體為例，當(dāng)t=0.006 4 s時(shí)，Dy達(dá)到最大值，在該時(shí)刻，彈體軸線的彎曲變形情況如圖7。對(duì)于其他長(zhǎng)細(xì)比的模型，當(dāng)出現(xiàn)最大彎曲時(shí)，彈體變形情況也是如此。

n6-n8截面處的Dy隨時(shí)間的變化如圖8所示。從圖中可見(jiàn)，彈體的橫向彎曲變形在運(yùn)動(dòng)開(kāi)始不久就達(dá)到極值，然后很快衰減。這與彈體所受外載荷的變化情況相符合。

彈體內(nèi)Dy的最大值與彈體長(zhǎng)細(xì)比的關(guān)系如圖9所示?？梢?jiàn)，隨著長(zhǎng)細(xì)比增加，彈體變得更細(xì)長(zhǎng)，從而使其彎曲變形有明顯增大的趨勢(shì)?？梢?jiàn)，較細(xì)長(zhǎng)的彈體飛行時(shí)的抖動(dòng)更嚴(yán)重，這將使射彈的運(yùn)動(dòng)更加不穩(wěn)定。

6 結(jié) 論

本文通過(guò)對(duì)射彈-流體的流固耦合仿真分析，得到如下結(jié)論：

（1）彈體的長(zhǎng)細(xì)比變大將使其速度響應(yīng)的頻率變小，彈體結(jié)構(gòu)振動(dòng)的頻率降低。

（2）彈體的第一主應(yīng)力、第三主應(yīng)力、最大剪切應(yīng)力和Von Mises等效應(yīng)力同時(shí)達(dá)到最大值，最大值出現(xiàn)在彈體中部。在運(yùn)動(dòng)的初期，各應(yīng)力達(dá)到極值，然后迅速衰減。彈體長(zhǎng)細(xì)比的變化對(duì)應(yīng)力最大值的影響較復(fù)雜，應(yīng)將長(zhǎng)細(xì)比分段考慮。

（3）彈體橫向變形最大的位置在彈體的中部。射彈運(yùn)動(dòng)開(kāi)始不久橫向變形達(dá)到最大值。隨著長(zhǎng)細(xì)比增加，彎曲變形程度明顯增加。

[1]曹 偉,王 聰,魏英杰等．自然超空泡形態(tài)特性的射彈試驗(yàn)研究[J].工程力學(xué),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]吳三靈,溫 波,于永強(qiáng)等.火炮動(dòng)力學(xué)實(shí)驗(yàn)[M].北京:國(guó)防工業(yè)出版社,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]孟慶昌,張志宏,顧建農(nóng)等.超空泡射彈尾拍分析與計(jì)算[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]李 政,金先龍,亓文果.流體-結(jié)構(gòu)耦合問(wèn)題的有限元并行計(jì)算研究[J].計(jì)算力學(xué)學(xué)報(bào),2007,24(6):727-732.