平頭錐型回轉(zhuǎn)體高速入水結(jié)構(gòu)強(qiáng)度數(shù)值分析*

黃志剛，孫鐵志，楊碧野，張桂勇,2,3，宗 智,2,3

（1.大連理工大學(xué)船舶工程學(xué)院遼寧省深海浮動(dòng)結(jié)構(gòu)工程實(shí)驗(yàn)室，遼寧 大連 116024；2.大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024；3.高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240）

回轉(zhuǎn)體從空中高速跨越空氣和水交界面時(shí)，受到水介質(zhì)對(duì)其作用的巨大沖擊載荷，從而引起回轉(zhuǎn)體發(fā)生形變和破壞，這就需要設(shè)計(jì)的回轉(zhuǎn)體強(qiáng)度滿足結(jié)構(gòu)安全要求。同時(shí)由于回轉(zhuǎn)體入水過程涉及氣-液-固三相耦合作用[1]，在高速條件下這種流固耦合作用變得更復(fù)雜，從而大大增加了研究高速回轉(zhuǎn)體入水問題的難度。因此，開展回轉(zhuǎn)體高速入水過程結(jié)構(gòu)強(qiáng)度研究對(duì)回轉(zhuǎn)體結(jié)構(gòu)設(shè)計(jì)具有重要意義。

入水問題研究長期以來一直受到廣泛關(guān)注。最早Worthington等[2]利用閃光攝影技術(shù)，觀測(cè)到了小球落水時(shí)的液體飛濺和空泡現(xiàn)象。在理論研究上，1929年，Von Karman[3]最早提出使用附加質(zhì)量法計(jì)算了入水沖擊載荷，并采用動(dòng)量守恒定律推導(dǎo)了相關(guān)的計(jì)算公式。隨后Wagner[4]在此基礎(chǔ)上運(yùn)用伯努利方程，提出了平板理論，此理論在研究結(jié)構(gòu)入水過程受力和流體動(dòng)能方面得到了廣泛應(yīng)用。Mayo[5]考慮了波浪對(duì)附加質(zhì)量、壓力分布、入水沖擊載荷的影響。May[6]研究了導(dǎo)彈在垂直和傾斜入水過程中的入水空泡形成、發(fā)展和潰滅的過程，同時(shí)研究了入水過程的流場(chǎng)變化情況。近些年來，隨著數(shù)值計(jì)算能力的提高，有限元將流體和結(jié)構(gòu)耦合解決了大型三維幾何入水問題。Anghileri等[7]利用有限元分析了剛性球在垂直入水過程中所受的載荷。Faltinsen等[8]建立了入水問題的數(shù)值和理論模型，并通過相關(guān)的入水實(shí)驗(yàn)驗(yàn)證其模型的準(zhǔn)確性。Korobkin等[9]運(yùn)用邊界元和有限差分法計(jì)算懸浮體入水沖擊過程中引起的液流非定常問題。Donguy等[10]對(duì)二維、三維、剛性和彈性體的入水耦合過程進(jìn)行了詳細(xì)分析，其采用的是一種變化的有限元耦合方法，通過耦合矩陣處理流體和結(jié)構(gòu)之間的作用。

