遠(yuǎn)場(chǎng)水下爆炸作用下平板的沖擊響應(yīng)仿真*

辛春亮,王新泉,涂 建,王 偉,葉志萍

(1 北京航天長(zhǎng)征飛行器研究所,北京 100076; 2中國(guó)航天員科研訓(xùn)練中心,北京 100094)

遠(yuǎn)場(chǎng)水下爆炸作用下平板的沖擊響應(yīng)仿真*

辛春亮1,王新泉1,涂 建1,王 偉1,葉志萍2

(1 北京航天長(zhǎng)征飛行器研究所,北京 100076; 2中國(guó)航天員科研訓(xùn)練中心,北京 100094)

為了對(duì)艦船水下防護(hù)結(jié)構(gòu)設(shè)計(jì)提供參考,采用LS-DYNA軟件中的聲固耦合法對(duì)遠(yuǎn)場(chǎng)水下爆炸作用下平板的沖擊響應(yīng)問(wèn)題進(jìn)行了模擬,計(jì)算結(jié)果與公認(rèn)的雙漸進(jìn)近似算法計(jì)算結(jié)果較為接近。計(jì)算模型中采用了空化聲單元和粘性非反射邊界,并對(duì)比了網(wǎng)格加密、不同水柱高度和是否考慮空化的影響。LS-DYNA的聲固耦合法計(jì)算速度快,計(jì)算準(zhǔn)確度接近公認(rèn)的雙漸進(jìn)近似算法,為國(guó)內(nèi)水下爆炸數(shù)值計(jì)算人員提供了一種新的思路和計(jì)算工具。

LS-DYNA軟件;遠(yuǎn)場(chǎng)水下爆炸;粘性邊界條件;數(shù)值模擬

0 引言

艦船在水下爆炸下的結(jié)構(gòu)響應(yīng)一直是艦船設(shè)計(jì)者最為關(guān)心的問(wèn)題,為此國(guó)內(nèi)外進(jìn)行了大量實(shí)船水下爆炸試驗(yàn)。水下爆炸試驗(yàn)耗費(fèi)大、組織時(shí)間長(zhǎng)、限制多,還會(huì)帶來(lái)很多環(huán)境問(wèn)題。隨著計(jì)算技術(shù)的發(fā)展,研究者們逐漸轉(zhuǎn)向采用數(shù)值計(jì)算來(lái)解決實(shí)船水下爆炸問(wèn)題。

1 遠(yuǎn)場(chǎng)水下爆炸數(shù)值計(jì)算方法

根據(jù)艦船結(jié)構(gòu)的受力狀態(tài),水下爆炸可分為近場(chǎng)和遠(yuǎn)場(chǎng)水下爆炸。在遠(yuǎn)場(chǎng)水下爆炸中艦船結(jié)構(gòu)只發(fā)生彈性變形,氣泡效應(yīng)通?？梢员缓雎?數(shù)值計(jì)算模型中可只考慮沖擊波和空化加載。LS-DYNA是LSTC公司推出的通用顯式動(dòng)力分析程序,該軟件有多種算法可用于水下爆炸數(shù)值模擬,如流固耦合法、LOAD_SSA、LS-DYNA/USA等[1]。

流固耦合法主要用于近場(chǎng)水下爆炸計(jì)算。由于計(jì)算耗費(fèi)和計(jì)算準(zhǔn)確度的原因,難以用于遠(yuǎn)場(chǎng)水下爆炸模擬。這種方法需要將從爆源至結(jié)構(gòu)、結(jié)構(gòu)附近很大范圍的流體劃分網(wǎng)格,爆源附近的網(wǎng)格需要尤其細(xì)密,否則會(huì)抹平?jīng)_擊波峰值,如此一來(lái),計(jì)算模型將極其龐大。此外,EULER和ALE算法本身的耗散效應(yīng)也會(huì)抹平?jīng)_擊波峰值,使之傳播失真。

LOAD_SSA可以用來(lái)模擬遠(yuǎn)場(chǎng)水下爆炸對(duì)結(jié)構(gòu)的響應(yīng)。這種方法可以考慮球面入射波、反射波、輻射波以及附帶水質(zhì)量帶來(lái)的影響。由于不需要建立流體網(wǎng)格,這種方法計(jì)算速度很快,缺點(diǎn)是附帶水質(zhì)量沒(méi)有考慮主沖擊波過(guò)后后續(xù)流體的空化特性。

