基于SPH-FEM 耦合方法的回轉(zhuǎn)體高速入水?dāng)?shù)值研究

高英杰，剛旭皓

(1.中國(guó)船級(jí)社青島分社，山東 青島 266071；2.中國(guó)船舶科學(xué)研究中心，江蘇 無(wú)錫 214082)

0 引 言

入水沖擊問(wèn)題具有廣泛的工程應(yīng)用背景，如船舶航行受浪砰擊，海洋工程結(jié)構(gòu)物下水安裝，救生艇自由降落釋放，導(dǎo)彈、魚(yú)雷入水等。隨著海洋工程計(jì)算和武器裝備的發(fā)展，入水沖擊問(wèn)題的研究受到持續(xù)關(guān)注。入水沖擊作為一個(gè)復(fù)雜的物理現(xiàn)象，伴隨著流固耦合作用，而高速入水大大加劇了流固耦合的復(fù)雜程度，同時(shí)結(jié)構(gòu)物撞水瞬間受到的砰擊載荷可能導(dǎo)致結(jié)構(gòu)運(yùn)動(dòng)變化和結(jié)構(gòu)變形失效。因此有必要考慮流固耦合作用，分析入水結(jié)構(gòu)的動(dòng)態(tài)響應(yīng)，為結(jié)構(gòu)物高速入水設(shè)計(jì)提供參考。

對(duì)于入水沖擊問(wèn)題，國(guó)內(nèi)外學(xué)者做了大量研究。1929 年，Von Karman首先引入附加質(zhì)量的概念，推導(dǎo)計(jì)算出二維楔形體入水沖擊載荷，奠定了入水沖擊理論計(jì)算的基礎(chǔ)。Wagner在1932 年，針對(duì)楔形體入水沖擊問(wèn)題提出了小傾角平板近似理論，通過(guò)求解伯努利方程計(jì)算出二楔形體入水過(guò)程中的壓力分布。后續(xù)學(xué)者在二者的研究基礎(chǔ)上，通過(guò)試驗(yàn)、數(shù)值計(jì)算等方式對(duì)入水沖擊問(wèn)題開(kāi)展了廣泛研究。Worthington采用高速攝影技術(shù)，最早記錄了小球入水的流場(chǎng)演變情況，分析了空泡的變化規(guī)律。MAY開(kāi)展了不同形狀參數(shù)和入水速度的試驗(yàn)，研究了結(jié)構(gòu)物入水過(guò)程中的載荷特征。顧建農(nóng)開(kāi)展了彈體高速入水試驗(yàn)，分析了形狀參數(shù)和入水初速度對(duì)彈道特征和運(yùn)動(dòng)特性的影響。張偉通過(guò)試驗(yàn)分析了不同頭型彈體入水的彈道穩(wěn)定性，得到了彈體入水速度衰減公式。

數(shù)值計(jì)算成本較低且不受場(chǎng)地、設(shè)備等現(xiàn)實(shí)條件限制，可較為準(zhǔn)確地預(yù)報(bào)入水過(guò)程中結(jié)構(gòu)和流場(chǎng)的變化情況，逐漸成為研究入水沖擊問(wèn)題的關(guān)鍵方法之一。陳震開(kāi)展了船舶在波浪中出入水的數(shù)值模擬，分析了入水運(yùn)動(dòng)參數(shù)對(duì)抨擊壓力的影響。何春濤采用有限體積法，建立了結(jié)構(gòu)物垂直入水的計(jì)算模型，對(duì)比試驗(yàn)結(jié)果分析了入水的運(yùn)動(dòng)軌跡和空泡閉合規(guī)律。馬慶鵬采用VOF 方法開(kāi)展了不同入水速度下回轉(zhuǎn)體的數(shù)值計(jì)算，研究了入水沖擊載荷變化規(guī)律以及侵水深度和空泡發(fā)展的關(guān)系。李強(qiáng)基于Ls-dyna 軟件計(jì)算了魚(yú)雷入水過(guò)程，分析了不同頭型參數(shù)對(duì)魚(yú)雷水下運(yùn)動(dòng)的影響。楊力采用ALE 方法，分析了平底結(jié)構(gòu)入水過(guò)程中質(zhì)量、剛度對(duì)砰擊壓力的影響。路麗?；赗ANS 方程采用動(dòng)網(wǎng)格技術(shù)，模擬了凹形體入水過(guò)程，分析了液體噴濺對(duì)角運(yùn)動(dòng)的影響。