鄭金偉等[11]利用LS-DYNA程序計(jì)算了三維剛體結(jié)構(gòu)傾斜入水過程的流體壓力和沖擊載荷。施紅輝等[12]、Shi等[13]通過鈍體入水實(shí)驗(yàn)測(cè)量了水下聲場(chǎng)的變化，獲得了激波的空間能量分布以及傳播速度變化，同時(shí)其還研究了鈍體入水自由液面變形問題和超空泡現(xiàn)象。王云等[14]開展了模型實(shí)驗(yàn)，利用高速攝像機(jī)拍攝了回轉(zhuǎn)體入水過程的空泡形態(tài)演變，分析了頭型、入水角度、入水速度和水下彈道的影響。潘光等[15]利用MSC.DYTRAN軟件計(jì)算了入水過程中魚雷所受到的沖擊壓力以及壓力分布情況。黃凱等[16]通過實(shí)驗(yàn)研究了回轉(zhuǎn)體頭型對(duì)反彈水花的影響，發(fā)現(xiàn)回轉(zhuǎn)體頂角越大，空泡越大。李佳川等[17]基于入水彈道學(xué)和超空泡理論，研究了不同擾動(dòng)角速度下的回轉(zhuǎn)體入水軌跡、速度、俯仰角的變化。張偉等[18]、郭子濤等[19]開展了水平回轉(zhuǎn)體入水實(shí)驗(yàn)，研究了回轉(zhuǎn)體入水過程的彈道穩(wěn)定性和空泡拓展特性，計(jì)算了入水阻力系數(shù)，并將實(shí)驗(yàn)中回轉(zhuǎn)體的速度衰減曲線與理論解進(jìn)行對(duì)比，兩者吻合良好。馬慶鵬等[20]通過球體入水實(shí)驗(yàn)，分析了球體入水過程中的入水空泡演變過程，得到了不同的入水速度以及表面黏濕狀態(tài)對(duì)入水空泡流場(chǎng)的影響。

目前針對(duì)入水問題已開展了大量的研究，并取得了豐碩的成果，然而在研究回轉(zhuǎn)體速度達(dá)到100 m/s及以上的高速入水問題時(shí)，多以入水空泡演化以及回轉(zhuǎn)體入水穩(wěn)定性作為研究重點(diǎn)。由于高速入水問題是沖擊作用時(shí)間極短的非線性問題，對(duì)高速回轉(zhuǎn)體入水過程結(jié)構(gòu)強(qiáng)度方面的研究仍顯不足。本文中采用數(shù)值計(jì)算的方法，分析入水過程中的沖擊載荷特點(diǎn)，并開展回轉(zhuǎn)體在不同殼體厚度情況下高速入水過程結(jié)構(gòu)強(qiáng)度研究，以期獲得的研究成果可為高速回轉(zhuǎn)體模型結(jié)構(gòu)設(shè)計(jì)提供參考，同時(shí)豐富入水過程機(jī)理內(nèi)涵。

1 數(shù)值計(jì)算模型設(shè)置

1.1 流固耦合ALE算法及狀態(tài)方程

基于LS-DYNA通用非線性有限元程序，使用ALE(arbitrary Lagrangian-Eulerian)方法模擬回轉(zhuǎn)體入水，計(jì)算入水沖擊載荷、平均壓力以及速度等的變化，著重考慮此過程中回轉(zhuǎn)體的形變情況。

在計(jì)算中，首先進(jìn)行Lagrange時(shí)步計(jì)算，單元網(wǎng)格隨材料流動(dòng)開始變形，然后進(jìn)行ALE時(shí)步計(jì)算：

(1) 保持變形后物體邊界條件不變，內(nèi)部單元進(jìn)行重分網(wǎng)格，稱為光滑步。

(2) 變形后的網(wǎng)格中的單元變量（密度、能量、應(yīng)力張量等）和節(jié)點(diǎn)速度矢量輸運(yùn)到重新劃分的新網(wǎng)格中，稱為對(duì)流步。

在LS-DYNA中，流體介質(zhì)的壓力由狀態(tài)方程進(jìn)行描述，本文中研究三維回轉(zhuǎn)體-水-空氣相互耦合求解過程。對(duì)于空氣和水，采用LS-DYNA中的NULL材料，對(duì)空氣介質(zhì)壓力，使用多項(xiàng)式狀態(tài)方程描述，在LS-DYNA中通過*EOS_LINEAR_POLYNOMIAL關(guān)鍵字施加，狀態(tài)方程壓力公式用下式表示：

式中：pa為氣體壓力，C0～C6為自定義常數(shù)，Va為相對(duì)體積， μa為體積變化率。對(duì)理想氣體，C0～C3、C6均為0，C4=C5=γa-1， 其中 γa為單位熱值率,初始?xì)怏w內(nèi)能Ea0=253.3 kJ

對(duì)于水介質(zhì)，其壓力使用Grüneisen狀態(tài)方程描述，在LS-DYNA中所使用的關(guān)鍵字為*EOS_Grüneisen，水壓力方程表示為：