目前流體動(dòng)力學(xué)計(jì)算軟件更多采用有限元和非反射邊界來(lái)模擬遠(yuǎn)場(chǎng)水下爆炸問(wèn)題。這種方法可稱之為聲固耦合法,流體為線性聲學(xué)介質(zhì),采用有限元模擬,在非反射邊界處波可以無(wú)反射地透射出去。例如常用的有限元和邊界元耦合法,在LS-DYNA中稱之為L(zhǎng)S-DYNA/USA模塊,這種方法在流體邊界上采用雙漸進(jìn)近似算法(DAA),可以考慮附帶水質(zhì)量的影響,計(jì)算準(zhǔn)確度很高,但這種方法計(jì)算附帶水質(zhì)量矩陣時(shí),需要不斷對(duì)非稀疏矩陣求逆,計(jì)算耗費(fèi)較大,計(jì)算時(shí)間較長(zhǎng)。目前國(guó)內(nèi)研究者仍然難以獲得LS-DYNA/USA模塊。

LS-DYNA中另一種非反射邊界為粘性邊界[2],這是一種局部近似非反射邊界,最早由Lysmer和Kuhlemeyer[3]提出用來(lái)解決土壤和結(jié)構(gòu)相互作用問(wèn)題,后來(lái)Cohen和Jennings[4]又做了進(jìn)一步的研究。對(duì)于入射壓力波,粘性邊界算法計(jì)算每個(gè)邊界節(jié)點(diǎn)的法向應(yīng)力和切向應(yīng)力,在每個(gè)節(jié)點(diǎn)匹配入射應(yīng)力,由此得到的反射應(yīng)力為零。

σnorm=-ρcdvnorm

τshear=-ρcsvtan

(1)

式中:σnorm是邊界法向應(yīng)力;τshear是邊界切向應(yīng)力;ρ是材料密度;cd是畸變波速;cs是膨脹波速;v是質(zhì)點(diǎn)速度。

為了模擬流體中的空化效應(yīng),水域采用空化聲單元(CAFE)。空化聲單元將流體域看做帶有雙線性狀態(tài)方程的聲學(xué)域,假定密度變化很小,只有畸變波可以在流體中傳播,這就意味著當(dāng)粘性邊界用于流體邊界時(shí)只需匹配σnorm。在流體中忽略粘性,流體的應(yīng)力狀態(tài)可以描述為流體的總壓力p,這樣當(dāng)粘性邊界施加在流體網(wǎng)格邊界時(shí),壓力可采用如下簡(jiǎn)化形式:

p=-ρcvnorm

(2)

式中c是流體聲速。

粘性邊界可較為精確的處理法向入射波[3],但當(dāng)入射角度從90°降到0°時(shí)精度急劇下降。此外,這種邊界不考慮附帶水質(zhì)量,在附帶水質(zhì)量顯著影響艦船響應(yīng)的問(wèn)題中計(jì)算精度不如雙漸進(jìn)近似算法。然而式(2)采用平面波假設(shè)形式,平面波假設(shè)均施加于水下爆炸早期階段,因此對(duì)于空化聲單元模型早期響應(yīng),粘性邊界與雙漸進(jìn)近似算法的計(jì)算結(jié)果差別不會(huì)很大,但計(jì)算速度要快得多。

2 Bleich-Sandler平板問(wèn)題

Bleich-Sandler平板問(wèn)題是一個(gè)非常經(jīng)典的考慮空化效應(yīng)的遠(yuǎn)場(chǎng)水下爆炸驗(yàn)證算例[5-6]。在該問(wèn)題中,平板浮于3.81 m高水柱的自由面上,平板和水域只允許垂向(Z向)運(yùn)動(dòng)。模型尺寸、材料參數(shù)見(jiàn)表1、表2和圖1。

采用8節(jié)點(diǎn)聲單元在水柱X、Y、Z三個(gè)方向分別劃分100、1、1個(gè)網(wǎng)格,即流體網(wǎng)格尺寸為邊長(zhǎng)0.038 1 m的正方體,采用聲學(xué)材料MAT_90來(lái)模擬流體。為了模擬無(wú)限水域,水域底部施加粘性邊界條件。采用一個(gè)Belyschko-Tsay殼單元來(lái)模擬平板,材料模型為線彈性。

表2 平板幾何尺寸和材料參數(shù)

圖1 Bleich-Sandler問(wèn)題示意圖

在LS-DYNA聲固耦合法中,耦合界面處的流體和結(jié)構(gòu)網(wǎng)格可以采用3種方式:1)共節(jié)點(diǎn);2)網(wǎng)格一一對(duì)應(yīng)耦合;3)網(wǎng)格不匹配耦合。這里采用第二種一一對(duì)應(yīng)的耦合方式。為了計(jì)算穩(wěn)定,聲固耦合時(shí),流體單元尺寸須滿足:

