王平平 張阿漫, 彭玉祥 孟子飛

* (哈爾濱工程大學船舶工程學院,哈爾濱 150001)

? (中山大學海洋工程與技術學院,廣東珠海 519000)

引言

近場水下爆炸主要包括沖擊波、氣泡脈動和射流以及結(jié)構(gòu)毀傷等.當炸藥在水中發(fā)生爆炸,對周圍水介質(zhì)產(chǎn)生劇烈的擠壓作用,造成水介質(zhì)的密度、壓力以及速度急劇上升,在水中產(chǎn)生強間斷沖擊波.炸藥引爆后形成的氣態(tài)爆轟產(chǎn)物被水包圍形成氣泡,在內(nèi)部高壓和外部流場壓力作用下發(fā)生膨脹和收縮,此過程被稱為氣泡脈動.氣泡脈動頻率與船體結(jié)構(gòu)的固有頻率往往非常接近,因此容易引起船舶結(jié)構(gòu)發(fā)生共振,對其總縱強度非常不利,甚至導致艦船折斷.氣泡在自由液面以及艦船結(jié)構(gòu)表面等邊界耦合作用下還會形成高速氣泡射流,其速度可達數(shù)十米甚至上百米每秒,使艦船結(jié)構(gòu)產(chǎn)生嚴重破壞.目前,水下武器的引爆多采用非觸發(fā)引信的形式,在艦船附近發(fā)生爆炸,因此研究近場水下爆炸問題對保護艦船生命力具有重要意義.

目前水下爆炸數(shù)值模擬大多基于商業(yè)軟件或自編程序的有網(wǎng)格方法,對于解決水下爆炸中的復雜多相界面捕捉和結(jié)構(gòu)大變形問題常常面臨較大困難.因此,許多學者將目光轉(zhuǎn)向了新興的無網(wǎng)格算法,其中以光滑粒子流體動力學(smoothed particle hydrodynamics,SPH) 方法最為成熟和常用.由于SPH 方法在求解強非線性大變形和動邊界問題時避免了網(wǎng)格畸變,省去了界面捕捉等復雜過程,因此在處理移動界面和大變形問題上具有先天優(yōu)勢.下面將從沖擊波、氣泡以及結(jié)構(gòu)毀傷三個方面對當前的SPH 研究現(xiàn)狀進行介紹.

針對水下爆炸沖擊波研究,2002 年Liu 等[1]首次將傳統(tǒng)人工黏性SPH 方法應用到水下爆炸問題的仿真中,模擬了二維有限方形域中水下爆炸沖擊波的傳播和反射等現(xiàn)象,較好地捕捉了沖擊波的幅值和位置.同樣基于傳統(tǒng)人工黏性SPH 方法,楊剛等[2]模擬了二維近自由面水下爆炸和炸藥爆轟過程,并研究了空氣對壓力場和密度場的影響.采用傳統(tǒng)人工黏性SPH 方法,明付仁[3]模擬了三維近自由面水下爆炸,考察了自由面下不同反射區(qū)沖擊波的波形,并探究了不同爆深時近自由面沖擊載荷特性.但需要指出的是,傳統(tǒng)人工黏性SPH 方法在處理流場間斷時常常導致劇烈的非物理振蕩或引起過大的數(shù)值耗散.為了更好地處理流場的間斷,獲得更加精確的數(shù)值結(jié)果,一些學者[4-7]將黎曼思想引入SPH中得到黎曼SPH 方法,并將其應用到水下爆炸模擬中.例如,蔡志文[8]采用一階黎曼SPH 方法模擬了二維水下爆炸問題,并探究了不同的近似黎曼求解器對壓力場的影響.可以看出,已有的水下爆炸沖擊波SPH 模擬多局限在二維框架下,且大多數(shù)研究采用了傳統(tǒng)的人工黏性力來實現(xiàn)對沖擊波的捕捉,但人工黏性力的大小常需要通過人工參數(shù)調(diào)節(jié),過大的人工黏性將會導致沖擊波峰值的降低,而過小的黏性力會引起劇烈的數(shù)值振蕩.雖然有少數(shù)研究采用了黎曼SPH 算法,但由于未采用重構(gòu)算法往往導致間斷附近出現(xiàn)較大的數(shù)值耗散和非物理振蕩.上述原因?qū)е乱延械腟PH 方法得到的水下爆炸沖擊波往往無法達到令人滿意的精度.