式中：pw為水壓力，cw為聲音在水中的傳播速度， μw為體積變化率，α為對(duì)Grüneisen系數(shù) γ0的一階體積修正，S1、S2、S3分別為us-up曲線斜率無量綱系數(shù)，us為沖擊波速度，up為流體質(zhì)點(diǎn)的速度，Ew0為材料初始內(nèi)能。模擬中水狀態(tài)方程參數(shù)：cw=1 647 m/s，S1=1.921，S2=-0.096,S3=0， γ0=0.35，Ew0=289.5 kJ，相對(duì)初始體積Vw0=1.0。

1.2 計(jì)算模型及網(wǎng)格劃分

建立計(jì)算域模型如圖1所示，計(jì)算域包括空氣域和水域?？諝庥虺叽鐬?.2 m×1.2 m×1.0 m，水域尺寸為1.2 m×1.2 m×1.5 m。水域和空氣域?qū)挾葹榛剞D(zhuǎn)體最大直徑的10倍。在模擬中采用了2種回轉(zhuǎn)體,其結(jié)構(gòu)形式分別為：(1)全回轉(zhuǎn)體等厚度；(b)回轉(zhuǎn)體頭部部分厚度為8 mm，后體區(qū)域賦予一種厚度。回轉(zhuǎn)體模型整體為平頭錐形結(jié)構(gòu)，最大直徑為120 mm，頭部直徑為100 mm，回轉(zhuǎn)體長度為0.508 m，其首端傾斜斜面與回轉(zhuǎn)體頭部平面所呈角度為104.7°，回轉(zhuǎn)體外形如圖2所示。

圖1 計(jì)算域Fig.1 Computational domain

圖2 回轉(zhuǎn)體幾何模型Fig.2 The geometrical model for a revolution body

在模型單元類型的選擇上，彈殼體采用4節(jié)點(diǎn)SHELL單元，水和空氣采用單元類型為LS-DYNA中模擬ALE物質(zhì)的SOLID_ALE六面體單元。在入水的初始時(shí)刻保證回轉(zhuǎn)體在空氣中垂直水面向下，頭部貼近水面位置。模擬中應(yīng)用了ALE多物質(zhì)算法，空氣和水為2種ALE物質(zhì)。結(jié)構(gòu)和流體通過*CONSTRAINED_LAGRANGE_IN_SOLID進(jìn)行耦合，同時(shí)使用*INITIAL_FRACTION_GEOMETRY關(guān)鍵字將回轉(zhuǎn)體的內(nèi)部物質(zhì)設(shè)置為空氣，這樣符合實(shí)際情況，模擬更準(zhǔn)確。

由于高速入水過程時(shí)間極短，重力對(duì)整個(gè)入水過程的影響很小，所以在模擬中忽略重力的作用。

在數(shù)值計(jì)算過程中，回轉(zhuǎn)體為彈塑性材料，外殼采用45鋼材，四邊形網(wǎng)格大小為0.005 m。水和空氣六面體網(wǎng)格大小為0.02 m，計(jì)算模擬中總網(wǎng)格量為520 000。45鋼材的密度為7 850 kg/m3，楊氏模量為210 GPa，泊松比為0.3，屈服應(yīng)力為355 MPa，塑性失效應(yīng)變?yōu)?.35。

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

2.1 數(shù)值方法驗(yàn)證

為了驗(yàn)證數(shù)值方法的準(zhǔn)確性，選取第1種工況下板厚為5 mm的回轉(zhuǎn)體進(jìn)行入水計(jì)算。在入水的瞬間回轉(zhuǎn)體頭部會(huì)受到巨大的沖擊載荷，其中心區(qū)域有最大的壓強(qiáng)峰值。通過數(shù)值模擬得到的回轉(zhuǎn)體中心區(qū)域在入水抨擊過程中的壓強(qiáng)曲線，如圖3所示，其抨擊壓強(qiáng)的最大值為189 MPa。參考王珂等[21]對(duì)彈性回轉(zhuǎn)體入水抨擊載荷峰值的預(yù)報(bào)，彈性回轉(zhuǎn)體入水抨擊壓強(qiáng)峰值與回轉(zhuǎn)體的厚度、入水速度以及材料的彈性模量有關(guān)，抨擊過程峰值壓強(qiáng)p可以表示為：

