單流體空化模型在水下爆炸誘導(dǎo)空化問題中的對比分析

金澤宇, 殷彩玉, 孔祥韶

(1. 武漢理工大學 綠色智能江海直達船舶與郵輪游艇研究中心, 武漢 430063; 2. 華中科技大學 船舶與海洋工程學院, 武漢 430074)

水下爆炸誘導(dǎo)空化主要包括自由面附近的片空化和結(jié)構(gòu)表面附近的局部空化。在水下爆炸后,沖擊波在水中傳播并與自由面或水下結(jié)構(gòu)相互作用,自由面或結(jié)構(gòu)的運動與變形使水中產(chǎn)生稀疏波。當入射波、反射波、稀疏波的壓力與靜水壓在某點的疊加壓力值小于飽和蒸汽壓,流體就會發(fā)生空化。由于片空化問題中水的重力作用和局部空化問題中結(jié)構(gòu)的減速作用,空化會潰滅,產(chǎn)生以沖擊波或者流體動量的形式加載到艦船或水下航行器的作用力。典型片空化和局部空化示意圖,如圖1所示。空化潰滅載荷是一種源自水下爆炸的典型的破壞載荷。

水下爆炸載荷與不同邊界相互作用產(chǎn)生的空化得到廣泛研究。Kennard[1]在雙線彈性流體中理論描述了空化。他的研究表明,當流體壓力小于空化極限時,會產(chǎn)生空化,空化的邊界或者是以超聲速運動的破水波,或者形成保持靜止的自由面,或者退化為以亞聲速運動的閉合波。根據(jù)雙線彈性流體模型,研究人員對水下沖擊波與懸浮在自由面的平板相互作用產(chǎn)生的空化效應(yīng)[2],水下沖擊波與線彈簧支撐平板[3]、厚彈性板[4]、夾芯結(jié)構(gòu)[5]相互作用產(chǎn)生的空化效應(yīng)建立了理論模型。Schiffer等[6]和Feng等[7]搭建了透明水激波管試驗?zāi)Ｐ?開展了系列試驗研究,記錄了近壁面空化的初生、擴展及潰滅過程。

數(shù)值模擬是研究空化的重要手段。典型的空化包括單流體空化模型和雙流體空化模型。雙流體模型假設(shè)汽相和液相在流體域中同時存在,由獨立的微分方程控制,其控制方程數(shù)量是單流體模型的兩倍。雙流體空化模型能夠考慮質(zhì)量/能量轉(zhuǎn)化、熱傳遞、表面張力等。但兩相流的一些初始條件很難直接獲取。典型的雙流體模型包括4方程模型[8-9]、6方程模型[10]和5方程模型[11]等。單流體模型假設(shè)汽液混合物宏觀表現(xiàn)為單一流體,遵循統(tǒng)一的狀態(tài)方程,因此只需要求解一套守恒方程,其難點是構(gòu)建合理的混合物狀態(tài)方程。典型的單流體模型包括截斷模型,修正的Schmidt模型[12],等熵空化流模型[13]等。單流體模型由于構(gòu)造簡單,在水下爆炸問題中得到廣泛應(yīng)用[14-16]。但在應(yīng)用單流體空化模型在水下爆炸問題時,不同模型究竟有多大差異仍存疑。

基于上述問題,本文建立經(jīng)典水下爆炸問題的數(shù)值模型,對比單流體空化模型的截斷模型,修正的Schmidt模型,等熵空化流模型的計算結(jié)果,得到不同空化模型的差異,為應(yīng)用單流體模型求解水下爆炸問題提供依據(jù)。

1 空化模型

1.1 常態(tài)水狀態(tài)方程

水下爆炸問題水通常采用Tait狀態(tài)方程描述

(1)