對于水下爆炸氣泡動力學的研究,盡管網(wǎng)格方法取得了很大成功[9],但其在處理復雜界面時仍然存在較大挑戰(zhàn),需要精細的網(wǎng)格分辨率和高精度的計算格式.考慮到SPH 方法在處理多相界面、流體大變形等問題上具有顯著優(yōu)勢,近年來一些學者將SPH 方法應用到水下爆炸氣泡數(shù)值模擬中.然而,由于水下爆炸氣泡的模擬往往涉及到強可壓縮、大密度比、粒子體積變化過大以及計算效率低等問題,因此相關的SPH 研究非常有限.Liu 等[10]采用傳統(tǒng)人工黏性SPH 方法模擬了二維水下爆炸氣泡在密閉室內(nèi)的運動,并研究了此過程中沖擊波的傳播以及氣泡尺寸和形狀的演化.Pineda 等[11]基于黎曼SPH 模擬了自由場和剛性壁面附近的二維氣泡坍塌,并分析了包括壓力波和射流等現(xiàn)象的主要特征,但所采用的黎曼SPH 數(shù)值耗散較大,且未對氣泡膨脹和收縮過程中粒子體積變化過大的問題做出處理.為了解決粒子體積變化過大導致的精度下降問題,Sun 等[12-13]在傳統(tǒng)人工黏性SPH 的框架下提出了粒子體積自適應方案,該方案通過對粒子進行分割和合并以降低粒子的體積變化,模擬了自由表面和剛性壁面附近氣泡的動力學特性.同樣為了解決粒子體積變化過大的問題,Peng 等[14]基于一階黎曼SPH 方法提出了粒子再生技術,并模擬了剛性壁面附近和自由面附近的二維氣泡坍塌.總的來看,目前基于無網(wǎng)格SPH 方法的水下爆炸氣泡研究大多局限于二維情況,且多數(shù)采用了傳統(tǒng)的人工黏性SPH 多相流模型.該模型對于水氣界面的密度間斷處理常帶來壓力不穩(wěn)定等問題.為了使水氣界面穩(wěn)定,一些研究對于氣相采用遠大于物理聲速的人工聲速,導致無法考慮真實的氣體可壓性,且穩(wěn)定時間步非常小,大大降低了計算效率.

對于近場水下爆炸的結(jié)構(gòu)毀傷模擬,近年來無網(wǎng)格SPH 方法也得到了初步應用.張阿漫等[15]采用基于體積近似的SPH 方法在二維框架下模擬了水下接觸爆炸對單層和雙層鋼板的毀傷.文獻[3,15-16]采用SPH-SPS 耦合算法模擬了水下接觸爆炸對簡單平板的毀傷,并探究了不同起爆方式對結(jié)構(gòu)毀傷程度的影響.此外還采用SPH-FEM 耦合算法模擬了接觸爆炸對板架結(jié)構(gòu)的毀傷.彭玉祥[17]采用SPH算法模擬了近場水下爆炸對圓板、雙層殼結(jié)構(gòu)以及舷側(cè)結(jié)構(gòu)的破壞.上述SPH 水下爆炸結(jié)構(gòu)毀傷模擬雖然取得了一些成果,但需要指出的是由于計算效率和精度的限制,這些研究主要考慮了水下爆炸前期沖擊波對結(jié)構(gòu)的破壞,而忽略了水下爆炸氣泡與結(jié)構(gòu)的相互作用,因此對水下爆炸結(jié)構(gòu)毀傷的研究并不完整.而且,這些研究的沖擊波載荷均基于傳統(tǒng)的人工黏性SPH 方法求解得到,導致載荷誤差較大.

針對上述近場水下爆炸SPH 模擬存在的不足,本文推導了適用于間斷流場的多相流黎曼SPH 數(shù)值模型,開發(fā)基于RKPM 的結(jié)構(gòu)損傷斷裂求解器,并應用法向通量邊界條件構(gòu)建SPH-RKPM 流固耦合數(shù)值模型,實現(xiàn)對近場水下爆炸沖擊波、氣泡以及結(jié)構(gòu)毀傷的高精度數(shù)值求解,并給出了水下爆炸載荷及其對結(jié)構(gòu)的毀傷特性,以期為近場水下爆炸載荷預報提供理論和基礎性技術支撐,為毀傷威力評估和艦船防護結(jié)構(gòu)設計提供參考.

1 基本理論

1.1 黎曼SPH 間斷流場模型

在SPH 方法[18]中,對于無黏流體,其拉格朗日形式的控制方程如下

其中,e,p,g,u,ρ分別代表流體單位質(zhì)量的內(nèi)能、壓力、重力加速度、速度和密度.通過核近似和粒子近似,可得到SPH 的離散控制方程如下

式中,Wij=W(|xi-xj|,hi) 為SPH 使用的核函數(shù),常需要滿足歸一性、緊支性以及狄拉克性質(zhì),核函數(shù)的有效范圍稱為支持域,支持域的半徑為kh,其中k是一個與所取核函數(shù)有關的常量,h稱為光滑長度.本文中,如不特別說明均采用Wendland 核函數(shù)[19],且光滑長度取為2 倍粒子間距,即h=2Δx.下標i和j分別代表目標粒子和其支持域內(nèi)的粒子.可以看出,核函數(shù)值的大小與粒子之間的距離以及光滑長度有關.

為了使控制方程(1)封閉,一般采用狀態(tài)方程來表述流體的壓力和密度、內(nèi)能之間的關系.例如本文2.3 節(jié)、3.2 節(jié)和3.3 節(jié)中認為水是正壓流體[20],采用如下Tait 狀態(tài)方程

其中,絕熱系數(shù) γw=7.15,參考密度 ρref=998.0 kg/m3.考慮大氣壓力時參考壓力pref=101325 Pa,否則pref=0.0 Pa.在3.1 節(jié)、3.4 節(jié)和3.5 節(jié)中需要考慮水內(nèi)能變化對其壓力的影響,即水為斜壓流體,采用Mie-Gruneisen 狀態(tài)方程,其具體形式和參數(shù)參考文獻[10].

對于氣體,在3.2 和3.3 節(jié)中認為其是正壓[21],采用Cole 狀態(tài)方程

其中,絕熱系數(shù) γa1=1.25,pref=101325 Pa,ρref=2.0 kg/m3.在2.1 節(jié)中認為其是斜壓氣體,采用理想氣體狀態(tài)方程[22],如下