圖3 回轉(zhuǎn)體頭部中心區(qū)域壓強(qiáng)曲線Fig.3 Pressure intensity curve in the central head region of the revolution body

式中：kd=0.00423d2+0.003758d+0.42，kE=ln(E/Pa)+15，kv=+5； ρw為水的密度；kd為厚度相關(guān)系數(shù)，d為回轉(zhuǎn)體厚度(m)；kE為彈性模量相關(guān)系數(shù)，E為材料彈性模量；kv為速度相關(guān)系數(shù)，v0為回轉(zhuǎn)體入水的初速度(m/s)。由式(4)可以求得在5 mm板厚下的回轉(zhuǎn)體入水沖擊載荷壓強(qiáng)峰值為173 MPa，和本文求得的數(shù)值解比較接近，但數(shù)值模擬結(jié)果略大于公式求得的抨擊壓強(qiáng)峰值。其原因是在數(shù)值模擬中未考慮空氣墊的作用，即高速時(shí)的空氣壓縮效應(yīng)，這也在一定程度上解釋了數(shù)值模擬所得壓強(qiáng)數(shù)值偏大的現(xiàn)象，空氣墊效應(yīng)是數(shù)值模擬過程中較難考慮的一個(gè)因素，在數(shù)值模擬方面還有很大的研究空間，亟待將來的研究。

除將入水時(shí)回轉(zhuǎn)體頭部抨擊壓力的峰值與擬合公式相對(duì)比外，還進(jìn)行了回轉(zhuǎn)體入水過程中速度衰減的對(duì)比驗(yàn)證，為了更好地對(duì)比速度衰減過程，在此數(shù)值模擬中延長了入水過程的求解時(shí)間，增大水域深度至3.5 m，計(jì)算時(shí)間為0.11 s。參考文獻(xiàn)[18]，根據(jù)牛頓第二定律，忽略重力效應(yīng)。

假定回轉(zhuǎn)體在入水過程中空泡內(nèi)外壓差保持不變，同時(shí)空化數(shù)為非定值，有：

式中：Δp為回轉(zhuǎn)體入水空泡內(nèi)外壓差， ρa(bǔ)和 ρw分別為空氣和水的密度，v0為回轉(zhuǎn)體入水初速度，vp為回轉(zhuǎn)體在水中的速度；Ca為氣流壓力降因數(shù)，Ca=5～15； σ0為初始空化數(shù)，σ0=0.006～0.018。

柱形平頭回轉(zhuǎn)體阻力因數(shù)和空化數(shù)關(guān)系式：

忽略入水過程中重力的影響，回轉(zhuǎn)體有以下運(yùn)動(dòng)方程：

式中：mp為回轉(zhuǎn)體的質(zhì)量，A0為回轉(zhuǎn)體頭部面積。

根據(jù)式(5)～(6)得到速度衰減和時(shí)間的關(guān)系式：

式中：衰減系數(shù)k=ρwA0Cp/2mp。依據(jù)式(8)計(jì)算回轉(zhuǎn)體入水速度衰減曲線并與數(shù)值模擬結(jié)果進(jìn)行對(duì)比，如圖4所示。從圖4可以看出，數(shù)值模擬結(jié)果與理論速度衰減曲線吻合較好，從而進(jìn)一步驗(yàn)證了所建立的數(shù)值方法的有效性。

圖4 回轉(zhuǎn)體速度衰減曲線Fig.4 Velocity attenuation of the revolution body over time

圖5 流-構(gòu)耦合力曲線Fig.5 Fluid-structure interaction force varying with time