目前針對(duì)結(jié)構(gòu)入水的研究多集中于剛體低速入水，主要關(guān)注入水運(yùn)動(dòng)軌跡和空泡演變規(guī)律，考慮流固耦合的高速入水研究較少。本文進(jìn)行回轉(zhuǎn)體高速入水流固耦合數(shù)值計(jì)算，分析回轉(zhuǎn)體高速入水的運(yùn)動(dòng)特性和結(jié)構(gòu)響應(yīng)，其結(jié)果可以為結(jié)構(gòu)物高速入水的預(yù)報(bào)提供參考。

1 流固耦合計(jì)算模型

采用光滑粒子流體動(dòng)力學(xué)（SPH）和有限元(FEM)耦合方法對(duì)高速入水過(guò)程進(jìn)行模擬，其中?；剞D(zhuǎn)體模型采用有限元Lagrange 網(wǎng)格，水域采用SPH 粒子。FEM 方法在計(jì)算結(jié)構(gòu)運(yùn)動(dòng)和動(dòng)力響應(yīng)方面,具有更高的精度和效率；SPH 方法中計(jì)算域使用粒子形式表示，無(wú)需網(wǎng)格劃分，在此計(jì)算介質(zhì)大變形，可以有效避免有限元方法中網(wǎng)格畸變等情況。FEM/SPH 耦合可發(fā)揮兩者優(yōu)勢(shì)，在計(jì)算入水沖擊等介質(zhì)大變形的動(dòng)力問(wèn)題上，具有較好的準(zhǔn)確性和較高計(jì)算效率。

1.1 SPH 方法基本理論

SPH 方法使用粒子分布表達(dá)計(jì)算域，采用核函數(shù)近似的方式將求解的偏微分方程轉(zhuǎn)換為積分形式，使用粒子近似的方式將核函數(shù)方程離散化，通過(guò)求解離散方程得到粒子場(chǎng)內(nèi)變量。

在SPH 方法中，給定任意一個(gè)函數(shù)，其積分表達(dá)形式如下：

式中：Ω為積分域；δ (-)為狄拉克函數(shù)。

由上式可知，任意連續(xù)函數(shù)()采用狄拉克函數(shù)可以表示為以下積分形式：

式中：(-,)為核函數(shù)或稱光滑函數(shù)，核函數(shù)的選取往往會(huì)影響計(jì)算精度和效率。

根據(jù)分布積分和散度定理，可以將上式進(jìn)行變換：

最后，將變換后的函數(shù)進(jìn)行離散化處理，得到域內(nèi)所有粒子疊加的離散形式：

式中：和為粒子編號(hào)；和為粒子的質(zhì)量和密度；為粒子總數(shù) 。

所以，在SPH 方法中對(duì)于基本控制方程可以表示為如下離散形式：

1.2 界面耦合模型

在入水過(guò)程中,實(shí)現(xiàn)結(jié)構(gòu)體和水中介質(zhì)耦合作用，需要在SPH-FEM 耦合方法中，將網(wǎng)格和粒子之間的界面進(jìn)行處理。結(jié)構(gòu)體網(wǎng)格和水粒子之間需要采用接觸算法進(jìn)行相應(yīng)的計(jì)算處理。

在接觸算法中，每一時(shí)間步內(nèi)計(jì)算有限元網(wǎng)格節(jié)點(diǎn)穿過(guò)SPH 粒子的距離，罰函數(shù)根據(jù)有限元網(wǎng)格和SPH 粒子間的穿透距離計(jì)算出兩者之間的接觸力，通過(guò)接觸力實(shí)現(xiàn)兩者之間的相互作用，接觸力將SPH 粒子拉出交界面，避免SPH 粒子穿越有限元網(wǎng)格，導(dǎo)致計(jì)算失效。

1.3 計(jì)算模型設(shè)置

在計(jì)算模型中，回轉(zhuǎn)體最大直徑 120 mm，頭部直徑 100 mm，回轉(zhuǎn)體長(zhǎng)度為0.5 m，結(jié)構(gòu)厚度為5 mm，其首端傾斜斜面與回轉(zhuǎn)體頭部平面所呈角度為 93.18°，入水速度為100 m/s，采用有限元網(wǎng)格進(jìn)行劃分。水域尺寸為 1.2 m×1.2 m×3.6 m，采用SPH 粒子表示，回轉(zhuǎn)體外形及流域如圖2 所示。

圖2 計(jì)算域設(shè)置Fig.2 Computational domain settings

