基于CEL法的破損邊界附近氣泡動力學特性研究

張 嘯, 崔 杰, 李明遠, 秦孜凱

(江蘇科技大學 船舶與海洋工程學院,江蘇 鎮江 212003)

現代戰爭中,船舶在受到攻擊破損后仍具備一定的抗沉性和破艙穩性,在戰場上仍具有作戰能力,當其受到二次打擊時,可喪失戰斗能力。Cui等[1-2]基于控制變量法,利用高速攝影技術和電火花氣泡試驗裝置研究了背空氣面的水下爆炸問題,發現一定的距離參數和破口參數,氣泡呈現中部凹陷并產生對射流的特殊氣泡運動特征。賀銘等[3]通過電火花試驗驗證了歐拉有限元法解決雙層結構附近水下爆炸問題的有效性,總結了不同破口參數、不同起爆位置和不同殼間水位條件下的耦和作用規律,發現了破口參數小于0.5時,內艙室會出現二次涌流現象。Huang等[4]采用了5點最小二乘方法來保證破口處自由液面的光滑和穩定的運動,討論了氣泡到破口的距離、破口的大小、破口的初始深度以及重力的變化對氣泡運動特性的影響。張之凡等[5]基于耦合歐拉-拉格朗日(coupled Eulerian-Lagrangian, CEL)法討論了不完整直角結構近自由面附近氣泡的動態特性,發現了在結構與自由面的共同作用下氣泡出現射流偏折、水冢一側凹陷和氣泡破碎等現象。盛振新等[6]通過小當量水下爆炸試驗對轟爆產物沖擊帶破口雙層板結構進行研究,總結了炸藥量、破口半徑和艙室寬度對氣泡運動和內板壁壓的影響規律。

在數值模擬方面,邊界元法、有限體積法、有限元法、CEL法、光滑粒子流體動力學法(smoothed particle hydrodynamics, SPH)等方法百家爭鳴。王瑩等[7]基于任意拉格朗日-歐拉法建立了爆破破冰數值模型,與試驗結果吻合較好。劉云龍等[8]基于歐拉有限元方法和聲學有限元場分離方法,自主編程實現了近、中、遠場下的水下爆炸對結構的響應。王平平等[9]基于SPH法實現了水下爆炸氣泡沖擊波、脈動、與結構相互作用的數值仿真。姜忠濤等[10]基于CEL法討論了近場水下爆炸射流在和對船體外板的沖擊過程,總結了射流的載荷壓力特性、速度分布及結構剪應力分布特性。目前,計算機技術在水下爆炸氣泡上的應用已趨于成熟,但由于破艙結構附近交界面復雜,呈現大變形、多介質等情況,計算精度和計算時間很難保證。在試驗方面,開展實船試驗可以獲得很好的研究效果,但是實船試驗的成本過于昂貴。目前主要通過電火花試驗,結合高速攝影技術進行研究,但電火花試驗由于氣泡尺度較小,其相似性還無法滿足。目前對于破損艙室遭受二次爆炸的研究有限,且主要局限于簡單的邊界條件,如平板、雙層板、圓柱體。破損壁面附近的氣泡為氣-液-固三相耦合,氣泡射流沖擊下結構的變形和損傷十分復雜。

針對以上問題,本文將破艙壁面簡化為具有圓形破口的彈性壁面,基于CEL法建立了水下爆炸數值模型,并結合模型試驗驗證了數值方法。研究了氣泡射流、氣泡體積、脈動周期,以及氣泡對壁面沖擊損傷的影響,探明了不同破口參數下的氣泡的脈動規律和射流穿透破口、撕裂等特征規律。為水下高能武器的研發和艦船的抗沖擊提供技術支撐。

1 CEL法

在氣泡動力學問題的研究中,與傳統的拉格朗日方法和歐拉方法相比,CEL法結合了拉格朗日網格與歐拉網格的優點,數值模擬中拉格朗日單元的節點隨物質點同步移動,而歐拉單元的節點在空間中固定不動;拉格朗日單元可以在歐拉網格內變形和運動,并且與含有材料的歐拉單元發生耦合,流體網格對結構網格施加壓力,結構網格則限制流體材料的流動。可以有效地解決水下爆炸問題中有關大變形和材料破壞等諸多問題。

在數值模擬計算過程中,將結構劃分為拉格朗日區域,流體劃分為歐拉區域,在兩種區域內分別采用歐拉有限元方法和拉格朗日有限元方法分別求解。歐拉有限元方法和拉格朗日有限元方法所使用的控制方程為流體力學中的質量守恒、動量守恒和能量守恒方程,如式(1)～(3)所示。兩者之間本質區別在于：歐拉方法是對計算域的離散,而拉格朗日方法是對物質的離散。