其中絕熱系數(shù) γa2取為1.4.此外,對于3.1 節(jié)、3.4 節(jié)和3.5 節(jié)中TNT 藥包爆轟產(chǎn)生的氣體,通過JWL 狀態(tài)方程描述其壓力和密度、比內(nèi)能之間的關系,詳見文獻[10].

在傳統(tǒng)的SPH 方法中,為了提供維持數(shù)值穩(wěn)定的耗散,常常在動量與能量方程中加入人工黏性力項,但需要指出的是,這種方法對于流場間斷的處理常常面臨較大困難.具體來說,在間斷附近,人工黏性力常常導致劇烈的非物理振蕩,而在遠離間斷處,可能引起過大的數(shù)值耗散.為了解決這一問題,一些學者將黎曼思想引入SPH 中得到黎曼SPH 方法,該方法認為每兩個相互作用的粒子對構(gòu)成一維黎曼問題,且粒子對的連線中點存在間斷,通過求解一維黎曼問題得到粒子之間的相互作用,進而得到整個計算域的解[23].基于上述思想,控制方程(2)變?yōu)槿缦滦问?/p>

式中,p*和u*分別為中間變量,由黎曼求解器得到.本文中采用了適用于大密度比多相流的PVRS 近似黎曼求解器[24],如下

式中,下標l和r分別表示黎曼問題的左右狀態(tài),c表示流體的聲速,Z為流體的阻抗,計算為Z=ρc.標量u*和中間變量u*的關系為

其中,eij=(xj-xi)/|xj-xi| 表 示由粒子i指向j的單位矢量.

不同于傳統(tǒng)人工黏性力SPH 方法,黎曼SPH 通過黎曼求解器的固有數(shù)值黏性而非顯式人工黏性來保證數(shù)值格式的穩(wěn)定性,因此避免了相關的可調(diào)人工參數(shù).更重要的是,由于采用黎曼求解器來確定粒子間的相互作用,黎曼SPH 更擅長處理間斷問題,因此對于大密度比強間斷多相流的模擬具有更高的精度.但同時,原始的黎曼SPH 方法直接將相互作用的粒子的物理量作為黎曼問題的初始左右狀態(tài),因此僅具有一階精度,常常引起過大的數(shù)值耗散,導致精度下降.為解決這一問題,一些學者提出了適用于SPH 方法的黎曼高階重構(gòu)算法,例如MUSCL 重構(gòu)[24-25]、WENO 重構(gòu)[26-28]以及TENO[22]重構(gòu).對于本文中模擬的水下爆炸問題,選取了MUSCL 重構(gòu)算法[25],該算法采用線性重構(gòu)得到黎曼問題左右初始狀態(tài),對于變量 φ 的重構(gòu)如下

其中,Δ φi和 Δ φj由下式求得

式中,常數(shù) β 本文中取為1.5,變量 φ 的變化量 Δ φi和Δφj由下式計算得到

其中,矢量rji計算為rji=xj-xi.

1.2 RKPM 結(jié)構(gòu)損傷斷裂模型

上節(jié)中介紹了基于SPH 方法的流體計算模型,對于結(jié)構(gòu)損傷斷裂模型,也可以采用SPH 方法來建立.但由于SPH 方法中采用的是強形式的控制方程,導致自由邊界無法自動滿足,對于存在大量自由邊界的復雜結(jié)構(gòu)模擬存在較大困難,同時也不利于處理結(jié)構(gòu)裂紋拓展中動態(tài)增加的內(nèi)部自由邊界.因此本文中采用基于弱形式控制方程的重構(gòu)核粒子法(reproducing kernel particle method,RKPM)來模擬結(jié)構(gòu)的彈塑性以及損傷斷裂,目前RKPM 已經(jīng)被廣泛應用到艦船結(jié)構(gòu)的爆炸毀傷響應預報當中[29-30].本文采用RKPM 殼單元來對結(jié)構(gòu)進行建模,采用退化實體幾何表述,僅需在殼中面布置一層粒子來對結(jié)構(gòu)進行離散,依據(jù)虛功率原理,得到殼體的控制方程如下

式中,ρ0為材料的初始密度,η 表示殼體的厚度,和分別為I粒子的速度和轉(zhuǎn)動角速度的時間導數(shù).MI和Q分別是RKPM 型函數(shù)和第一P-K (Piola-Kirchhoff)應力張量,和 ω 分別為結(jié)構(gòu)牽引力以及單位質(zhì)量體積力.表示型函數(shù)MI與殼厚度方向的參數(shù)坐標 ψ 的乘積,即=ψMI.zI表示I粒子的纖維矢量,V0為I粒子支持域的殼體,λ0為對應的殼體參考面,為牽引力邊界.同時針對復雜結(jié)構(gòu)的大變形損傷,需要結(jié)構(gòu)的大變形接觸檢測以及接觸力計算模型,參考文獻[30],針對結(jié)構(gòu)裂紋萌生及拓展,相關的材料彈塑性損傷本構(gòu)模型,以及粒子的裂紋萌生拓展算法,可以參考文獻[31].

1.3 流固耦合算法

在流固耦合界面上,需要滿足動力連續(xù)和運動連續(xù)條件,也就是說界面上壓力以及沿界面法向的速度要連續(xù).為了滿足上述條件,本文將結(jié)構(gòu)粒子作為流體粒子的移動邊界,并依據(jù)Leffe 等[32]提出的法向通量邊界條件,來計算結(jié)構(gòu)對流體的作用.此時,考慮結(jié)構(gòu)作用的流體控制方程變?yōu)?/p>