(3)

式中:ρ是流體密度;D是流體聲單元在垂直于結(jié)構(gòu)濕面方向上的厚度;ρs是結(jié)構(gòu)密度;ts是結(jié)構(gòu)單元厚度。

輸入載荷為以指數(shù)形式衰減的平面入射沖擊波,見(jiàn)式(4),輸入點(diǎn)沖擊波參數(shù)如表3所列。零時(shí)刻計(jì)算程序根據(jù)輸入載荷和式(4)初始化流體壓力,沖擊波波前位于距平板濕面一個(gè)流體單元的節(jié)點(diǎn)處。這樣做的優(yōu)點(diǎn)是不用計(jì)算入射波在流體中的傳播,避免了沖擊波傳播失真和能量耗散。

表3 Bleich-Sandler問(wèn)題中入射波參數(shù)

(4)

式中:p0是輸入點(diǎn)沖擊波峰值超壓;τ是時(shí)間常數(shù)。

3 粘性邊界算法與雙漸進(jìn)近似算法計(jì)算結(jié)果對(duì)比

圖2是平板垂向速度計(jì)算結(jié)果,圖中還疊加了采用雙漸進(jìn)近似邊界[2]時(shí)的計(jì)算曲線,該計(jì)算結(jié)果已經(jīng)過(guò)驗(yàn)證。與雙漸進(jìn)近似邊界條件相比,采用粘性邊界條件時(shí),平板最大速度略低,空化加載前負(fù)向最大速度也存在差異,除此以外兩條曲線其他部分吻合較好。產(chǎn)生這些差異的原因是雙漸進(jìn)近似算法邊界條件考慮了附帶水質(zhì)量而粘性邊界條件未考慮。

圖2 分別采用雙漸進(jìn)近似和粘性邊界條件時(shí)平板垂向速度

圖2中,采用粘性邊界條件時(shí)計(jì)算曲線峰值過(guò)后存在輕微偽震蕩,偽震蕩與在有限元模型考慮空化效應(yīng)有關(guān)。圖3是不同網(wǎng)格密度的計(jì)算結(jié)果,可以看出網(wǎng)格加密后曲線更為光滑,噪音很低,而峰值幾乎沒(méi)有變化,但空化加載時(shí)間提前。網(wǎng)格加密到2和4倍原始網(wǎng)格后,曲線B和C幾乎重合,這說(shuō)明,網(wǎng)格加密到一定程度后,計(jì)算結(jié)果趨于收斂。

圖3 采用不同網(wǎng)格密度時(shí)平板垂向速度

在遠(yuǎn)場(chǎng)水下爆炸問(wèn)題中,船體的整體運(yùn)動(dòng)受附帶水質(zhì)量的影響。由于非反射邊界不考慮附帶水質(zhì)量效應(yīng),這就需要流體域足夠大,能夠包含全部附帶水質(zhì)量。圖4是計(jì)算模型中考慮不同水柱高度時(shí)的計(jì)算結(jié)果。在這些計(jì)算模型中,流體長(zhǎng)度方向網(wǎng)格數(shù)量分別為50、100、200和400,即流體網(wǎng)格尺寸不變。由圖可見(jiàn),水柱長(zhǎng)度越小,與雙漸進(jìn)近似算法計(jì)算結(jié)果偏差就越大。曲線C和D對(duì)應(yīng)的水柱長(zhǎng)度分別為7.62 m和15.24 m,這兩條曲線幾乎完全重合,與雙漸進(jìn)近似算法計(jì)算曲線也非常接近。這說(shuō)明,隨著水柱長(zhǎng)度的加大,粘性邊界條件不考慮附帶水質(zhì)量的缺陷逐漸被抵消,計(jì)算結(jié)果趨于收斂。

圖4 不同水柱高度時(shí)平板垂向速度

空化加載對(duì)艦船結(jié)構(gòu)影響較大,如圖5所示。由于平板最大速度由主沖擊波決定,所以是否考慮空化對(duì)此沒(méi)有影響。而不考慮空化時(shí)平板速度在峰值過(guò)后很快衰減為零,且沒(méi)有負(fù)向速度。

圖5 考慮和不考慮空化時(shí)平板垂向速度

圖6是水域空化時(shí)空?qǐng)D,對(duì)應(yīng)計(jì)算模型的水柱長(zhǎng)

