基于虛擬介質方法數值模擬應力波穿透水艙

史汝超， 孫曉旺

(1.江蘇海洋大學 理學院，江蘇 連云港 222005；2.連云港中復連眾復合材料集團有限公司，江蘇 連云港 222023； 3. 南京理工大學 機械工程學院，南京 210094)

水艙是艦艇的常見裝置，可分為壓載水艙、補重水艙、間歇水艙等多種專用水艙。應力波穿透水艙是一個典型的流固耦合問題，常見于艦艇的防護設計和抗爆抗沖擊研究[1-3]。目前，數值模擬該類問題較常使用的方法是SPH(smoothed particle hydrodynamics)[4]和ALE(arbitrary-Lagrangian-Eulerian)方法[5]。這兩種方法適用性很廣，因此，大部分主流商業軟件，如ANSYS、Abaqus等，都是基于這兩種方法開發的。然而，如果想非常精確地求解并刻畫固體中的應力波，很多時候需要使用數值精度較高的有限差分法，選擇特定的差分格式，并設以特定的時間步長和空間網格。目前，完全使用有限差分法處理流固耦合問題的計算方法是虛擬介質方法[6-7]。

虛擬介質方法可進一步細分為OGFM法則(original ghost fluid method)[8]、MGFM法則(modified ghost fluid method)[9]和RGFM法則(real ghost fluid method)[10]。收斂性較好的是MGFM法則和RGFM法則，均已初步推廣至水下沖擊波沖擊流固界面問題，可獲得很好的模擬效果[11-13]。本文進一步以MGFM法則為基礎，將虛擬介質方法推廣至固體應力波在流固界面處的折射問題，并以此為基礎，數值模擬應力波穿透水艙。

1 虛擬介質方法

1.1 流體控制方程

二維軸對稱歐拉方程可寫為

(1)

其中，

(2)

(3)

(4)

(5)

式中:ρ為密度；p為壓力；u，v，w為x，y，z方向的速度;r為任一點至y軸的距離；單位質量總能E=ρ(u2+v2+w2)/2+ρe。對于一維問題，式(1)可簡化為

(6)

1.2 固體控制方程

取壓應力為負值、拉應力為正值。固體的二維軸對稱彈性動力學方程可寫為

(7)

其中，

(8)

(9)

(10)

(11)

(12)

1.3 近似黎曼解算器

流固界面在數學上屬于黎曼問題。圖1 給出了求解黎曼問題的特征線法示意圖。在流固界面處沿特征線差分，可得到界面處的差分近似解。圖中固體設定在界面左側，故取右特征線方程

dσ-ρcdu=0

(13)

式中，c為聲速。流體在界面右側，故取左特征線方程

dp-ρcdu=0

(14)

式(11)、式(12)分別向后和向前差分

(15)

式中：下標“I”為界面處的參數；上標“n”，“n+1”為時間步。

在流固界面處：當固體處于受壓狀態，界面受力方向從流體指向固體，pI=-σI；當固體處于受拉狀態，界面受力方向從固體指向流體，pI=σI。

水的密度可由等熵狀態方程確定

(16)

式中：N為水的絕熱指數；ρ0=1.0×103kg/m3；A=1.0×103Pa；B=3.31×108Pa。

在流固界面變形較小的情況下，求解固體一側的網格坐標

(17)

圖1 特征線法求解黎曼問題

1.4 MGFM賦值法則

(18)

(19)

式中，dis為點與面之間的最短距離。

圖2 MGFM賦值法則

在數值模擬中，由于界面的移動，流體虛擬節點有時候會轉變為真實節點。根據等熵修正思想，為避免數值振蕩，用相鄰真實節點的狀態代替虛擬節點的狀態(見圖3)。

圖3 賦值修正步

1.5 數值模擬步驟

選取無量綱參考態：壓力pref=0.1 MPa,密度ρref=1.0×103kg / m3，長度lref=1.0 m，速度uref=10 m/s，時間tref=0.1 s。用Zwas格式[15]離散固體控制方程；用五階WENO格式[16]離散流體控制方程；用三階Runge-Kutta方法離散時間項。使用歐拉網格離散流體區域，并求解網格節點的狀態值；使用歐拉網格求解固體區域并求解網格中心的值。數值模擬步驟如下：

步驟1用近似黎曼解算器成對求解界面兩側的流體網格節點和固體網格中心的狀態；

步驟2用1.4節中的賦值法則對虛擬區域賦值；

步驟3用當前時刻的狀態值求解下一時刻流體網格節點(包括虛擬節點)和固體網格中心(包括虛擬網格中心)的值；

步驟4用式(17)求解界面位置，并根據新界面重新劃分真實區域和虛擬區域；

步驟5回到步驟1進行一下時間步的求解。

2 數值測試

圖4 固體側的拉伸波與流體側的稀疏波示意圖