式中:ρ為流體的密度;p為流體的壓力;參考密度ρ0= 1 000 kg/m3;常數(shù)A= 1.0 ×105Pa;常數(shù)B= 3.31 ×108Pa; 常數(shù)N= 7.15。Tait狀態(tài)方程在壓力小于20 000 atm都可以描述水的狀態(tài)。

在水的壓力較小時,可將Tait狀態(tài)方程線性化,可得到線性方程描述水

(2)

作為Tait狀態(tài)方程的簡化模型,線性方程具有求解簡單,計算效率高等優(yōu)點,可用于關(guān)于遠場沖擊波輸入的問題。其中ρ0、p0和c0為初始給定密度、壓力和相應(yīng)的聲速。此方程具有在水發(fā)生空化前聲速為常數(shù)的特點。

1.2 空化流狀態(tài)方程

截斷模型可以表示成

(3)

式中,飽和蒸汽壓psat可取2 068.5 Pa。

修正Schmidt模型可以表示成

p=

(4)

式中:飽和蒸汽壓psat=2 068.5 Pa; 正壓參數(shù)pε= 1×10-9Pa;α為孔隙分數(shù)?；旌狭黧w中:水蒸氣密度ρg= 0.08 kg/m3; 水蒸氣聲速ag=190 m/s2; 液態(tài)水密度ρl= 999.958 6 kg/m3; 液態(tài)水聲速al=1 538.2 m/s2。pgl可由下式確定

(5)

等熵空化流模型可以表示成

(6)

(7)

(8)

將不同空化模型運用MATLAB編程,并畫出壓力小于飽和蒸汽壓的壓力密度曲線,如圖2所示。當計算壓力小于飽和蒸汽壓時,截斷模型(Cut Off)的壓力保持在飽和蒸汽壓,而修正的Schmidt模型(M-Schmidt)和等熵空化流模型(Isentropic)隨著密度降低壓力迅速下降。從圖2可知,M-Schmidt模型在較小的密度變化條件下,下降得更加迅速,而Isentropic模型表現(xiàn)出較為平滑的曲線。在壓力小于飽和蒸汽壓且大于修正Schmidt模型的截斷壓力pε時,相同密度條件下,運用M-Schmidt模型獲得的壓力小于Isentropic模型。

2 計算實例

本章開展水下沖擊波與浮在自由面平板相互作用,水下沖擊波與彈性支撐平板相互作用以及下沖擊波與水背襯平板相互作用模擬。歐拉方程采用3階龍格庫塔間斷伽遼金(RKDG)方法求解,其中空間離散采用基于間斷伽遼金(DG)方法,時間離散采用龍格庫塔(RK)方法,數(shù)值通量采用一階中心(FORCE)通量,而非線性斜率限制器選取用全變差有界(TVB)檢測的von Leer 全變差下降(TVD)限制器(在歐拉方程的特征量上使用),結(jié)構(gòu)的運動方程采用3階龍格庫塔法數(shù)值求解,流固耦合基于虛擬流體法,具體的求解過程可以參考Jin等,求解程序基于Fortran編譯器編制。

2.1 水下沖擊波與浮在自由面平板相互作用

水下沖擊波與浮在自由面平板相互作用問題,首先由Bleich等用特征線法求解該問題。由于重力引起水下落過程時間較長,直接研究自由面與沖擊波相互作用不利于分析,因此在自由面上懸浮一平板,可以使水中空化潰滅時間提前,其實質(zhì)特征與沖擊波與自由面相互作用類似。水下沖擊波與浮在自由面平板相互作用問題示意圖,如圖3所示。

由于該問題需考慮重力,因此需對水的控制進行修正?？紤]重力的水的控制方程為

Ut+F(U)x=S(U)

(9)

其中:

(10)

由于沖擊波壓力不大,線性化水的狀態(tài)方程。水下沖擊波與浮在自由面平板相互作用問題初始條件取ρ0= 989 kg/m3,c0= 1 451 m/s,p0= 0.101 MPa,g= 9.81 kg m/s2,pmax= 0.71 MPa,θ= 0.995 8 ms,m1= 143.405 kg/m2。計算域取[-10, 0.2]m,網(wǎng)格大小選取10 mm。