(1)

(2)

(3)

式中：fj為單位質量的張力;e為單位質量的內能;σij為結構應力張量;qi為熱通量;ρ為參考結構密度。

1.1 狀態方程

在數值模擬中,材料的各項屬性通常使用狀態方程來描述,本文中所使用到的材料有炸藥(PETN)、水、Q235碳素鋼。其中,炸藥采用JWL狀態方程,水的狀態方程采用線性Us-Up方程進行描述,整個模型采用體積分數法對氣泡進行材料劃分和定義。

1.1.1 炸藥

本文采用的炸藥為PETN,炸藥的狀態方程用JWL方程描述

(4)

式中：η=ρ/ρo,A、B、R1、R2、ω為炸藥的各個常數;ρo為炸藥密度;E為炸藥能量密度;ρ為水密度,取1 024 kg/m3,如表1所示。

表1 PETN炸藥的狀態參數

1.1.2 水

本文采用的不可壓縮黏性流體為水,其狀態方程用線性Us-Up方程描述

(5)

式中：Co為流體中的聲速;Γo為Gruneisen常數;ρo為流體的密度;S為聲速和沖擊波的傳播速度之間的線性系數;Em為流體單位質量的內能。水狀態參數如表2所示。此外,水的比熱為4 130 kJ/(kg·℃);動力黏性系數η為0.001 Pa·s。

表2 水的狀態參數

1.1.3 破損壁面

破損壁面為彈性鋼板,各參數取值如表3所示。

表3 破損壁面參數取值

1.2 CEL法有效性驗證

本文基于ABAQUS軟件,建立了水下爆炸CEL數值模型,其中,歐拉域尺寸為12 m×7 m×12 m;鋼板大小為厚度0.02 m,長0.8 m的正方形板。炸藥放置在水下1 m,大于三倍氣泡半徑,因此忽略自由液面對氣泡的影響,歐拉域僅包含水、炸藥、鋼板三種材料。正方形鋼板垂直于破口方向的邊被剛性固定,歐拉域的邊界上設置邊界條件無反射以模擬無限水域,設置無流入、速度為0,以模擬無窮遠處邊界條件。為了驗證本文建立的數值方法正確性,與水下爆炸模型試驗[11]進行了對比。

圖1為破口參數Rh=0.25,距離參數γw=0.7下水下爆炸試驗與數值模擬氣泡運動特性對比圖,上方一組視圖(a)為水下爆炸試驗,下方一組試圖(b)為數值模擬結果。其中,氣泡脈動時間標注在每幀圖片的左上方,單位為ms。觀察發現,數值模擬結果與水下爆炸試驗在各個階段保持較好的一致性。本文將基于與試驗相同的距離參數研究破口對氣泡動力學行為的影響。

(a) 模型試驗

(b) 數值模擬

試驗氣泡最大等效半徑為0.267 m,第一周期為45.7 ms,數值模擬計算的氣泡最大等效半徑為0.262 m,誤差為4.81%,周期為43.5 ms,誤差為1.87%。氣泡表面頂點位移對比圖,如圖2所示。由圖2可知,數值計算結果與模型試驗吻合度較高,本文建立的數值方法具有較好的精度。

圖2 氣泡表面頂點位移對比圖

2 破口參數對氣泡脈動、射流以及壁面損傷研究

2.1 數值模型

實際海戰中,影響氣泡運動的因素諸多,為了使研究具有普遍性,本文將破口簡化為圓形破口,研究破口參數Rh對其運動及損傷影響。

圖3為水下爆炸數值模擬中各參數定義,Rh為破口半徑;Dw為氣泡中心點到破口中心點距離,Df為氣泡中心點與自由液面距離,無量綱過程如式(6)所示

γw=D/Rmax

(6)

式中：D為有量綱距離;Rmax為氣泡最大體積時刻的等效半徑。

圖3 試驗示意圖

2.2 氣泡動態特性

水下爆炸氣泡在不同破口參數下會呈現不同的脈動特性,為了研究破口參數對于氣泡的射流、運動特性、氣泡體積、周期等因素的影響,本文設置了14組工況,破口半徑(無量綱)依次為0、0.05、0.10、0.15、0.20、0.25、0.30、0.35、0.40、0.50、0.60、0.70、0.90、1.10。

2.2.1 小型破口下(Rh≤0.10)氣泡脈動特性

