張勁生,張嘉鐘,王 聰,魏英杰
(哈爾濱工業(yè)大學(xué) 航天學(xué)院,哈爾濱 150001)
超空泡減阻是目前實(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)和彎曲變形的影響。
射彈模型前部為平頭圓錐體,平頭直徑為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ā)生碰撞。
由于多相流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ò)程。
流體域長(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ù)。
射彈模型采用拉格朗日算法,流體域采用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ì)。
射彈的長(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)行了比較分析。
對(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值變小。
對(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)系還是可以確定的。
建立射彈模型的局部坐標(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)定。
本文通過(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.