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

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

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

引言

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

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

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

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

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

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

1 基本理論

1.1 黎曼SPH 間斷流場(chǎng)模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

其中,矢量rji計(jì)算為rji=xj-xi.

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

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

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

1.3 流固耦合算法

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

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

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

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

2 數(shù)值驗(yàn)證

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

2.1 激波管

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

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

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

圖1 黎曼SPH 方法得到的“激波-密度波相互作用問(wèn)題”在t=1.8時(shí)刻的密度和速度分布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 圓管受壓屈曲

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

圖2 圓管受壓動(dòng)力屈曲計(jì)算模型Fig.2 Computational model of the dynamic buckling of the circular tube under compression

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

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

圖3 圓管受壓動(dòng)力屈曲過(guò)程Fig.3 Dynamic buckling process of the circular tube under compression

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

圖4 圓管屈曲變形數(shù)值結(jié)果(上)與文獻(xiàn)[37]中的實(shí)驗(yàn)結(jié)果(下)對(duì)比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é)中,通過(guò)模擬高速水射流沖擊彈性平板算例,以驗(yàn)證當(dāng)前采用的流固耦合算法在強(qiáng)沖擊問(wèn)題中的精度.其計(jì)算模型如圖5 所示,一個(gè)半徑為1.5 mm 的球頭柱形水射流,其長(zhǎng)度為12 mm,以570 m/s 的初始速度沖擊一個(gè)四周剛固的方形板中心,方形板的長(zhǎng)度和厚度分別是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 水射流沖擊彈性板計(jì)算模型Fig.5 Computational model of the water jet impacting elastic plate

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

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

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

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

圖8 分別采用人工黏性SPH 方法(左)和當(dāng)前的黎曼SPH 方法(右)得到的三個(gè)典型時(shí)刻的流場(chǎng)壓力分布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 方法(左)和當(dāng)前的黎曼SPH 方法(右)得到的三個(gè)典型時(shí)刻的流場(chǎng)壓力分布(續(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 給出了得到的兩個(gè)測(cè)點(diǎn)的壓力時(shí)歷對(duì)比,并與Geers-Hunter 模型[40](簡(jiǎn)稱GH 模型)得到的解進(jìn)行對(duì)比.可以發(fā)現(xiàn),兩種SPH 方法得到的曲線在波陣面附近存在明顯的壓力爬升期,且在峰值后呈現(xiàn)不同程度的壓力振蕩.但相對(duì)而言,黎曼SPH 的爬升期更小,壓力振蕩更弱一些,更重要的是壓力峰值精度明顯提升,與GH 模型吻合良好.

圖9 不同方法得到的測(cè)點(diǎn)壓力時(shí)歷Fig.9 Pressure histories of measuring points obtained by different methods

3.2 水下爆炸氣泡模擬

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

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

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

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

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

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

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

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

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

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

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

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

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

圖14 圓板中心點(diǎn)的Mises 應(yīng)力和位移時(shí)歷曲線Fig.14 Time histories of the Mises stress and displacement of the center point of the circular plate

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

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

致謝

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