SPH 粒子數(shù)為1 536 000，使用 Gruneisen 狀態(tài)方程控制，回轉(zhuǎn)體采用4 節(jié)點(diǎn)殼單元，網(wǎng)格數(shù)量為520 000，材料模型考慮彈性和塑形，材質(zhì)為45 號(hào)鋼材，材料性能參數(shù)見(jiàn)表1。

表1 材料模型參數(shù)Tab.1 Material parameters

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

2.1 運(yùn)動(dòng)特性分析

首先進(jìn)行回轉(zhuǎn)體入水速度變化分析并與理論公式對(duì)比，驗(yàn)證計(jì)算方法的準(zhǔn)確性。假設(shè)物體在入水過(guò)程中空泡內(nèi)部和外部壓力差保持一致，在空化數(shù)不為定值的情況下，存在如下規(guī)律：

式中：Δ為空泡內(nèi)外壓力差；C為氣流壓力降系數(shù)；σ為初始空化數(shù)；σ=0.006-0.018 ；ρ為空氣密度；ρ為水密度；v為入水初速度；v為入水后速度。

阻力系數(shù)與空化數(shù)關(guān)系：

忽略重力的情況下，入水運(yùn)動(dòng)方程為：

式中：m為物體的質(zhì)量；A為觸水面積。

通過(guò)式(9)和式(10)，可得到速度衰減公式：

式中：衰減系數(shù)k=ρ A C/ 2m。

圖3 為衰減公式與數(shù)值計(jì)算得到的速度曲線?？梢钥闯?，觸水瞬間回轉(zhuǎn)體高速撞擊水面受到極大的瞬間沖擊載荷，速度衰減劇烈，隨著入水進(jìn)程加深，速度衰減減弱，最后近似線性衰減。計(jì)算得到的速度曲線與公式曲線具有好的一致性，驗(yàn)證了SPH-FEM 耦合計(jì)算模型在計(jì)算高速入水方面的有效性。

圖3 回轉(zhuǎn)體速度衰減曲線對(duì)比圖Fig.3 Curve of velocity attenuation comparison

2.2 動(dòng)力特性及結(jié)構(gòu)動(dòng)力響應(yīng)分析

王珂擬合了回轉(zhuǎn)體入水抨擊壓強(qiáng)峰值，入水過(guò)程中回轉(zhuǎn)體頭部中心處峰值壓強(qiáng)可以表示為：

式中：k=0.004 23×d+0.003 758×d+0.42，ln+15，k=+5；ρ為水的密度；為厚度；為材料彈性模量；為入水的初速度。

由式(11)可以計(jì)算出回轉(zhuǎn)體入水頭部中心錯(cuò)誤壓力峰值為1.73×10Pa，采用SPH-FEM 耦合計(jì)算得到的回轉(zhuǎn)體頭部中心處壓力情況如圖4 所示，其峰值為1.92×10Pa。可以看出，數(shù)值計(jì)算得到的回轉(zhuǎn)體頭部中心處壓力峰值與文獻(xiàn)公式計(jì)算結(jié)果接近，數(shù)值計(jì)算結(jié)果稍大。在數(shù)值計(jì)算過(guò)程中未充分考慮高速入水過(guò)程中壓縮空氣的影響，導(dǎo)致結(jié)果略微偏大。高速入水是一個(gè)涉及氣-液-固三相耦合的瞬時(shí)沖擊問(wèn)題，空氣高速壓縮帶來(lái)的空氣墊效應(yīng)一直作為數(shù)值計(jì)算中的難點(diǎn)，在未來(lái)具有廣闊的研究空間。

圖4 回轉(zhuǎn)體中心單元壓力圖Fig.4 Pressure intensity curve in the central head region of the revolution body

圖5 給出了入水過(guò)程中回轉(zhuǎn)體垂直方向受力情況。從圖中可以看出，撞水瞬間回轉(zhuǎn)體頭部受到的沖擊力急劇增加，在毫秒量級(jí)時(shí)間達(dá)到峰值，峰值為1.37×10N，頭部觸水后，峰值迅速下降，隨著回轉(zhuǎn)體完全進(jìn)入水面，沖擊力呈小幅度震蕩衰減，逐漸趨于穩(wěn)定。

圖5 回轉(zhuǎn)體沿軸向受力圖Fig.5 Curve of the axial force of the revolution body