其中,下標b代表流體粒子i鄰域內(nèi)的結(jié)構(gòu)粒子,nb和Sb分別為結(jié)構(gòu)粒子的單位法向矢量和面積.

另一方面,對于結(jié)構(gòu)粒子,通過對周圍流體粒子壓力的SPH 插值得到其受到的流體載荷pfb,具體插值計算公式如下[30]

上述算法可以很好地實現(xiàn)流體粒子與單層殼結(jié)構(gòu)粒子的流固耦合,從而建立起SPH-RKPM 流固耦合算法,對于處理水下爆炸瞬態(tài)流固耦合問題具有良好的精度[32-33].

2 數(shù)值驗證

為驗證所建立的流場和結(jié)構(gòu)數(shù)值模型以及流固耦合算法的數(shù)值精度,本節(jié)分別模擬了激波管、平板大變形撕裂和水射流沖擊彈性板三個算例,并將得到的結(jié)果與解析解或其他數(shù)值解進行對比.

2.1 激波管

首先,為了驗證所建立的黎曼SPH 方法對激波的捕捉精度,本文首先求解了經(jīng)典的一維“激波-密度波相互作用問題”[34-35].此算例的初始條件為

可以發(fā)現(xiàn),此算例中初始時刻壓力、密度和速度在x=1處均呈間斷分布,這將導致密度場中產(chǎn)生馬赫數(shù)為3 的激波和正弦波,因此可以被用來檢驗數(shù)值模型捕捉間斷的能力.

為驗證當前SPH 方法的收斂性和精度,采用0.05,0.025 和0.0125 三種粒子分辨率對計算域進行離散,得到的黎曼SPH 結(jié)果如圖1 所示,可以看到無論是密度還是速度分布曲線,隨著粒子分辨率的提升,SPH 的結(jié)果與文獻[36]中的經(jīng)典五階WENO方法得到的參考解的差別減小.當粒子分辨率為0.0125時,SPH 結(jié)果與參考解高度吻合,從而證明了當前提出的黎曼SPH 方法對于此類可壓縮間斷流動問題,具有良好的精度和收斂性,為后續(xù)模擬水下爆炸問題奠定了基礎.

圖1 黎曼SPH 方法得到的“激波-密度波相互作用問題”在t=1.8時刻的密度和速度分布Fig.1 Density and velocity distributions of the shock-density wave interaction at t=1.8 obtained by the present Riemann SPH scheme

2.2 圓管受壓屈曲

針對RKPM 結(jié)構(gòu)計算模型,采用一個外壓作用下的圓柱管的非線性動力屈曲算例來進行驗證.在圖2 中給出了結(jié)構(gòu)的幾何參數(shù)、邊界條件設置以及相應的材料參數(shù).計算中對于材料的彈塑性響應采用各向同性線性硬化假設,將粒子間距設置為0.02 m,共采用5940 個粒子對圓管進行幾何離散.

圖2 圓管受壓動力屈曲計算模型Fig.2 Computational model of the dynamic buckling of the circular tube under compression

文獻[37]中Kyriakide 的實驗結(jié)果給出了圓管會在第三種屈曲模式中被壓潰.通過屈曲分析得到的解析解表明圓管在該模式下壓潰所需要的外界壓力為30.43 MPa,在實驗中測得的臨界壓力為28.26 MPa.采用RKPM 計算得到的該模式下屈曲所需的外界壓力載荷為29.64 MPa.

在計算中,圓管中間部分的加載壓力在1.0 ms內(nèi)線性增加,然后維持在29.64 MPa,在1.5 ms 對圓柱殼的圓周方向坐標加入一個縮放擾動1+cos(3θ)/20,在第三種屈曲模式下對粒子坐標進行擾動.最終計算得到圓管的動力屈曲過程如圖3 所示,從圖中可以看到計算在2.0 ms 圓柱殼初步產(chǎn)生了第三種屈曲模式的應力分布,然后在2.2 ms 結(jié)構(gòu)就發(fā)生了動力屈曲,到2.5 ms 圓管就被完全壓潰,殼體壓縮到發(fā)生接觸,之后結(jié)構(gòu)的變形變得穩(wěn)定.

圖3 圓管受壓動力屈曲過程Fig.3 Dynamic buckling process of the circular tube under compression

在圖4 中將RKPM 計算得到圓管最終屈曲變形構(gòu)型與實驗結(jié)果進行了對比,從圖中可以看到數(shù)值計算得到屈曲模式與實驗結(jié)果吻合良好,驗證了RKPM 殼體大變形數(shù)值計算模型的有效性.實驗中圓管一端結(jié)構(gòu)出現(xiàn)了斷裂,但數(shù)值模擬當中未考慮材料的損傷,因此殼體沒有出現(xiàn)斷裂.

圖4 圓管屈曲變形數(shù)值結(jié)果(上)與文獻[37]中的實驗結(jié)果(下)對比Fig.4 Comparison of the buckling deformation of the circular tube between numerical results (upper) and the experimental results in Ref.[37] (lower)

2.3 水射流沖擊彈性板

在本節(jié)中,通過模擬高速水射流沖擊彈性平板算例,以驗證當前采用的流固耦合算法在強沖擊問題中的精度.其計算模型如圖5 所示,一個半徑為1.5 mm 的球頭柱形水射流,其長度為12 mm,以570 m/s 的初始速度沖擊一個四周剛固的方形板中心,方形板的長度和厚度分別是30 mm 和5.9 mm,其材料參數(shù)為:密度1170 kg/m3,彈性模量3.3 GPa,泊松比0.36.分別采用本文提出的SPHRKPM 方法和軟件ABAQUS 中基于網(wǎng)格的耦合歐拉-拉格朗日(coupled Eulerian-Lagrangian,CEL)方法模擬此算例,其中粒子間距或網(wǎng)格大小均取為0.05 mm.