圖4為三種小破口參數下的氣泡脈動特性,分別為完整壁面工況(a)Rh=0、工況(b)Rh=0.05、工況(c)Rh=0.10。為方便分析,將氣泡運動特征形態采用數字幀的形式標記在圖片上方(如1～7),時間標注在每張圖片的左上方,單位為ms。為方便觀察射流,對部氣泡形態進行透視處理。

(a) Rh=0

(b) Rh=0.05

(c) Rh=0.10

第1～3幀為氣泡的膨脹階段,在小型破口參數下各工況脈動形態相似,觀察體積最大時刻(第3幀),發現了從右向左的射流。這一現象的出現,是由于破口周圍的氣泡,在缺少壁面阻攔的情況下受到水流的沖擊,此時氣泡內的壓力最小,遠低于破口右側流場域壓力,故射流方向由右向左。破口參數Rh=0.05這一工況尤為特殊,主要表現在該射流細小,依靠氣泡左側朝向壁面收縮方能穿透氣泡左表面。該現象出現的原因是破口較小,受到的沖擊力有限。第4～6幀為氣泡的收縮階段,第5幀均出現了從左向右的射流,該射流是由于壁面的阻擋在氣泡左側形成高壓區,同時受重力的影響,形成斜射流。為了方便分析,定義由右向左的射流為射流1(第4幀),由左向右的為射流2(第5幀)。在射流2和射流1的共同作用下,會在氣泡左側“撕裂”出部分氣泡,如圖4(c)破口參數Rh=0.10第6幀所示。

2.2.2 中型破口下(Rh=0.15～0.50)氣泡脈動特性

圖5為中型破口下(Rh=0.15～0.50)氣泡脈動特性,破口參數Rh依次為0.15、0.20、0.25、0.30、0.35、0.40、0.50,上述破口參數下破口的存在對壁面Bjerknes力產生了削弱作用,但壁面對氣泡仍存在一定的吸引力。第1～4幀為氣泡的膨脹階段,第2幀氣泡近壁面端向破口內產生“凸進”現象。“凸進”是指在破損壁面的影響下,當氣泡內氣體壓力發生變化時,氣泡脈動過程向壁面方向擠壓,相對于完整壁面,其破口位置處對氣泡脈動過程限制較小,使得近破口位置處氣泡向破口方向運動,甚至穿過破口結構,呈現獨特的尖嘴狀。研究發現,當破口尺寸較小時,破口近似于完整壁面,對氣泡運動的影響很小,尖嘴狀形態不明顯;反之當破口參數越大,氣泡右表面產生的“凸進”越明顯,即涌入破口的氣泡體積越大。氣泡“凸進”在穿過破口后,“凸進”內壓力得到釋放,依靠慣性繼續向右運動,破損壁面阻礙氣泡主體運動,在第3幀出現了“撕裂”的現象。在較小破口參數工況中,氣泡“凸進”還未穿透破口,在較大破口參數工況中,“凸進”較為圓潤,較難產生“撕裂”,因此“撕裂”現象僅存在于破口參數Rh=0.20～0.40的工況中。破口參數越大,“撕裂”發生的時間越晚。氣泡右表面在撕裂完成后開始收縮,收縮主要集中于破口附近的氣泡局部區域,第4幀形成了從右向左的射流1。第5～6幀為氣泡的收縮階段,第5幀觀察到從左向右的射流2,對比不同破口參數下的兩種射流發現,射流1的半徑略小于破口半徑。第一周期末,射流2會再次貫穿射流1導致的環形氣泡,將氣泡再次“撕裂”。第7幀為氣泡第二周期的脈動情況,在破口參數較大工況中,氣泡體積相對較小,這是因為較大破口工況“撕裂”部分氣泡會依靠慣性穿過破口脫離主體向右運動,使得氣泡二次膨脹體積減小。

2.2.3 大型破口下(Rh=0.70～1.10)的氣泡脈動特性

圖6為大型破口下(Rh=0.70～1.10)氣泡脈動特性,破口參數Rh依次為0.70、0.90、1.10,大型破口工況由于破口參數較大,壁面對于氣泡產生的Bjerknes力相對較弱。第1～5幀為氣泡的第一周期,在破口參數超過0.7后,氣泡周圍充滿流體,破口難以限制壁面右側的流體形成集中載荷沖擊氣泡產生射流1。破口參數Rh=0.90的工況(b)34.8 ms產生了向左的射流,該射流無法沖擊氣泡左表面,在氣泡坍塌最小時形成環形氣泡。在Rh=1.1的工況(c)中,第一周期沒有觀察到射流1的存在,第一周期末期形成了與自由場工況類似的斜向右上方的射流2。

(a) Rh=0.15

(b) Rh=0.20

(c) Rh=0.25

(d) Rh=0.30

(e) Rh=0.35