高速入水過(guò)程中，撞水初期會(huì)產(chǎn)生瞬時(shí)巨大沖擊力，可能導(dǎo)致結(jié)構(gòu)變形甚至失效，因此有必要研究入水沖擊的結(jié)構(gòu)強(qiáng)度問(wèn)題。圖6 為回轉(zhuǎn)體頭部中心單元的等效應(yīng)力曲線，圖7 給出了回轉(zhuǎn)體頭部中心單元的塑形應(yīng)變曲線。本次計(jì)算考慮流固耦合效應(yīng)，物體撞水后會(huì)產(chǎn)生變形緩沖抨擊，變形呈反復(fù)狀態(tài)。在應(yīng)力曲線上表現(xiàn)為：出現(xiàn)瞬時(shí)峰值后，應(yīng)力衰減呈現(xiàn)波動(dòng)震蕩形式，隨著入水進(jìn)程加深，震蕩的頻率和幅度減小，伴隨著整個(gè)入水進(jìn)程。在100 m 垂直入水情況下，回轉(zhuǎn)體頭部中心處應(yīng)力峰值迅速達(dá)到材料的屈服應(yīng)力 355 MPa，材料發(fā)生塑形應(yīng)變，塑性應(yīng)變值0.146，雖然未出現(xiàn)破壞現(xiàn)象，但是在設(shè)計(jì)過(guò)程中有必要加強(qiáng)結(jié)構(gòu)物觸水位置處的結(jié)構(gòu)強(qiáng)度，避免結(jié)構(gòu)失效。

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

圖7 回轉(zhuǎn)體頭部單元塑形應(yīng)變曲線Fig.7 Plastic strain curve of the central element of the revolution body head

2.3 入水空泡特性分析

入水空泡形態(tài)是研究空泡問(wèn)題的關(guān)鍵，研究入水空泡動(dòng)力學(xué)問(wèn)題在基于一定的假設(shè)條件下可以推導(dǎo)出在入水過(guò)程中的侵徹距離與空泡半徑的函數(shù)關(guān)系。此時(shí)假定入水過(guò)程中的能量流失轉(zhuǎn)化為回轉(zhuǎn)體周圍空泡區(qū)域的流體動(dòng)能和空泡介質(zhì)勢(shì)能。利用分布點(diǎn)源理論推導(dǎo)出回轉(zhuǎn)體周圍流體動(dòng)能為：

式中：為空泡半徑；為無(wú)量綱系數(shù)，在忽略重力作用下，對(duì)有限長(zhǎng)度的平頭回轉(zhuǎn)體，其空泡內(nèi)勢(shì)能E可以表示為：

在忽略入水初始階段入水速度的影響下，同時(shí)根據(jù)能量守恒原理有：

通過(guò)式(12)～式(16)，同時(shí)應(yīng)用入水速度關(guān)系式v=(-)/(-)，最終可以求得在入水過(guò)程中的空泡半徑隨入水侵徹距離的函數(shù)關(guān)系：

式中：為時(shí)刻的位移值，為回轉(zhuǎn)體頭部的空泡半徑長(zhǎng)度?？梢钥闯?，回轉(zhuǎn)體入水過(guò)程中的空泡尺寸與回轉(zhuǎn)體的形狀、經(jīng)驗(yàn)系數(shù)、阻力系數(shù)以及空化數(shù)等有關(guān)，對(duì)高速入水問(wèn)題學(xué)者提出經(jīng)驗(yàn)系數(shù)取值一般小于3。

在得到圖8 各個(gè)時(shí)刻回轉(zhuǎn)體入水空泡體積分?jǐn)?shù)的前提下，提取回轉(zhuǎn)體入水的最終狀態(tài)即0.005 s 時(shí)刻的空泡氣-液交界面輪廓，同時(shí)利用式(17)的空泡半徑與侵徹距離公式關(guān)系，取變量數(shù)值=2。對(duì)比數(shù)值模擬和理論計(jì)算得到的空泡尺寸輪廓如圖9 所示。從圖中可以看出2 條曲線表征的空泡尺寸吻合較好，證明了數(shù)值模擬中計(jì)算空泡的準(zhǔn)確性。

圖8 回轉(zhuǎn)體入水開(kāi)空泡過(guò)程Fig.8 Cavity feature during the water-entry process

圖9 回轉(zhuǎn)體泡輪廓與理論模型對(duì)比Fig.9 Comparison of the cavity model with the simulation