圖5 水射流沖擊彈性板計算模型Fig.5 Computational model of the water jet impacting elastic plate

圖6 給出了CEL 和SPH-RKPM 得到的幾個典型時刻的結(jié)果,可以發(fā)現(xiàn),在撞擊到方板后,水柱頭部被壓平,形成一個緊貼方板的圓形水面,隨著其不斷朝向板運動,水柱長度不斷縮短,而圓形水面半徑迅速增大.可以發(fā)現(xiàn),CEL 方法只在模擬后期捕捉到少許飛濺的水滴,這是因為其采用了流體體積分數(shù)(volume of fluid,VOF)方法來重構(gòu)界面,當網(wǎng)格分辨率不足時,較小的飛濺水滴難以精細捕捉,而SPH方法由于其天然的拉格朗日特性,采用自由運動的粒子來代表流體,能更好地捕捉流體飛濺,此算例很好地體現(xiàn)了SPH 方法在模擬流體飛濺現(xiàn)象上的優(yōu)勢.

圖6 幾個典型時刻SPH-RKPM(左)和CEL(右)結(jié)果對比Fig.6 Comparison of SPH-RKPM (left) and CEL (right) results at several typical moments

為了定量地比較兩種方法的結(jié)果,圖7 給出了方板中心點位置處的流場壓力以及板的位移時歷曲線,并與文獻[38]中的實驗結(jié)果進行對比.可以發(fā)現(xiàn),當射流沖擊平板時,中心點的壓力迅速增大,其幅值約為0.9 GPa,隨后壓力降低并逐漸穩(wěn)定在約0.2 GPa.兩種數(shù)值方法得到的壓力峰值和變化趨勢與實驗結(jié)果均較為接近,但CEL 方法得到的壓力曲線在后期出現(xiàn)較為明顯的振蕩,而由于當前的SPH 采用了黎曼方法,因此得到的壓力變化更加平穩(wěn).由圖7(b)可以看出,在沖擊過程中中心點的位移不斷增大,在 0～2 μs 和2～12 μs 呈分段線性變 化,2 μs 之前增速較大,而 2～12 μs 增速變緩.大體而言,兩種方法得到的位移曲線較為吻合.通過圖7 給出的流場壓力和結(jié)構(gòu)響應驗證了當前采用的SPH-RKPM 流固耦合方法具有較高的精 度.

圖7 不同方法得到的圓板中心點位置處的(a)流場壓力和(b)結(jié)構(gòu)位移時歷曲線Fig.7 Time histories of (a) flow field pressure and (b) structural displacement at the central point of circular plate obtained by different methods

3 算例結(jié)果

3.1 水下爆炸沖擊波模擬

首先,基于所建立的黎曼SPH 流體動力學模型,對三維水下爆炸沖擊波傳播過程進行模擬.計算模型為邊長0.7 m 的方形水域,水域中心布置一個半徑為0.05 m 的球形TNT 藥包,爆點位于藥包中心,在距藥包中心0.3 m 和0.35 m 處分別布置有兩個測點A和B來記錄壓力時歷.整個計算域的粒子分辨率為6.25 mm.此算例中水采用Mie-Gruneisen 狀態(tài)方程,TNT 藥包采用JWL 狀態(tài)方程[39].為了驗證本文采用的黎曼SPH 方法在激波捕捉方面的優(yōu)勢,本節(jié)中分別采用傳統(tǒng)的人工黏性SPH 方法和黎曼SPH 方法,得到的三個典型時刻的壓力場分布如圖8所示.可以發(fā)現(xiàn)TNT 藥包起爆后,產(chǎn)生迅速向外傳播的壓力波.從壓力場看,兩種方法得到的結(jié)果存在以下幾點區(qū)別:首先,傳統(tǒng)人工黏性SPH 得到的壓力場在水氣界面存在壓力不連續(xù),且波陣面附近存在多個峰值,而黎曼SPH 方法得到的壓力場較為光順;其次,相比于傳統(tǒng)人工黏性SPH 方法,黎曼SPH 得到的波陣面壓力峰值略大,說明其數(shù)值耗散相對較小;最后,黎曼SPH 得到的波陣面更窄,具有更好的間斷特性.

圖8 分別采用人工黏性SPH 方法(左)和當前的黎曼SPH 方法(右)得到的三個典型時刻的流場壓力分布Fig.8 Pressure distribution of the flow field at three typical time instants,obtained by using artificial viscous SPH method (left) and the present Riemann SPH method (right),respectively

圖8 分別采用人工黏性SPH 方法(左)和當前的黎曼SPH 方法(右)得到的三個典型時刻的流場壓力分布(續(xù))Fig.8 Pressure distribution of the flow field at three typical time instants,obtained by using artificial viscous SPH method (left) and the present Riemann SPH method (right),respectively (continued)

為了定量地比較兩種方法,圖9 給出了得到的兩個測點的壓力時歷對比,并與Geers-Hunter 模型[40](簡稱GH 模型)得到的解進行對比.可以發(fā)現(xiàn),兩種SPH 方法得到的曲線在波陣面附近存在明顯的壓力爬升期,且在峰值后呈現(xiàn)不同程度的壓力振蕩.但相對而言,黎曼SPH 的爬升期更小,壓力振蕩更弱一些,更重要的是壓力峰值精度明顯提升,與GH 模型吻合良好.