水的密度、速度、壓力賦值方程為

p(x)=pmaxe-(t-x/c0)/τ+p0+m1g-ρ0gx

u(x)=pmaxe-(t-x/c0)/θ/ρ0c0

(11)

2.1.1 飽和蒸汽壓取2 068.5 Pa時不同模型對比

當飽和蒸汽壓psat取2 068.5 Pa時,不同空化模型計算空化區(qū)域隨著時間變化曲線和平板速度時間歷程曲線的結(jié)果對比如圖4所示。從圖4(a)可知,空化2.4 ms以前,不同空化模型計算向遠處傳播的空化并無差別,而在2.4 ms以后,截斷模型在空化遠端開始向近端方向潰滅,而修正的Schmidt模型與等熵空化流模型的空化則繼續(xù)向遠處傳播一定距離后才開始潰滅。修正的Schmidt模型空化傳播距離最遠,其次為等熵空化流模型。修正的Schmidt模型空化潰滅的速度比等熵空化流模型更快,等熵空化流模型在遠端潰滅過程較為平緩。在近結(jié)構(gòu)端為空化潰滅邊界,逐漸向遠離結(jié)構(gòu)方向移動,修正的Schmidt模型和等熵空化流模型計算結(jié)果基本一致,而截斷模型潰滅速度較慢,與另外兩個模型有一定差別。截斷模型空化完全閉合發(fā)生的深度小于另外兩個模型,發(fā)生時刻也晚于另外兩個模型。修正的Schmidt模型和等熵空化流模型完全閉合深度與發(fā)生時刻基本一致。從圖4(b)可知,不同空化模型計算結(jié)果對結(jié)構(gòu)速度曲線影響很小,由此可以推斷不同空化模型計算的空化閉合載荷的差別對結(jié)構(gòu)響應(yīng)影響不大。

截斷模型、修正的Schmidt模型和等熵空化流模型在飽和蒸汽壓取2 068.5 Pa,t= 10.4 ms時刻壓力場和速度場計算結(jié)果對比如圖5所示。根據(jù)圖4(a)的結(jié)果,在t= 10.4 ms時,空化完全潰滅。從圖5可以看,截斷模型計算的場結(jié)果與修正的Schmidt模型、等熵空化流模型相比有顯著差異。空化潰滅時刻晚于另外兩種模型。而等熵空化流模型與修正的Schmidt模型在靠近結(jié)構(gòu)端也出現(xiàn)細微差異。修正的Schmidt模型潰滅載荷向結(jié)構(gòu)方向傳播的距離比等熵空化流模型潰滅載荷向結(jié)構(gòu)方向傳播的距離遠。而在潰滅點遠離結(jié)構(gòu)方向,兩種模型計算結(jié)果無顯著差異。

(a) 壓力場

2.1.2 飽和蒸汽壓取0.1 Pa時不同模型對比

Bleich和Sandler研究此問題時,將流體假設(shè)為雙線性水,當流體壓力小于零時用零替代。由于修正Schmidt模型和等熵空化流模型需要將飽和蒸汽壓設(shè)置為正數(shù),否則無法求解,本文取飽和蒸汽壓力0.1 Pa。飽和蒸汽壓的敏感性分析表明取飽和蒸汽壓力0.1 Pa接近飽和蒸汽壓取零的結(jié)果,可與Bleich和Sandler解析結(jié)果對比。