將計算域劃分為400個網格，網格步長Δx=0.005。水的絕熱指數N=7.15；固體密度ρ=7.8、彈性模量E=2.06×106、泊松比ν=0.25。圖5表明數值結果和理論結果吻合很好且沒有數值振蕩。其中，理論結果的求解可參考文獻[17]。稀疏波的理論解在左右兩端各有一個折點；而稀疏波的數值解，由于數值誤差，兩端折點不明顯。表1給出了固體區域和流體區域的數值誤差，本文給出的賦值法則具有一階收斂精度。

圖5 拉伸波與稀疏波的數值結果(t=1.243 429×10-3)

表1 數值測試結果

3 一維數值模擬

水艙殼體中的應力波一般由水或者空氣中的壓縮波透射產生,如圖6所示。壓縮波透射產生固體應力波后，才能進一步在艙內產生另一道透射壓縮波。本節分別對水和空氣中的壓縮波穿透問題數值驗證。以下計算中，取網格步長Δx=0.005；空氣絕熱指數γ=1.38；水的絕熱指數N=7.15；固體密度ρ=7.8、彈性模量E=2.06×106、泊松比ν=0.25。

3.1 空氣中的壓縮波穿透水艙

圖6 空氣/水中的壓縮波穿透水艙

圖7 空氣中的壓縮波進入水艙殼體生成固體應力波(t=5.892 452×10-3)

圖8 固體應力波進入水艙生成水中壓縮波(3.1節)(t=8.081 076×10-3)

3.2 水中壓縮波穿透水艙

3.1節與3.2節中的數值結果均與理論解吻合一致。雖然虛擬介質方法中的流固界面在數值模擬時是連續的，但是在數學上仍然是間斷的。這使得在數值模擬時，流體與固體必須分開求解，從而導致數值結果與理論解在局部區域有一些較小的差異。但是，這些誤差是收斂的，并沒有隨著計算的進行而逐漸增長。

圖9 水中的壓縮波進入水艙殼體生成固體應力波(t=4.208 894×10-3)

圖10 固體應力波進入水艙生成水中壓縮波(3.2節)(t=6.734 230×10-3)

4 二維軸對稱數值模擬

圖11給出了一個水艙中心計算域的示意圖。計算域底面半徑0.5，高度1.0，殼板背面壓力分別取pw=100 kPa,200 kPa,400 kPa,800 kPa，取該中心區域的中軸面(虛線區域)建立二維軸對稱坐標系，網格步長為Δx=Δy=0.005。水的絕熱指數N=7.15。殼板無量綱參數為：密度ρ=7.85、彈性模量E=2.06×106、泊松比ν=0.28、屈服強度4 000、抗拉強度4 800。

在底面中心A點按圖12加載應力，加載范圍為直徑0.01的圓形區域。圖13(a)描繪了應力波在水艙殼體內部的傳播。如圖13(b)所示，在t=1.122×10-3時，應力波波陣面已經到達水艙殼體背面。圖14(a)、圖14(b)分別給出了由虛擬介質方法和ALE方法得到的背面中心B點的動態響應。由于數值耗散，背面中心B點的應力大約在t=7.0×10-4時開始增加。但波陣面大約于t=1.0×10-3時，到達背面中心，此時，應力顯著增加。ALE算法在流固界面處采用弱耦合形式，只使用迎風面(殼體側)的狀態參數，界面應力在峰值之前上升略早，峰值時間略晚；虛擬介質方法在流固界面處采用強耦合形式，即同時使用界面兩側的狀態參數，界面應力在峰值之前上升略晚，峰值時間略早。兩種方法得到的數值結果差異較小，應力變化趨勢相近。表2表明，當水艙內的環境壓力增加后，應力峰值也相應的增加，峰值時間近似保持不變；同時也驗證了虛擬介質方法與ALE方法所得到的應力峰值和峰值時間相近。圖15表明水艙壓力越大，透射波強度越大；同時，由于4個算例(pw=100 kPa,200 kPa,400 kPa,800 kPa)均使用相同的應力加載，殼體背面的應力波折射規律相近。

圖11 水艙中心計算域示意圖

圖12 應力加載曲線

圖13 固體中的應力波傳播(固體壓強ps=∑|σi|/3)

圖14 水艙殼體背面動態響應

表2 峰值數值結果

圖15 應力波在殼板背面的折射 (t=2.486×10-3，固體壓強ps=∑|σi|/3)

5 結 論

本文將虛擬介質方法應用到應力波在流固界面處的折射問題，并給出了對應的計算法則；數值測試表明，提出的計算法則具有一階精度。本文對水艙殼體中的應力波穿透水艙的過程進行了數值模擬和理論分析，一維數值結果與理論解能很好地吻合。多維數值模擬驗證了虛擬介質方法與ALE方法得到數值結果相近。對應力波穿透水艙的初步研究表明，水艙加壓后，透射波的強度也相應地增加；在不改變初始應力加載的情況下，應力波在水艙殼體背面的折射規律也對應的近似不變。