(f) Rh=0.40

(g) Rh=0.50

(a) Rh=0.70

(b) Rh=0.90

(c) Rh=1.10

2.3 水流沖擊氣泡產生的射流現象

在目前的水下爆炸氣泡試驗中,由于氣泡內部的結構較難觀察到,僅對于射流2進行了描述[11],然而研究發現射流1仍然會對氣泡的形態和破損壁面產生影響。本節總結了射流1在不同破口參數下的形態特征,如圖7所示,其中,破口參數標注在每個工況的左上方,時間標注在每個工況的左下方,單位為ms。

圖7 不同工況射流現象

當破口參數0

圖8和圖9描述了不同破口參數下的射流速度,圖8描述了破口影響下右側水流沖擊氣泡形成的從右向左的射流1和在第一周期末期氣泡左側高壓而形成的從左向右的射流2在不同破口參數下的射流速度,其中針劑型射流頂端“撕裂”導致針頭狀射流并不能代表射流的真實速度,選取針頭末端為監測點進行測量。射流1速度一般在6～12 m/s,破口參數越大,氣泡受影響的區域越大,射流1速度越快。射流2的速度一般在40～62 m/s,因此,氣泡左側高壓所導致的射流2仍然是壁面造成損傷的主要原因。此外,射流1的速度和作用區間對于射流2存在削弱作用,在破口參數Rh>0.5的工況中,射流2的方向轉化為由下向上甚至消失。圖9描述了不同破口參數下,射流1頂端速度的時間歷程曲線。射流1頂端“撕裂”形成的針頭狀射流速度顯著高于破損壁面右側流體沖擊氣泡形成的射流主體,使射流1變得復雜,最大速度一般出現在即將貫穿氣泡時刻。

圖8 不同破口參數下的射流速度對比圖

圖9 不同破口參數下的射流1頂端速度時歷曲線

2.4 破口參數Rh對氣泡和壁面損傷的影響

為了研究破口參數對于氣泡動態特性的影響、分別給出氣泡體積、第一周期、第二周期體積峰值、氣泡周期隨時間變化曲線,具體如圖10、圖11和圖12所示。

從圖10和圖11中可以看出,氣泡體積峰值在破口參數的影響下,整體呈現在大型破口和小型破口工況體積峰值大,中型破口工況體積峰值小的特點。這一特點在第一周期尤其明顯,最小體積出現在破口參數Rh=0.35的工況中,究其原因是右側流體沖擊氣泡形成的射流使右表面提前收縮,削弱了氣泡后期膨脹到最大體積,減少了氣泡體積峰值。

圖10 不同破口參數下氣泡體積時歷曲線

圖12 不同破口參數下氣泡體積極值發生時間

國內外學者對比了多組試驗,發現了壁面作用下的氣泡脈動周期更長[12]。從圖12中觀察到：破口參數越大,氣泡周期越短,且周期變短的趨勢減弱。這說明破口會削弱壁面的影響,從而增大破口參數會縮短氣泡脈動周期。

圖13描述了不同破口參數下的破口邊緣處變形峰值和發生的時間。從圖13可以發現：破口參數較小時,含有破口壁面變形較大,加劇壁面損傷;破口參數較大時,壁面變形值較小,壁面的最大變形出現在破口參數Rh=0.1的工況中。此外,破口參數增加,會滯后破損壁面最大位移到來的時間。

圖13 不同破口參數下破口邊緣變形峰值和時間

3 結 論

本文基于CEL法建立了破損壁面附近水下爆炸模型,將數值結果與水下爆炸試驗進行對比,驗證了數值方法的有效性。對不同破口參數下的破損壁面附近氣泡進行數值計算,結論總結如下：

(1) 受含有破口壁面的影響,破口附近氣泡運動朝向破口甚至穿過破口呈現尖嘴狀的 “凸進”現象;當破口尺寸較小時,對氣泡運動的影響很小,尖嘴狀形態不明顯;反之當破口尺寸較大時,含有破口壁面對氣泡膨脹限制較弱, “凸進”現象通常出現在破口Rh=0.10～0.90的工況中。

(2) 氣泡受到破口右側水流沖擊而產生的從右向左的射流,射流出現在Rh=0.05～1.10的工況中,射流寬度隨破口參數增大而增大,略小于破口半徑,按照其形狀可總結成四種射流現象,分別為直線型、燈塔型、針劑型、沙堆型。

(3) 破口參數對氣泡脈動和射流影響較大,氣泡體積峰值隨著破口參數的增加先增大后減小;破口參數較小時,含有破口壁面變形較大,加劇壁面損傷;氣泡射流2的速度大于射流1的速度。