圖6 水域空化時(shí)空?qǐng)D

度為7.62 m、流體長(zhǎng)度方向網(wǎng)格數(shù)量為400。空化形成和閉合時(shí)間大約在0.4 ms和10.1 ms。圖6(a)形狀與圖6(b)文獻(xiàn)[6-7]計(jì)算結(jié)果大致吻合。圖6(a)中,空化時(shí)空區(qū)上邊界非常清晰,由于氣泡和水混雜在一起,下邊界較難明確分辨。在空化出現(xiàn)的時(shí)間段內(nèi),空化并未延伸到平板底部,平板和空化帶之間始終存在未空化水的堆積區(qū)。空化帶上、下邊界收斂閉合在一起時(shí)空化潰滅,對(duì)應(yīng)的時(shí)空位置為(10.1 ms,2.2 m)。圖7是距離平板底部2 m處水單元壓力曲線,曲線中第一個(gè)壓力峰值是入射波將要抵達(dá)平板底部時(shí)該單元的壓力,第二個(gè)峰值是衰減后的入射波和平板反射波的疊加形成的,隨后由于平板運(yùn)動(dòng)導(dǎo)致水中壓力急劇下降為空化壓力,9.9 ms時(shí)刻由于空化閉合,壓力又迅速回升。

圖7 距平板底部2 m處水單元壓力曲線

4 結(jié)論

文中介紹了LS-DYNA中的聲固耦合法,并采用該方法對(duì)Bleich-Sandler平板問(wèn)題進(jìn)行了模擬。聲固耦合法采用粘性非反射邊界,主要考慮沖擊波和空化加載,不計(jì)及氣泡脈動(dòng)和附帶水質(zhì)量。當(dāng)計(jì)算水域取值較大時(shí),LS-DYNA中的聲固耦合法能夠較為準(zhǔn)確的預(yù)測(cè)遠(yuǎn)場(chǎng)水下爆炸下平板的響應(yīng)。聲固耦合法內(nèi)嵌于LS-DYNA軟件中,適合于遠(yuǎn)場(chǎng)水下爆炸分析,計(jì)算效率高,這種分析方法可為國(guó)內(nèi)水下爆炸仿真人員和艦船設(shè)計(jì)者提供參考。

[1] 辛春亮,秦健,劉科種,等.基于LS-DYNA軟件的水下爆炸數(shù)值模擬研究 [J].彈箭與制導(dǎo)學(xué)報(bào),2008,28(3):156-158.

[2] KLENOW B,BROWN A.Assessment of non-reflecting boundary conditions for application in far-field UNDEX finite element models[C]// 77th Shock and Vibration Symposium.Monterey:SAVIAC,2006.

[3] LYSMER J,KUHLEMEYER R L.Finite dynamic model for infinite media [J].Journal of Engineering Mechanics Division,1969,95(4):859-878.

[4] CHOEN M.Silent boundary methods for transient wave analysis:PB82-201831[R].1983.

[5] MAIR H U.Benchmarks for submerged structure response to underwater explosions[J].Shock and Vibration,1999,6(4):169-181.

[6] BLEICH H H,SANDLER I S.Interaction between structures and bilinear fluids[J].International Journal of Solids and Structures.1970,6(5):617-639.

[7] SPRAGUE M A,GEERS T L.Computational treatment of cavitation effects in near-free-surface underwater shock analysis [J].Shock and Vibration,2001,8(2):105-122.

NumericalSimulationofPlateResponseSubjecttoFar-fieldUnderwaterExplosion

XIN Chunliang1,WANG Xinquan1,TU Jian1,WANG Wei1,YE Zhiping2

(1 Beijing Institute of Space Long March Vehicle,Beijing 100076,China; 2 China Astronaut Research and Training Center,Beijing 100094,China)

Acoustic-structure interaction method in LS-DYNA software was used to solve plate response subject to far-field underwater explosion.Cavitation acoustic finite element (CAFE) and non-reflecting boundary were used in UNDEX models.Results with different fluid mesh density,different length of fluid column,with or without CAFE effects were compared with DAA results.This kind of method could offer significant savings in computational time over DAA without greatly affecting accuracy of results.

LS-DYNA software; far-field underwater explosion; viscous boundary condition; numerical simulation

10.15892/j.cnki.djzdxb.2017.02.019

2016-06-12

辛春亮(1973-),男,山東日照人,研究員,博士,研究方向:戰(zhàn)斗部設(shè)計(jì)和數(shù)值模擬。

TG156

A