飽和蒸汽壓力取0.1 Pa時不同空化模型計算空化區(qū)域隨著時間變化曲線和板速度時間歷程曲線的結(jié)果對比,如圖6所示。從圖6(a)可知,截斷模型、修正的Schmidt模型和等熵空化流模型在飽和蒸汽壓力取0.1 Pa計算得到空化區(qū)域隨著時間變化結(jié)果無顯著差別。參考圖4可以推斷,不同模型空化區(qū)域結(jié)果與飽和蒸汽壓力取值有很大關(guān)聯(lián)。不同空化模型與解析結(jié)果存在一定差別,但差別不大。差別可能來源于初始條件、流體狀態(tài)方程無法完全一致引起。從圖6(b)可知,不同空化模型計算平板速度時間歷程曲線結(jié)果無顯著差別,同時與Bleich和Sandler解析結(jié)果也差別不大。與Bleich和Sandler解析結(jié)果的對比也可以印證空化模型的正確性。

(a) 空化區(qū)域演化

截斷模型、修正的Schmidt模型和等熵空化流模型在飽和蒸汽壓取0.1 Pa,t= 10.2 ms時刻壓力場和速度場計算結(jié)果對比如圖7所示。根據(jù)圖6(a)的結(jié)果,在t= 10.2 ms時,空化完全潰滅。從圖7可知,截斷模型、修正的Schmidt模型和等熵空化流模型計算得到的流場結(jié)果無顯著差別。對比圖5也可以推斷,不同模型空化區(qū)域結(jié)果與飽和蒸汽壓力取值有很大關(guān)聯(lián)。

(a) 壓力場

2.2 水下沖擊波與彈性支撐平板相互作用

圖8 水下沖擊波與彈性支撐平板相互作用模型示意圖

計算域取[-4,0.01]m,網(wǎng)格尺寸0.2 mm,m= 10 kg/m2,ψ= 2.5,κ= 0.2,pmax= 10 MPa, 狀態(tài)方程中取p0= 0.101 MPa。由于理論結(jié)果空化壓力取0,因此為了與理論結(jié)果對比,取飽和蒸汽壓psat= 0.01 Pa。截斷模型與Schiffer理論計算結(jié)果對比如圖9所示。從圖9可知,截斷模型計算空化區(qū)域隨著時間變化的無量綱曲線、彈性支撐板濕表面無量綱壓力時程曲線以及彈性支撐板無量綱速度時程曲線與Schiffer理論計算結(jié)果均較好吻合,說明計算模型的正確性。

(a) 空化區(qū)域演化

2.2.1 輸入峰值壓力10 MPa不同模型對比

圖9為與理論結(jié)果對比,取飽和蒸汽壓0.01 Pa。為更接近實際結(jié)果,下面研究取飽和蒸汽壓2 068.5 Pa,對比不同空化模型對空化載荷以及結(jié)構(gòu)響應(yīng)的影響。取輸入峰值壓力pmax= 10 MPa。對比截斷模型、修正的Schmidt模型和等熵空化流模型計算空化區(qū)域隨著時間變化無量綱曲線、彈性支撐板濕表面無量綱壓力時程曲線以及彈性支撐板無量綱速度時程曲線的結(jié)果,如圖10所示。計算結(jié)果表明在峰值壓力pmax= 10 MPa,飽和蒸汽壓取2 068.5 Pa時,截斷模型、修正的Schmidt模型和等熵空化流模型計算空化區(qū)域隨著時間變化、濕表面壓力、彈性支撐板速度結(jié)果并無顯著差別。

(a) 空化區(qū)域演化

(a) 壓力場

2.2.2 輸入峰值壓力0.3 MPa不同模型對比

當輸入峰壓pmax= 10 MPa時,不同空化模型計算結(jié)果區(qū)別較小,這可能由于流體壓力相比空化壓力較大引起的,下面取輸入峰壓pmax= 0.3 MPa,空化壓力取2 068.5Pa時,開展不同空化模型對比研究。