圖9 不同方法得到的測點壓力時歷Fig.9 Pressure histories of measuring points obtained by different methods

3.2 水下爆炸氣泡模擬

本節(jié)對自由場中的水下爆炸氣泡進行模擬.該水下爆炸實驗由哈爾濱工程大學在江西實驗中心的小型爆炸水箱中完成.實驗中1 g 的黑索金炸藥在1.7 m 水深處引爆.由文獻[20]可確定氣泡相應的初始條件為r0=0.02 m,初始壓力為122 MPa.對于水面上方的空氣,其初始壓力設置為101325 Pa.計算域的初始粒子間距取為r0/8.考慮到氣泡脈動過程中氣體粒子體積變化劇烈,導致與周圍粒子水粒子的分辨率差異十分顯著,引起SPH 數(shù)值精度大大降低,甚至會導致計算崩潰[41],因此,本文采用了Sun 等[13]提出的自適應粒子分割與合并算法,該算法以粒子的體積和粒子之間的距離作為判據(jù),對滿足預設條件的水粒子和氣泡粒子進行自動的分割與合并,以確保粒子的體積變化在一個較小的范圍內(nèi).此外,為了維持整個系統(tǒng)的守恒性,依據(jù)質(zhì)量、動量和能量守恒,得到粒子分割前后的物理量之間的關系.通過上述處理,有效地克服了傳統(tǒng)SPH方法在模擬高壓氣泡時存在的精度低、穩(wěn)定性差的難題.

圖10 顯示了實驗和SPH 得到的幾個典型時刻的壓力場.可以看出,在內(nèi)外壓差的作用下,氣泡迅速發(fā)生膨脹.在膨脹初期,向外輻射沖擊波.在第一個膨脹過程中氣泡呈球形,但在氣泡收縮過程中,由于重力的作用,氣泡下方壓力要大于上方壓力,導致氣泡下半部分的收縮速度要大于上半部分,因此氣泡上浮,同時產(chǎn)生向上的氣泡射流.此過程中氣泡不再呈圓形.在35.5 ms 左右,射流穿透氣泡上表面形成環(huán)狀氣泡,此時氣泡體積最小.由于高速射流的穿透,氣泡在垂向上被拉長,在回彈過程中呈梨形,如圖10(f)所示.由圖10 可以看出,SPH 獲得的氣泡運動過程與實驗結(jié)果基本吻合,進一步驗證了模型的有效性和精度.

圖10 水下爆炸實驗與SPH 數(shù)值結(jié)果對比Fig.10 Comparison between underwater explosion experiment and SPH numerical results

圖11 給出了SPH 得到的氣泡的等效半徑時歷曲線,其中等效半徑req由下式計算得到

圖11 水下爆炸氣泡等效半徑時歷Fig.11 Time history of equivalent radius of the underwater explosion bubble

式中Vb代表氣泡粒子的體積.可以看出,氣泡半徑經(jīng)歷了先增大再減小然后再變大的歷程.在第一個脈動過程中,氣泡最大等效半徑和脈動周期分別為0.204 m 和35.34 ms.將SPH 得到的曲線與實驗結(jié)果以及Keller 模型解[42]進行對比可以發(fā)現(xiàn),與實驗數(shù)據(jù)相比,SPH 結(jié)果得到的氣泡最大等效半徑略小,但周期吻合非常好.SPH 的結(jié)果與Keller 模型解相比,氣泡最大等效半徑幾乎完全相同,但氣泡脈動周期略微偏大.

3.3 水下爆炸氣泡與圓板的相互作用

在近場水下爆炸過程中,結(jié)構(gòu)不僅會受到?jīng)_擊波載荷的作用,同時氣泡的脈動和射流也會造成結(jié)構(gòu)的變形和毀傷.目前已有的水下爆炸SPH 模擬中僅考慮前期沖擊波載荷的作用,而由于計算精度和效率的限制,對于后期的氣泡脈動和射流對結(jié)構(gòu)的作用未見相關研究.基于本文建立的SPH-RKPM 流固耦合方法,本節(jié)進一步模擬了水下爆炸高壓氣泡對背空平板的作用.

該算例的計算模型如下:一個初始半徑為0.025 m 的球形氣泡位于水深0.3 m 處,其初始壓力為95 MPa.在水面上有一個半徑為0.5 m 的圓板,其上方為空氣,圓板的厚度取為0.01 m,四周剛性固定.重力大小為g=9.8 m/s2,方向豎直向下.圓板采用Q235 鋼材,密度 為 ρ0=7850 kg/m3,楊氏模量E=211 GPa,泊松比 ν=0.3,采用的是Johnson-Cook 強化彈塑性模型,具體參數(shù)見文獻[43].計算域的粒子分辨率取為5 mm.