另一方面，理論公式為剛體入水計(jì)算結(jié)果，而本次數(shù)值計(jì)算考慮到材料的彈塑性，更加接近于實(shí)際情況。由于回轉(zhuǎn)體在入水過(guò)程中存在變形，可以看出2 條空泡輪廓曲線存在一定的差別。彈性體入水的首端空泡尺寸相比剛體結(jié)果較小，回轉(zhuǎn)體頭部為入水過(guò)程中變形最明顯的區(qū)域，在入水過(guò)程中，頭部的變形使流體進(jìn)入內(nèi)凹區(qū)域，導(dǎo)致彈性體的入水空泡形狀首部區(qū)域不再是平面，使其首部的直徑變小，隨著入水進(jìn)程加深，后體區(qū)域的尺寸差異隨空泡直徑的增大逐漸減小。

2.4 入水速度對(duì)入水進(jìn)程的影響

為進(jìn)一步探究入水速度的影響，開(kāi)展回轉(zhuǎn)體以60 m/s，80 m/s，100 m/s 速度條件下的計(jì)算。表2 給出了不同入水速度下回轉(zhuǎn)體頭部中心處壓力峰值，圖10為不同入水初速度下的歸一化速度曲線?？梢钥闯?，回轉(zhuǎn)體入水初期受到極大的沖擊力，速度衰減劇烈，之后沖擊力迅速降低，回轉(zhuǎn)體速度衰減程度放緩明顯，最后速度衰減表現(xiàn)為類似線性變化?；剞D(zhuǎn)體速度衰減情況受初速度影響，初速度越大速度衰減程度越大，隨著入水進(jìn)程加深，不同初速度下的回轉(zhuǎn)體速度衰減程度接近，趨于穩(wěn)定。

圖1 有限元網(wǎng)格和SPH 粒子耦合界面Fig.1 Coupling interface of FEM and SPH

表2 回轉(zhuǎn)體頭部中心壓力峰值Tab.2 Results of pressure peak of cylinders head central point

圖10 歸一化速度衰減曲線Fig.10 Curve of normalizes velocity attenuation

回轉(zhuǎn)體入水阻力系數(shù)計(jì)算公式為：

其中：為阻力值；ρ為水密度；為入水速度；為特征面積。

圖11 為初速下回轉(zhuǎn)體入水的阻力系數(shù)曲線?？梢钥闯觯诓煌胨俣认禄剞D(zhuǎn)體阻力系數(shù)變化情況趨于一致，阻力系數(shù)峰值接近，入水速度越大，峰值持續(xù)時(shí)間越短，隨著入水進(jìn)程發(fā)展，流動(dòng)穩(wěn)定，3 種情況的阻力系數(shù)值保持接近。由此可見(jiàn)，入水初速度對(duì)回轉(zhuǎn)體阻力系數(shù)的影響不明顯。

圖11 阻力系數(shù)曲線Fig.11 Curve of Cd with different velocity

3 結(jié) 語(yǔ)

本文采用SPH-FEM 耦合方法對(duì)回轉(zhuǎn)體高速入水過(guò)程進(jìn)行數(shù)值分析，得到以下結(jié)論：1）通過(guò)與速度衰減和回轉(zhuǎn)頭部壓力峰值的理論公式對(duì)比，表明SPHFEM 耦合方法可以較為準(zhǔn)確地預(yù)報(bào)回轉(zhuǎn)體高速入水過(guò)程中的運(yùn)動(dòng)及動(dòng)力特征，可以較好地解決數(shù)值計(jì)算過(guò)程中由于流體介質(zhì)不連續(xù)變形帶來(lái)的計(jì)算困難等問(wèn)題。

2）回轉(zhuǎn)體觸水瞬間產(chǎn)生應(yīng)力峰值，出現(xiàn)時(shí)間極短，受變形影響，隨后應(yīng)力衰減呈現(xiàn)波動(dòng)震蕩形式，其頭部中心處應(yīng)力峰值超過(guò)了材料的屈服應(yīng)力，材料發(fā)生塑形應(yīng)變，在設(shè)計(jì)上可考慮對(duì)頭部結(jié)構(gòu)進(jìn)行加強(qiáng)。

3）數(shù)值計(jì)算中得到的入水空泡輪廓尺寸數(shù)據(jù)與剛性回轉(zhuǎn)體的理論空泡尺寸結(jié)果總體吻合良好，由于彈性變形的影響，在空泡輪廓頭部存在一定的差異，后體區(qū)域的尺寸差異隨空泡直徑的增大逐漸減小。

4）回轉(zhuǎn)體速度衰減受入水初速度影響，初速度越大速度衰減程度越大，隨著入水進(jìn)程加深，不同初速度下的回轉(zhuǎn)體速度衰減程度接近，入水速度對(duì)阻力系數(shù)變化影響不大。