在入水?dāng)?shù)值驗(yàn)證工作的同時(shí)，圖5給出了當(dāng)前回轉(zhuǎn)體結(jié)構(gòu)形式下入水過程中回轉(zhuǎn)體頭部表面作用力曲線圖，可以進(jìn)一步分析回轉(zhuǎn)體入水過程受力特征。從圖5可以看出，回轉(zhuǎn)體頭部受到的作用力在入水瞬間急劇上升，達(dá)到峰值后迅速下降，最后以小幅度震蕩趨勢(shì)變化?？梢娫谌胨查g產(chǎn)生了巨大沖擊力，所以在設(shè)計(jì)回轉(zhuǎn)體模型時(shí)要重點(diǎn)考慮入水瞬間的結(jié)構(gòu)強(qiáng)度。下面將對(duì)2種結(jié)構(gòu)形式的回轉(zhuǎn)體進(jìn)行強(qiáng)度的分析計(jì)算。

2.2 全回轉(zhuǎn)體等厚度工況模擬

在模擬中，首先進(jìn)行回轉(zhuǎn)體等厚度工況分析，分別計(jì)算回轉(zhuǎn)體厚度為2、3、4和5 mm的情況。計(jì)算中發(fā)現(xiàn)：只有回轉(zhuǎn)體厚度為5 mm時(shí)，該種結(jié)構(gòu)的回轉(zhuǎn)體是安全的；此種結(jié)構(gòu)的回轉(zhuǎn)體發(fā)生破壞時(shí)，破壞部位均發(fā)生在回轉(zhuǎn)體頭部和錐形斜面連接處。對(duì)回轉(zhuǎn)體厚度為2 mm的計(jì)算結(jié)果進(jìn)行分析，其中1.1 ms和3.0 ms時(shí)刻的應(yīng)力云圖分別如圖6和圖7所示。由圖6～7可以清晰地看出，破壞部位為頭部和錐形斜面連接處，回轉(zhuǎn)體頭部中心處未發(fā)生破壞現(xiàn)象。提取回轉(zhuǎn)體頭部和錐形斜面連接處的塑性應(yīng)變和單元應(yīng)力曲線，分別如圖8和圖9所示。由這2條曲線可以得到：在回轉(zhuǎn)體入水過程中，等效應(yīng)力達(dá)到了材料的屈服應(yīng)力355 MPa，塑性應(yīng)變也達(dá)到了材料的塑性失效應(yīng)變0.35，因此回轉(zhuǎn)體出現(xiàn)破壞現(xiàn)象?？梢姡剞D(zhuǎn)體頭部和錐形斜面連接處應(yīng)予以加強(qiáng)，來承受入水瞬間高數(shù)量級(jí)的沖擊載荷。同時(shí)因?yàn)槿胨^程中回轉(zhuǎn)體頭部為主要的沖擊載荷作用面，所以對(duì)回轉(zhuǎn)體頭部進(jìn)行加強(qiáng)是十分必要的。

圖6 在t=1.1 ms時(shí)回轉(zhuǎn)體的應(yīng)力分布Fig.6 Stress distribution in the revolution body at t=1.1 ms

圖7 在t=3.0 ms時(shí)回轉(zhuǎn)體的應(yīng)力分布Fig.7 Stress distribution in the revolution body at t=3.0 ms

圖8 回轉(zhuǎn)體頭部邊緣單元應(yīng)變曲線Fig.8 Strain-time curve of the edge element for the head of the revolution body

圖9 回轉(zhuǎn)體頭部邊緣單元等效應(yīng)力曲線Fig.9 Equivalent stress-time curve of the edge element for the head of the revolution body

2.3 回轉(zhuǎn)體頭部壁厚的影響