模擬得到的幾個典型時刻的流場粒子分布與壓力分布如圖12 所示,相應時刻的板的Mises 應力和等效塑性應變分布如圖13 所示.由于氣泡初始壓力遠大于周圍的水,氣泡迅速膨脹并向外輻射強烈的沖擊波,沖擊波向四周傳播,到達平板后發(fā)生反射.沖擊波的作用使圓板出現(xiàn)了明顯的應力,如圖13(a)所示.可以發(fā)現(xiàn),當沖擊波反射過后,氣泡膨脹仍然對周圍的水有擠壓作用,導致平板與氣泡中間形成一個高壓力區(qū),如圖12(b)所示,在此高壓區(qū)的作用下,平板的Mises 應力不斷增大,在2 ms 左右其幅值達到約256 MPa,板的中心位置發(fā)生了塑性變形.隨著氣泡的膨脹,泡內(nèi)壓力迅速降低,小于周圍水的壓力.由于運動慣性,氣泡仍繼續(xù)膨脹,但速度明顯降低,因此周圍水受到的擾動變小,壓力也變小,因而板的Mises 應力值減小,但等效塑性應變分布保持不變.在約15 ms 時,氣泡體積達到最大值,此時氣泡的膨脹速度為0,在周圍流場擠壓下開始發(fā)生收縮.由于氣泡下方水的壓力大于氣泡上方,因此氣泡下方的坍塌速度也要高于上方,生成向上的高速射流.在約45.5 ms 時,射流穿透氣泡上表面,形成中空的環(huán)形氣泡,如圖12(e)所示.當氣泡射流到達板下方時,會導致圓板中心區(qū)域產(chǎn)生較大的流場載荷,約為8 MPa,如圖12(f)所示.在氣泡射流的沖擊下,圓板中心的Mises 應力顯著增大,最大約為261 MPa,同時等效塑性應變也略有增大,如圖13(f)所示.

圖12 幾個典型時刻的流場壓力(左)和粒子分布(右)Fig.12 Flow field pressure (left) and particle distributions (right) at several typical time instants

圖13 幾個典型時刻圓板的Mises 應力和等效塑性應變分布Fig.13 Mises stress and equivalent plastic strain distributions of the circular plate at several typical time instants

圖14 定量給出了該算例中平板中心點的Mises應力以及位移的時歷曲線.可以看到,與圖13 一致,在模擬初期沖擊波作用下圓板中心點應力迅速增大,峰值約為253 MPa,且隨著應力波的反射,中心點應力出現(xiàn)較為劇烈的振蕩.在氣泡膨脹過程中,中心點應力在100 MPa 呈緩慢下降趨勢,在氣泡收縮的后期(43～46 ms),圓板中心點的應力迅速降低.隨著氣泡射流的產(chǎn)生,中心點的應力峰值再次迅速增大,其峰值約為261 MPa.由位移時歷曲線可以發(fā)現(xiàn),整個模擬過程中位移主要有兩個峰值,其大小分別約為0.01165 m 和0.0169 m.第一個峰值的產(chǎn)生與前期的沖擊波載荷有關,而第二個峰值由氣泡射流載荷引起.可以發(fā)現(xiàn),氣泡射流導致的結(jié)構(gòu)應力和變形均大于前期沖擊波.

圖14 圓板中心點的Mises 應力和位移時歷曲線Fig.14 Time histories of the Mises stress and displacement of the center point of the circular plate

3.4 水下爆炸對舷側(cè)結(jié)構(gòu)的毀傷模擬

在上一節(jié)采用SPH-RKPM 流固耦合模型對爆炸氣泡與簡單平板的耦合進行了模擬,本節(jié)對接觸爆炸作用下的艦船舷側(cè)結(jié)構(gòu)的毀傷進行研究.計算模型的設置如圖15 所示,整個結(jié)構(gòu)的殼體結(jié)構(gòu)厚度設置為8 mm,材料采用945 鋼,材料的初始屈服應力設置為 σy0=440 MPa.本文采用了線性硬化模型,硬化模量Ep設置為1000 MPa.對于本文中的水下爆炸毀傷算例,材料的應變率效應不能忽略,其采用Cowper-Symonds 模型[30,45]描述,此時動態(tài)屈服應力σyd計算為

圖15 艦船舷側(cè)防護結(jié)構(gòu)受接觸爆炸作用示意圖,灰色粒子為裝藥Fig.15 Schematic diagram of the ship side protective structure subjected to contact explosion,with grey particles as charge

式中,r代表各向同性強化的狀態(tài)變量,是塑性應變率.材料常數(shù)設為C=9870,p=2.43.材料損傷演化由Lemaitre 損傷模型[44]得到,損傷判據(jù)參數(shù)設為Dc=0.32,臨界損傷應變=0.16,斷裂應變=0.28.計算中采用均勻的粒子分布,粒子間距設置為0.02 m,TNT 當量為10 kg,流體和結(jié)構(gòu)的時間步長均設置為1.0 × 10-6s.

圖16 給出了0.25 ms 時的水中壓力云圖、爆炸氣體產(chǎn)物的形態(tài)以及結(jié)構(gòu)的變形,在結(jié)構(gòu)的變形構(gòu)型上繪制出了殼結(jié)構(gòu)厚度方向最外層積分點上的Mises 等效應力云圖.從圖中可以看出,在計算開始之后,隨著爆炸產(chǎn)物的膨脹,與裝藥接觸的船體板立刻產(chǎn)生了局部的凹陷大變形,并在炸藥接觸的中心部分結(jié)構(gòu)發(fā)生了破裂,同時加強筋和面板的交叉部位產(chǎn)生了應力集中現(xiàn)象.

圖16 水下接觸爆炸作用下艦船舷側(cè)防護結(jié)構(gòu)在 t=0.25 ms 時的流場壓力云圖以及結(jié)構(gòu)的Mises 應力云圖Fig.16 Pressure fields of the flow field and Mises stress fields of the warship side protective structure under underwater contact explosion att=0.25 ms