對比截斷模型、修正的Schmidt模型和等熵空化流模型計算空化區(qū)域隨著時間變化無量綱曲線、彈性支撐板濕表面無量綱壓力時程曲線以及彈性支撐板無量綱速度時程曲線的結(jié)果,如圖12所示。圖12(a)計算得到的空化區(qū)域隨著時間變化無量綱曲線相比圖10(a)空化區(qū)域顯著減小,這是由于入射壓力變小,導(dǎo)致結(jié)構(gòu)速度變小,從而引起稀疏波變小,這使得流體更難發(fā)生空化,因此空化區(qū)域變小。而隨著輸入峰壓降低,可以看到修正的Schmidt模型計算遠結(jié)構(gòu)端空化域的傳播比截斷模型和等熵空化模型計算結(jié)果快。而近結(jié)構(gòu)端空化潰滅邊界的傳播,不同模型計算結(jié)果差距不大。從圖12(b)和圖12(c)可知,即使輸入峰壓取0.3 MPa,不同模型計算濕表面壓力和板速度結(jié)果差別很小。

(a) 空化區(qū)域演化

(a) 壓力場

2.3 水下沖擊波與水背襯平板相互作用

水下沖擊波與水背襯平板相互作用問題,Schiffer等通過試驗方法開展過相關(guān)研究。通過該問題的分析,可以深化不同空化模型對由于水背襯結(jié)構(gòu)運動而產(chǎn)生的局部空化的影響的認識。水下沖擊波與水背襯平板相互作用模型示意圖如圖14所示。為了與Schiffer試驗結(jié)果對比,取流體密度ρw= 1 000 kg/m3, 流體聲速cw= 1 053 m/s,計算域[-4,2]m,網(wǎng)格尺寸0.5 mm,m= 111.1 kg/m2,ψ= 0.91,pmax= 15.1 MPa,θ= 0.096 ms, 狀態(tài)方程中取p0= 0.101 MPa。空化壓力取2 068.5Pa計算結(jié)果與Schiffer試驗結(jié)果對比,如圖15所示。從圖15可知,空化邊界和潰滅邊界的計算結(jié)果與試驗結(jié)果吻合較好。出現(xiàn)的誤差可能是來自觀測誤差以及平面波假設(shè)不完全(試驗中透明管子直徑為27 mm)。

圖14 水下沖擊波與水背襯平板相互作用模型示意圖

圖15 水下沖擊波與水背襯平板相互作用問題試驗和截斷

三種空化模型分別計算空化區(qū)域、面板前壓力、面板后壓力和面板速度隨時間變化的結(jié)果對比,如圖16所示。從圖16(a)可知,空化邊界在t= 0.85 ms以后不同模型出現(xiàn)了不同曲線,修正的Schmidt模型繼續(xù)以原來的傳播速度向遠處傳播,等熵空化流模型則以相對緩慢的速度繼續(xù)向遠處傳播,而截斷模型則向結(jié)構(gòu)端方向潰滅,近結(jié)構(gòu)端的潰滅邊界三種模型并無差別。這可能由于數(shù)值誤差引起的。由于間斷波在傳播過程中,通常在間斷波前存在較小的波動,當遠端壓力低到與間斷波前的波動量級接近時,會表現(xiàn)出較大數(shù)值誤差,而空化區(qū)域?qū)O小的波動敏感,從而影響圖16(a)中的曲線。圖16(a)中在遠端的空化區(qū)域可能會因為數(shù)值求解格式、數(shù)值通量、求解格式精度等因素不同而存在差別。但實際上這對結(jié)構(gòu)和流體整體響應(yīng)結(jié)果影響極小。從圖16(b)、圖16(c)和圖16(d)三種模型計算結(jié)構(gòu)響應(yīng)結(jié)果并無差別。

(a) 空化區(qū)域演化

水下沖擊波與水背襯平板相互作用問題在t= 1.5 ms時刻不同空化模型計算壓力場和速度場結(jié)果,如圖17所示。在t= 1.5 ms時刻,參考圖16(a),在沖擊波后,空化發(fā)生后,擴展到較大規(guī)模,且三種模型計算遠端的空化邊界具有很大不同?？栈瘍?nèi)流體速度仍朝向結(jié)構(gòu)方向。并且越靠近結(jié)構(gòu)空化區(qū)域的流體速度越大。從圖17可知,不同空化模型計算壓力場和速度場在空化區(qū)域附近無明顯差別。