由2.2節(jié)中計(jì)算結(jié)果可知，回轉(zhuǎn)體頭部和錐形斜面連接處強(qiáng)度脆弱，因此對(duì)回轉(zhuǎn)體頭部進(jìn)行加強(qiáng)，考慮到工程安全系數(shù)要求，將頭部賦予8 mm的厚度，回轉(zhuǎn)體其余板厚度分別賦予2.0、2.5、3.0、3.5和4.0 mm。在初速度相同的條件下分別進(jìn)行計(jì)算，由計(jì)算結(jié)果可知：在頭部厚度為8 mm、后體區(qū)域厚度為2.0 mm時(shí)，回轉(zhuǎn)體無法滿足強(qiáng)度要求；在后體區(qū)域厚度為2.5 mm及以上時(shí)，結(jié)構(gòu)強(qiáng)度滿足要求。選取后體區(qū)域厚度為2.5 mm的情況加以分析?；剞D(zhuǎn)體入水過程中0.89 ms和2.00 ms時(shí)刻的應(yīng)力云圖分別如圖10和圖11所示，同時(shí)圖12給出了回轉(zhuǎn)體頭部邊緣和頭部中心單元塑性應(yīng)變曲線，圖13和圖14分別給出了回轉(zhuǎn)體頭部和頭部中心單元的等效應(yīng)力變化曲線。由應(yīng)力、應(yīng)變變化曲線可知，在回轉(zhuǎn)體入水過程中，無論其邊緣還是中心單元的等效應(yīng)力很快達(dá)到屈服應(yīng)力355 MPa，隨后等效應(yīng)力下降并表現(xiàn)為波動(dòng)變化?；剞D(zhuǎn)體頭部邊緣的塑性應(yīng)變由0迅速到達(dá)0.14，并最終趨近0.16，回轉(zhuǎn)體頭部中心的塑性應(yīng)變有同樣的變化規(guī)律并最終達(dá)到0.25。這個(gè)塑性應(yīng)變小于材料的塑性失效應(yīng)變0.35，因此，材料雖然進(jìn)入塑性階段，但還未達(dá)到破壞，此時(shí)回轉(zhuǎn)體結(jié)構(gòu)滿足強(qiáng)度要求。

圖10 t=0.89 ms時(shí)刻回轉(zhuǎn)體的應(yīng)力分布Fig.10 Stress distribution in the revolution body at t=0.89 ms

圖11 t=2.00 ms回轉(zhuǎn)體應(yīng)力圖Fig.11 Stress distribution in the revolution body at t=2.00 ms

圖12 回轉(zhuǎn)體頭部邊緣及中心單元應(yīng)變隨時(shí)間的變化Fig.12 Strain-time curves of the edge and central elements for the head of the revolution body

圖13 回轉(zhuǎn)體頭部邊緣單元等效應(yīng)力曲線Fig.13 Equivalent stress-time curve of the edge element for the head of the revolution body

圖14 回轉(zhuǎn)體頭部中心單元等效應(yīng)力曲線Fig.14 Equivalent stress-time curve of the central element for the head of the revolution body

3 結(jié) 論

采用有限元方法對(duì)平頭錐頭回轉(zhuǎn)體在高速入水過程中的結(jié)構(gòu)強(qiáng)度進(jìn)行了數(shù)值分析，得到的主要結(jié)論如下：

(1) 在回轉(zhuǎn)體入水的瞬間，其頭部受到高速?zèng)_擊載荷，且作用時(shí)間極短，在研究回轉(zhuǎn)體高速入水時(shí)要重點(diǎn)考慮入水瞬間的作用力。

(2) 在回轉(zhuǎn)體等厚度結(jié)構(gòu)形式下，等效應(yīng)力和塑性應(yīng)變均較大，在入水初速度為100 m/s、回轉(zhuǎn)體厚度小于5.0 mm時(shí)，回轉(zhuǎn)體應(yīng)力達(dá)到了材料的屈服應(yīng)力355 MPa，其塑性應(yīng)變也達(dá)到失效應(yīng)變0.35，此種形式下回轉(zhuǎn)體厚度低于5.0 mm時(shí)無法保證回轉(zhuǎn)體的強(qiáng)度要求。

(3) 當(dāng)回轉(zhuǎn)體入水初速度為100 m/s、其頭部厚度為8 mm、其后體壁厚為2.5 mm及以上時(shí)，結(jié)構(gòu)未發(fā)生破壞，此時(shí)回轉(zhuǎn)體可以滿足強(qiáng)度要求，但材料同樣進(jìn)入了塑性階段，建議對(duì)實(shí)際結(jié)構(gòu)還應(yīng)該考慮一定的安全系數(shù)。