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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

式中：pw為水壓力，cw為聲音在水中的傳播速度， μw為體積變化率，α為對Grüneisen系數(shù) γ0的一階體積修正，S1、S2、S3分別為us-up曲線斜率無量綱系數(shù)，us為沖擊波速度，up為流體質(zhì)點的速度，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，相對初始體積Vw0=1.0。

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

建立計算域模型如圖1所示，計算域包括空氣域和水域。空氣域尺寸為1.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 計算域Fig.1 Computational domain

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

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

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

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

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

2.1 數(shù)值方法驗證

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

圖3 回轉(zhuǎn)體頭部中心區(qū)域壓強曲線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)體入水沖擊載荷壓強峰值為173 MPa，和本文求得的數(shù)值解比較接近，但數(shù)值模擬結(jié)果略大于公式求得的抨擊壓強峰值。其原因是在數(shù)值模擬中未考慮空氣墊的作用，即高速時的空氣壓縮效應(yīng)，這也在一定程度上解釋了數(shù)值模擬所得壓強數(shù)值偏大的現(xiàn)象，空氣墊效應(yīng)是數(shù)值模擬過程中較難考慮的一個因素，在數(shù)值模擬方面還有很大的研究空間，亟待將來的研究。

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

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

式中：Δp為回轉(zhuǎn)體入水空泡內(nèi)外壓差， ρa和 ρ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)體有以下運動方程：

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

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

式中：衰減系數(shù)k=ρwA0Cp/2mp。依據(jù)式(8)計算回轉(zhuǎn)體入水速度衰減曲線并與數(shù)值模擬結(jié)果進(jìn)行對比，如圖4所示。從圖4可以看出，數(shù)值模擬結(jié)果與理論速度衰減曲線吻合較好，從而進(jì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ù)值驗證工作的同時，圖5給出了當(dāng)前回轉(zhuǎn)體結(jié)構(gòu)形式下入水過程中回轉(zhuǎn)體頭部表面作用力曲線圖，可以進(jìn)一步分析回轉(zhuǎn)體入水過程受力特征。從圖5可以看出，回轉(zhuǎn)體頭部受到的作用力在入水瞬間急劇上升，達(dá)到峰值后迅速下降，最后以小幅度震蕩趨勢變化。可見在入水瞬間產(chǎn)生了巨大沖擊力，所以在設(shè)計回轉(zhuǎn)體模型時要重點考慮入水瞬間的結(jié)構(gòu)強度。下面將對2種結(jié)構(gòu)形式的回轉(zhuǎn)體進(jìn)行強度的分析計算。

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

在模擬中，首先進(jìn)行回轉(zhuǎn)體等厚度工況分析，分別計算回轉(zhuǎn)體厚度為2、3、4和5 mm的情況。計算中發(fā)現(xiàn)：只有回轉(zhuǎn)體厚度為5 mm時，該種結(jié)構(gòu)的回轉(zhuǎn)體是安全的；此種結(jié)構(gòu)的回轉(zhuǎn)體發(fā)生破壞時，破壞部位均發(fā)生在回轉(zhuǎn)體頭部和錐形斜面連接處。對回轉(zhuǎn)體厚度為2 mm的計算結(jié)果進(jìn)行分析，其中1.1 ms和3.0 ms時刻的應(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)象?？梢?，回轉(zhuǎn)體頭部和錐形斜面連接處應(yīng)予以加強，來承受入水瞬間高數(shù)量級的沖擊載荷。同時因為入水過程中回轉(zhuǎn)體頭部為主要的沖擊載荷作用面，所以對回轉(zhuǎn)體頭部進(jìn)行加強是十分必要的。

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

圖7 在t=3.0 ms時回轉(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é)中計算結(jié)果可知，回轉(zhuǎn)體頭部和錐形斜面連接處強度脆弱，因此對回轉(zhuǎn)體頭部進(jìn)行加強，考慮到工程安全系數(shù)要求，將頭部賦予8 mm的厚度，回轉(zhuǎn)體其余板厚度分別賦予2.0、2.5、3.0、3.5和4.0 mm。在初速度相同的條件下分別進(jìn)行計算，由計算結(jié)果可知：在頭部厚度為8 mm、后體區(qū)域厚度為2.0 mm時，回轉(zhuǎn)體無法滿足強度要求；在后體區(qū)域厚度為2.5 mm及以上時，結(jié)構(gòu)強度滿足要求。選取后體區(qū)域厚度為2.5 mm的情況加以分析?；剞D(zhuǎn)體入水過程中0.89 ms和2.00 ms時刻的應(yīng)力云圖分別如圖10和圖11所示，同時圖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(zhuǎn)體頭部邊緣的塑性應(yīng)變由0迅速到達(dá)0.14，并最終趨近0.16，回轉(zhuǎn)體頭部中心的塑性應(yīng)變有同樣的變化規(guī)律并最終達(dá)到0.25。這個塑性應(yīng)變小于材料的塑性失效應(yīng)變0.35，因此，材料雖然進(jìn)入塑性階段，但還未達(dá)到破壞，此時回轉(zhuǎn)體結(jié)構(gòu)滿足強度要求。

圖10 t=0.89 ms時刻回轉(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)變隨時間的變化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é) 論

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

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

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

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