(a) 壓力場

雖然圖16(a)計算的遠端的空化邊界存在差別,但在場計算時并無顯著體現(xiàn),說明圖16(a)中不同空化模型計算結(jié)果的差別可能來自不同模型的數(shù)值誤差,空化區(qū)域不同的差別并不會引起場結(jié)果及結(jié)構(gòu)響應(yīng)的顯著差別。

3 結(jié) 論

本文建立了空化截斷、修正的Schmidt和等熵空化流三種空化模型,并將其結(jié)合在龍格庫塔間斷伽遼金數(shù)值框架下。開展三種不同空化模型應(yīng)用在水下沖擊波與浮在自由面平板相互作用問題,水下沖擊波與彈性支撐平板相互作用問題及水下沖擊波與水背襯平板相互作用問題的對比研究。研究結(jié)果表明:

(1) 當計算壓力小于飽和蒸汽壓時,截斷模型的壓力保持在飽和蒸汽壓,而修正的Schmidt模型和等熵空化流模型隨著密度降低壓力迅速下降。修正的Schmidt模型在較小的密度變化條件下,下降得更加迅速,而等熵空化流模型表現(xiàn)出較為平滑的曲線。在壓力小于飽和蒸汽壓且大于修正Schmidt模型的截斷壓力pε時,相同密度條件下,修正的Schmidt模型獲得的壓力小于等熵空化流模型。

(2) 當流體環(huán)境壓力較大時,由于飽和蒸汽壓遠遠小于流體壓力,因此,不同空化模型的差別難以對整體計算結(jié)果起到顯著影響。當流體環(huán)境壓力不大,接近大氣壓時,截斷模型計算的密度結(jié)果在空化中心處略大于修正的Schmidt模型和等熵空化流模型的結(jié)果,不同空化模型計算的壓力結(jié)果在流體壓力下降到空化壓力之后產(chǎn)生差別。截斷模型的壓力保持為恒定的壓力,修正的Schmidt模型在稀疏波過后壓力即快速下降,而等熵空化流模型在低壓區(qū)域比修正的Schmidt模型更小。總體來講,修正的Schmidt模型與等熵空化流模型計算結(jié)果差別很小,而截斷模型與另外兩個模型計算結(jié)果在空化區(qū)域附近略有差別,但并不顯著。

(3) 影響不同空化模型計算結(jié)果的主要因素包括飽和蒸汽壓和峰值壓力。當飽和蒸汽壓很小,不同空化模型計算結(jié)果差別很小。

(4) 當飽和蒸汽壓較大且峰值壓力很大時,不同空化模型計算結(jié)構(gòu)響應(yīng)差別很小,而空化區(qū)域演化會由于時間步長較多時累積的數(shù)值誤差出現(xiàn)差別。對于水背襯空化邊界,修正的Schmidt模型繼續(xù)以原來的傳播速度向遠處傳播,等熵空化流模型則以相對緩慢的速度繼續(xù)向遠處傳播,而截斷模型則向結(jié)構(gòu)端方向潰滅,近結(jié)構(gòu)端的潰滅邊界三種模型并無差別。

(5) 當飽和蒸汽壓較大且峰值壓力很小時,不同空化模型計算結(jié)構(gòu)響應(yīng)差別很小,而空化區(qū)域演化和場輸出結(jié)果可能會由于數(shù)值誤差出現(xiàn)差別。修正的Schmidt模型計算遠結(jié)構(gòu)端空化域的傳播以及靠近結(jié)構(gòu)的潰滅邊界向遠處傳播速度比截斷模型和等熵空化模型計算結(jié)果快。