平頭錐型回轉體高速入水結構強度數值分析*

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

（1.大連理工大學船舶工程學院遼寧省深海浮動結構工程實驗室，遼寧 大連 116024；2.大連理工大學工業裝備結構分析國家重點實驗室，遼寧 大連 116024；3.高新船舶與深海開發裝備協同創新中心，上海 200240）

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

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

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

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

1 數值計算模型設置

1.1 流固耦合ALE算法及狀態方程

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

在計算中，首先進行Lagrange時步計算，單元網格隨材料流動開始變形，然后進行ALE時步計算：

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

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

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

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

對于水介質，其壓力使用Grüneisen狀態方程描述，在LS-DYNA中所使用的關鍵字為*EOS_Grüneisen，水壓力方程表示為：

式中：pw為水壓力，cw為聲音在水中的傳播速度， μw為體積變化率，α為對Grüneisen系數 γ0的一階體積修正，S1、S2、S3分別為us-up曲線斜率無量綱系數，us為沖擊波速度，up為流體質點的速度，Ew0為材料初始內能。模擬中水狀態方程參數：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 計算模型及網格劃分

建立計算域模型如圖1所示，計算域包括空氣域和水域。空氣域尺寸為1.2 m×1.2 m×1.0 m，水域尺寸為1.2 m×1.2 m×1.5 m。水域和空氣域寬度為回轉體最大直徑的10倍。在模擬中采用了2種回轉體,其結構形式分別為：(1)全回轉體等厚度；(b)回轉體頭部部分厚度為8 mm，后體區域賦予一種厚度。回轉體模型整體為平頭錐形結構，最大直徑為120 mm，頭部直徑為100 mm，回轉體長度為0.508 m，其首端傾斜斜面與回轉體頭部平面所呈角度為104.7°，回轉體外形如圖2所示。

圖1 計算域Fig.1 Computational domain

圖2 回轉體幾何模型Fig.2 The geometrical model for a revolution body

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

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

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

2 計算結果與分析

2.1 數值方法驗證

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

圖3 回轉體頭部中心區域壓強曲線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為厚度相關系數，d為回轉體厚度(m)；kE為彈性模量相關系數，E為材料彈性模量；kv為速度相關系數，v0為回轉體入水的初速度(m/s)。由式(4)可以求得在5 mm板厚下的回轉體入水沖擊載荷壓強峰值為173 MPa，和本文求得的數值解比較接近，但數值模擬結果略大于公式求得的抨擊壓強峰值。其原因是在數值模擬中未考慮空氣墊的作用，即高速時的空氣壓縮效應，這也在一定程度上解釋了數值模擬所得壓強數值偏大的現象，空氣墊效應是數值模擬過程中較難考慮的一個因素，在數值模擬方面還有很大的研究空間，亟待將來的研究。

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

假定回轉體在入水過程中空泡內外壓差保持不變，同時空化數為非定值，有：

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

柱形平頭回轉體阻力因數和空化數關系式：

忽略入水過程中重力的影響，回轉體有以下運動方程：

式中：mp為回轉體的質量，A0為回轉體頭部面積。

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

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

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

圖5 流-構耦合力曲線Fig.5 Fluid-structure interaction force varying with time

在入水數值驗證工作的同時，圖5給出了當前回轉體結構形式下入水過程中回轉體頭部表面作用力曲線圖，可以進一步分析回轉體入水過程受力特征。從圖5可以看出，回轉體頭部受到的作用力在入水瞬間急劇上升，達到峰值后迅速下降，最后以小幅度震蕩趨勢變化。可見在入水瞬間產生了巨大沖擊力，所以在設計回轉體模型時要重點考慮入水瞬間的結構強度。下面將對2種結構形式的回轉體進行強度的分析計算。

2.2 全回轉體等厚度工況模擬

在模擬中，首先進行回轉體等厚度工況分析，分別計算回轉體厚度為2、3、4和5 mm的情況。計算中發現：只有回轉體厚度為5 mm時，該種結構的回轉體是安全的；此種結構的回轉體發生破壞時，破壞部位均發生在回轉體頭部和錐形斜面連接處。對回轉體厚度為2 mm的計算結果進行分析，其中1.1 ms和3.0 ms時刻的應力云圖分別如圖6和圖7所示。由圖6～7可以清晰地看出，破壞部位為頭部和錐形斜面連接處，回轉體頭部中心處未發生破壞現象。提取回轉體頭部和錐形斜面連接處的塑性應變和單元應力曲線，分別如圖8和圖9所示。由這2條曲線可以得到：在回轉體入水過程中，等效應力達到了材料的屈服應力355 MPa，塑性應變也達到了材料的塑性失效應變0.35，因此回轉體出現破壞現象。可見，回轉體頭部和錐形斜面連接處應予以加強，來承受入水瞬間高數量級的沖擊載荷。同時因為入水過程中回轉體頭部為主要的沖擊載荷作用面，所以對回轉體頭部進行加強是十分必要的。

圖6 在t=1.1 ms時回轉體的應力分布Fig.6 Stress distribution in the revolution body at t=1.1 ms

圖7 在t=3.0 ms時回轉體的應力分布Fig.7 Stress distribution in the revolution body at t=3.0 ms

圖8 回轉體頭部邊緣單元應變曲線Fig.8 Strain-time curve of the edge element for the head of the revolution body

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

2.3 回轉體頭部壁厚的影響

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

圖10 t=0.89 ms時刻回轉體的應力分布Fig.10 Stress distribution in the revolution body at t=0.89 ms

圖11 t=2.00 ms回轉體應力圖Fig.11 Stress distribution in the revolution body at t=2.00 ms

圖12 回轉體頭部邊緣及中心單元應變隨時間的變化Fig.12 Strain-time curves of the edge and central elements for the head of the revolution body

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

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

3 結 論

采用有限元方法對平頭錐頭回轉體在高速入水過程中的結構強度進行了數值分析，得到的主要結論如下：

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

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

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