圖17 給出了幾個典型時刻,水下接觸爆炸作用下艦船舷側(cè)防護結(jié)構(gòu)的變形構(gòu)型以及應力云圖,在每個時刻的應力云圖的右側(cè)都單獨顯示了弧形吸能板的變形.從圖中可以看出,t=0.5 ms 時在爆炸作用下與裝藥緊挨著的面板發(fā)生了局部凹陷和局部斷裂,同時應力波傳播到了內(nèi)層艙壁上,之后在t=1.0 ms 時結(jié)構(gòu)變形加劇,并且外層艙壁產(chǎn)生的斷裂破片碰撞到了內(nèi)部的艙壁上,內(nèi)部弧形吸能板產(chǎn)生了初步的變形,之后防護結(jié)構(gòu)進一步變形,最外側(cè)兩層防護艙壁均產(chǎn)生了較大的斷裂破口,并與內(nèi)部艙壁發(fā)生了接觸碰撞,到t=4.0 ms 時結(jié)構(gòu)的斷裂破口趨于穩(wěn)定,最靠近爆點的弧形吸能板被完全壓潰,產(chǎn)生了動力屈曲.

圖17 水下接觸爆炸載荷作用下雙層底幾個典型時刻的Mises 應力云圖Fig.17 Mises stress fields at several typical time instants of the double bottom under underwater contact explosion loading

3.5 水下爆炸對潛艇結(jié)構(gòu)的毀傷模擬

最后,基于本文建立的SPH-RKPM 流固耦合模型,模擬了近場水下爆炸對潛艇結(jié)構(gòu)的毀傷.計算模型縱剖面如圖18 所示,潛艇長度為110.5 m,耐壓殼半徑為5.6 m.在潛艇中部附近下方水域中有一個半徑為0.65 m 的球形TNT 藥包,藥包中心距離潛艇外殼1.0 m.潛艇所用材料及參數(shù)與上節(jié)相同.整個模型采用均勻的粒子分布,粒子間距設置為0.2 m,流體和結(jié)構(gòu)的時間步長均設置為2.0 × 10-7s.

圖18 水下爆炸對潛艇結(jié)構(gòu)的毀傷計算模型Fig.18 Computational model of submarine structure damaged by underwater explosion

圖19 給出了幾個典型時刻的流場壓力分布和結(jié)構(gòu)的Mises 應力分布.可以發(fā)現(xiàn)炸藥起爆后變?yōu)楦邷馗邏簹怏w,在水中形成強烈沖擊波,在沖擊波和爆轟氣體共同作用下,潛艇外板迅速產(chǎn)生變形屈服并形成局部凹陷和破口.隨著爆轟氣體的膨脹和沖擊波的傳播,破口半徑不斷增大,屈服區(qū)域迅速向四周擴展,但由于徑向上艙段間的環(huán)形肋板的存在,相鄰艙段的應力和變形明顯降低,減緩了破口的增大,說明其可以有效地增強潛艇的抗水下爆炸毀傷能力.

圖19 水下爆炸對潛艇結(jié)構(gòu)的毀傷SPH-RKPM 數(shù)值結(jié)果Fig.19 SPH-RKPM numerical results of submarine structure damage caused by underwater explosion

4 結(jié)論

本文針對近場水下爆炸流固耦合問題,在前期研究的基礎上,利用無網(wǎng)格粒子法的天然拉格朗日優(yōu)勢,開發(fā)了適用于間斷流場求解的黎曼SPH 方法和結(jié)構(gòu)損傷斷裂模擬的RKPM 方法,并應用法向通量邊界建立了高精度SPH-RKPM 瞬態(tài)流固耦合模型.首先通過幾個典型算例驗證了上述算法的精度,在此基礎上,對水下爆炸沖擊波傳播、氣泡脈動和射流以及結(jié)構(gòu)毀傷等現(xiàn)象進行了數(shù)值模擬,得到以下結(jié)論.

(1)相對于傳統(tǒng)人工黏性SPH 方法,本文建立的二階黎曼SPH 方法在模擬水下爆炸沖擊波時,在多相界面處可以得到更加連續(xù)的壓力場,以及更加尖銳的波陣面和更加準確的壓力峰值;通過引入粒子自適應分割和合并算法,能有效地克服氣泡脈動過程粒子體積變化導致的精度下降問題,從而實現(xiàn)水下爆炸氣泡的動力學特性的高精度求解.

(2)基于RKPM 方法,采用退化實體幾何表述,并應用Lemaitre 損傷模型,開發(fā)了基于粒子法的裂紋萌生拓展算法,可以較好地模擬殼結(jié)構(gòu)的大變形動力屈曲以及損傷斷裂行為.

(3)應用法向通量邊界條件和SPH 插值,建立了SPH 流體和RKPM 單層殼結(jié)構(gòu)的瞬態(tài)強非線性流固耦合模型,該模型對于高速水射流和水下爆炸等瞬態(tài)沖擊問題具有良好的數(shù)值精度.

(4)依據(jù)SPH-RKPM 流固耦合模型模擬了水下爆炸氣泡與圓板的相互作用,氣泡在脈動后期產(chǎn)生了指向圓板的高速射流載荷,前期的沖擊波載荷和后期的射流載荷均會對結(jié)構(gòu)造成較大的沖擊;此外模擬了水下接觸爆炸對艦船舷側(cè)防護結(jié)構(gòu)以及潛艇結(jié)構(gòu)的毀傷,給出了結(jié)構(gòu)的變形和毀傷特性.

致謝

感謝哈爾濱工程大學博士研究生薛冰在論文撰寫中提供的幫助.