李上明 屈明

摘要：

基于ALE方法分析楔形體入水問(wèn)題，楔形體采用拉格朗日網(wǎng)格離散，空氣和水采用ALE網(wǎng)格離散。將楔形體視為剛體，空氣和水的力學(xué)行為分別采用Gamma定律和GRUNEISEN狀態(tài)方程模擬。討論ALE流固耦合關(guān)鍵字中罰函數(shù)罰因子的取值方法，提出相應(yīng)的建議原則；分析水域截?cái)噙吔鐚?duì)楔形體響應(yīng)的影響，給出模擬無(wú)限水域的截?cái)噙吔缥恢玫慕ㄗh值；分析楔形體表面壓力振蕩的原因，提出楔形體表面壓力獲取方法。通過(guò)試驗(yàn)結(jié)果對(duì)比分析，驗(yàn)證方法的合理性。

關(guān)鍵詞：

入水； ALE方法； 截?cái)噙吔纾?流固耦合； 無(wú)限水域； 罰函數(shù)

中圖分類號(hào)： O347.5

文獻(xiàn)標(biāo)志碼： B

Wedge water entry analysis based on ALE method

LI Shangming， QU Ming

（Institute of Systems Engineering， China Academy of Engineering Physics， Mianyang 621999， Sichuan， China）

Abstract：

The wedge water entry problem is analyzed based on ALE method. The wedge is dispersed using Lagrange meshes. The air and water are dispersed using ALE meshes. The wedge is considered as a rigid body. The mechanics behaviors of air and water are simulated by Gamma law and GRUNEISEN state equation respectively. The method to get the penalty factors of penalty function of ALE fluidsolid coupling keywords is discussed， then accordingly principle is proposed. The influence of the truncated boundary of water area on the wedge response is analyzed. The position of truncated boundary is proposed for modelling infinite water area. The reason of pressure vibration on wedge surface is explained， and a method is proposed to obtain the pressure on the wedge surface. The rationality of the method is verified by the comparison and analysis of test results.

Key words：

water entry； ALE method； truncated boundary； fluidsolid coupling； infinite water area； penalty function

0 引 言

隨著人們對(duì)飛機(jī)水上迫降、返回艙水上回收、空投魚(yú)雷入水等問(wèn)題的進(jìn)一步關(guān)注，越來(lái)越多的研究開(kāi)始重視結(jié)構(gòu)沖擊入水的動(dòng)態(tài)響應(yīng)特性。響應(yīng)評(píng)估的手段有理論計(jì)算、數(shù)值模擬和試驗(yàn)等。與理論計(jì)算和試驗(yàn)相比，數(shù)值模擬具有適用范圍廣、成本低等特點(diǎn)，經(jīng)常被研究人員和工程師采用。

目前，常用于結(jié)構(gòu)入水的數(shù)值模擬方法有基于VOF的CFD方法[15]、SPH方法[69]、ALE方法[1015]等。CFD方法大多集中于剛體結(jié)構(gòu)入水分析方面，其獲取的結(jié)構(gòu)入水壓力波動(dòng)較??；SPH方法適合模擬大變形問(wèn)題，常用于結(jié)構(gòu)入水分析中水域的模擬，但其邊界處理還存在較大的提升空間，且壓力波動(dòng)較大。ALE方法雖然在流固耦合面上的壓力波動(dòng)也較大，但因?yàn)槠湟呀?jīng)被集成在LSDYNA中，所以獲得廣泛的應(yīng)用。調(diào)整相關(guān)參數(shù)可以有效控制ALE壓力的波動(dòng)，提高分析的合理性。

對(duì)結(jié)構(gòu)入水流固耦合模擬問(wèn)題，LSDYNA中的ALE方法基于統(tǒng)一的連續(xù)介質(zhì)動(dòng)量守恒方程，采用類似固體之間的接觸算法模擬界面力，利用顯式有限元法統(tǒng)一計(jì)算，固體采用本構(gòu)關(guān)系模擬，流體采用狀態(tài)方程模擬，適合模擬結(jié)構(gòu)大變形和以層流特征為主的流體大位移等情況。根據(jù)ALE顯式動(dòng)力學(xué)有限元法的計(jì)算特點(diǎn)，該方法在模擬快速運(yùn)動(dòng)（即慣性效應(yīng)較強(qiáng)的情況）比較有效，但因其采用罰函數(shù)法進(jìn)行流固耦合載荷的傳遞，故流固耦合界面力會(huì)出現(xiàn)振蕩。本文通過(guò)參數(shù)分析，討論罰函數(shù)法相關(guān)罰因子的選擇方法和壓力波動(dòng)的原因，進(jìn)而確定結(jié)構(gòu)入水分析的方法。

1 流固耦合建模方法

基于某二維楔形體結(jié)構(gòu)入水問(wèn)題，討論LSDYNA流固耦合建模方法。該楔形體為等腰楔形體，底邊長(zhǎng)為1.20 m，底角為25°。楔形體頂點(diǎn)

左、右兩邊水和空氣的寬為1.85 m，水深為1.00 m，空氣高為0.40 m。楔形體視為剛體，采用平面應(yīng)變單元離散；水域和空氣域采用ALE單元離散?？諝夂退淖杂梢后w界面網(wǎng)格共節(jié)點(diǎn)處理。楔形體網(wǎng)格與空氣、水的離散網(wǎng)格重疊。楔形體附近水域最小網(wǎng)格為2 mm，相對(duì)遠(yuǎn)處區(qū)域的網(wǎng)格逐步放大，以減小網(wǎng)格數(shù)量。楔形體網(wǎng)格大小約2 mm，其網(wǎng)格離散和壓力測(cè)點(diǎn)布置見(jiàn)圖1。

1.1 物理參數(shù)

楔形體密度為467 kg/m3，彈性模量為500 MPa，泊松比為0，采用*mat_rigid模擬。水和空氣采用*mat_null和狀態(tài)方程模擬，其中：水采用GRUNEISEN狀態(tài)方程模擬，見(jiàn)式（1）；空氣通過(guò)線性多項(xiàng)式狀態(tài)方程模擬，氣體的Gamma定律狀態(tài)方程見(jiàn)式（2）。

1.2 流固耦合關(guān)鍵字

為有效進(jìn)行結(jié)構(gòu)入水流固耦合分析，通常要求在LSDYNA的輸入文件k文件中包括如下關(guān)鍵字：*ALE_MULTIMATERIAL_GROUP關(guān)鍵字，用以指定并跟蹤水和空氣材料及其邊界；*CONTROL_ALE關(guān)鍵字，用以控制ALE單元的相關(guān)算法和部分邊界處理方法；*CONSTRAINED_LAGRANGE_IN_SOLID關(guān)鍵字，用以描述流體與固體間的耦合作用，在流固耦合面上實(shí)現(xiàn)載荷傳遞。本文采用罰函數(shù)法進(jìn)行流固耦合分析。

罰函數(shù)的物理意義相當(dāng)于在結(jié)構(gòu)節(jié)點(diǎn)（從節(jié)點(diǎn)）與流體單元表面（被穿透表面、主面）之間放置一個(gè)法向彈簧，以限制結(jié)構(gòu)節(jié)點(diǎn)與流體單元面之間的穿透[13]，其耦合力

式中：K和V分別為被穿透流體單元的體積模量和體積；A為連接從節(jié)點(diǎn)的結(jié)構(gòu)單元平均面積；p為罰函數(shù)的罰因子。p用于調(diào)整流固耦合接觸剛度，其取值的大小影響流固耦合面的相互作用力。式（3）和（4）中，除p為輸入量外，其余參數(shù)均由

LSDYNA軟件自動(dòng)計(jì)算。

2 算例分析

針對(duì)某二維楔形體入水，主要討論左、右兩側(cè)截?cái)噙吔绾土P函數(shù)的罰因子對(duì)分析結(jié)果的影響。邊界條件設(shè)為：水域底部固定垂直（y向）位移，左、右兩側(cè)水域固定水平（x向）位移。初始條件為：楔形體以5.05 m/s初始速度撞擊水面，其余均為靜止。外力條件為：考慮靜水壓力和重力加速度。因該問(wèn)題具有對(duì)稱性，故在分析中楔形體只考慮y向位移，x向位移約束為0。

截?cái)噙吔缬懻摰闹饕康氖翘剿鹘財(cái)噙吔邕h(yuǎn)近對(duì)計(jì)算結(jié)果的影響，以評(píng)估模擬無(wú)限水域時(shí)截?cái)噙吔缇嚯x的選取原則。罰函數(shù)罰因子討論的主要目的是評(píng)估其大小的選取原則，以便確保模擬結(jié)果的合理性。

2.1 截?cái)噙吔缬懻?/p>

考慮楔形體頂點(diǎn)兩側(cè)水域3種不同截?cái)嚅L(zhǎng)度下的楔形體入水響應(yīng)，即：工況一為水域側(cè)面單邊長(zhǎng)為1.20 m，即楔形體半邊長(zhǎng)度的2.0倍；工況二為水域側(cè)面單邊長(zhǎng)1.85 m，即楔形體半邊長(zhǎng)度的3.1倍；工況三為水域側(cè)面單邊長(zhǎng)3.20 m，即楔形體半邊長(zhǎng)度的5.3倍。

在不同工況下楔形體的速度差見(jiàn)圖2。在0.10 s內(nèi)，水域側(cè)面長(zhǎng)從1.20 m增大到3.20 m時(shí)，最大速度差不到12 μm/s。該數(shù)據(jù)表明：當(dāng)水域截?cái)嚅L(zhǎng)度遠(yuǎn)大于結(jié)構(gòu)尺寸時(shí)，截?cái)噙吔绲倪h(yuǎn)近對(duì)結(jié)構(gòu)入水沖擊的整體速度影響較小。因此，在下面的分析中，水域截?cái)嚅L(zhǎng)度均取1.85 m。

2.2 罰因子討論

在不同罰因子下楔形體入水的計(jì)算結(jié)果見(jiàn)圖3。由此可知：不同的罰因子p對(duì)速度和流固耦合面的合力影響很小。當(dāng)p=0.01時(shí)，楔形體在0.10 s時(shí)刻的速度最大，為0.795 m/s；當(dāng)p=1.00時(shí)，楔形體在0.10 s時(shí)刻的速度最小，為0.794 m/s。罰因子從0.01變化到1.00，0.10 s時(shí)刻的剩余速度只相差0.001 m/s，相對(duì)誤差不到0.13%。因此，罰因子的變化對(duì)剩余速度的影響可忽略不計(jì)。p=0.10時(shí)，楔形體滑移能最小。在純彈性碰撞中，滑移能應(yīng)完全轉(zhuǎn)化為動(dòng)能，因此分析時(shí)應(yīng)盡可能控制滑移能在很小的范圍內(nèi)，越接近0越好，所以建議p取0.10。

2.3 楔形體壓力討論

1 khz以上高頻濾波后楔形體表面測(cè)點(diǎn)壓力的響應(yīng)曲線見(jiàn)圖4。由于測(cè)點(diǎn)選擇在楔形體的壁面上，其壓力存在較大的數(shù)值振蕩。

經(jīng)反復(fù)數(shù)值試驗(yàn)發(fā)現(xiàn)，當(dāng)某時(shí)間步的流固耦合界面剛好超過(guò)流體網(wǎng)格的邊界節(jié)點(diǎn)且超過(guò)量較小時(shí)，該流體節(jié)點(diǎn)會(huì)將流固耦合面附近的流體網(wǎng)格“穿透”，形成向內(nèi)“穿刺”現(xiàn)象，“穿刺”距離未超過(guò)1個(gè)單元格，見(jiàn)圖5。這可能是引起邊界壓力波動(dòng)的根本原因。

為進(jìn)一步探索壓力波動(dòng)的原因，采用數(shù)值試驗(yàn)方法分析測(cè)點(diǎn)1和測(cè)點(diǎn)6附近不同位置處的壓力。距離楔形體測(cè)點(diǎn)6法向距離分別為0、1和2個(gè)網(wǎng)格處的壓力曲線見(jiàn)圖6a）。3個(gè)位置到楔形體壁面法向的距離分別為0、2和4 mm。計(jì)算中心區(qū)域網(wǎng)格大小為2 mm（對(duì)角線長(zhǎng)約為2.8 mm）。由圖6a）可以看出：距測(cè)點(diǎn)6越遠(yuǎn)的點(diǎn)，壓力振蕩越?。慌c測(cè)點(diǎn)6距離為0和1個(gè)網(wǎng)格處的點(diǎn)，因未避免流體向內(nèi)

個(gè)網(wǎng)格處的點(diǎn)壓力振蕩幾乎消失。

距離楔形體測(cè)點(diǎn)1法向距離為3和4 mm處的壓力響應(yīng)曲線見(jiàn)圖6b）。由此可以看出，這2個(gè)點(diǎn)的壓力振蕩很小，壓力響應(yīng)基本一致。這說(shuō)明在法向距離大于1個(gè)網(wǎng)格對(duì)角線距離的位置壓力較平穩(wěn)（壓力曲線自身小的振蕩是由小體積應(yīng)變引起的），可避免ALE網(wǎng)格“穿刺”引起的振蕩。為獲得相對(duì)穩(wěn)定的壓力值，耦合界面壓力記錄位置應(yīng)該選擇在距離流固耦合界面法向1個(gè)網(wǎng)格對(duì)角線處。在此位置上，盡管會(huì)存在一定的壓力損失，但所獲得的壓力是穩(wěn)定的，能相對(duì)真實(shí)地反映實(shí)際壓力。

3 試驗(yàn)結(jié)果對(duì)比分析

為說(shuō)明上述方法的合理性，將模擬結(jié)果與相應(yīng)的試驗(yàn)結(jié)果進(jìn)行對(duì)比。

文獻(xiàn)[16]楔形體速度的試驗(yàn)數(shù)據(jù)與本文數(shù)值模擬結(jié)果對(duì)比見(jiàn)圖7。由此可知：數(shù)值模擬的結(jié)果與ZHAO模型的結(jié)果一致性最好，并介于von Karman模型和Wagner模型結(jié)果之間；與試驗(yàn)結(jié)果相比，數(shù)值模擬結(jié)果中間一段存在一定差異，其原因可能是本文計(jì)算簡(jiǎn)化成平面應(yīng)變問(wèn)題，只考慮水池長(zhǎng)度方向上的邊界影響，未考慮水池寬度方向的邊界影響。

試驗(yàn)與數(shù)值模擬的測(cè)點(diǎn)壓力對(duì)比見(jiàn)圖8。在壓力波形上，數(shù)值模擬與試驗(yàn)較一致。靠近楔形體頂面的測(cè)點(diǎn)數(shù)值模擬壓力曲線存在較小的振蕩。該振蕩可能是由水池長(zhǎng)度方向的邊界反射和數(shù)值計(jì)算中流固耦合界面接近自由表面引起的。

隨著測(cè)點(diǎn)到楔形體頂點(diǎn)距離的增加，試驗(yàn)中測(cè)點(diǎn)壓力上升的起始時(shí)間越來(lái)越早于數(shù)值模擬，在最遠(yuǎn)的測(cè)點(diǎn)12上升時(shí)間差達(dá)到約5 ms。這可能是由試驗(yàn)的三維效應(yīng)和數(shù)值模擬中所記錄的位置與其法向存在3 mm距離共同導(dǎo)致的。隨著試驗(yàn)時(shí)間的增加，因楔形體厚度為1.2 m而水池在該方向上只有2.0 m厚，其壓力反射效果逐漸明顯，進(jìn)而導(dǎo)致峰值前移。

部分測(cè)點(diǎn)的試驗(yàn)與數(shù)值模擬的壓力峰值及其相對(duì)偏差見(jiàn)表1。由此可以看出：測(cè)點(diǎn)1的試驗(yàn)與數(shù)值模擬結(jié)果誤差最大，達(dá)到21.9%；誤差最小點(diǎn)為測(cè)點(diǎn)5，在給定有效數(shù)字范圍內(nèi)，其峰值正好一致。對(duì)沖擊問(wèn)題而言，該偏差在工程分析中完全可以接受。

4 結(jié) 論

基于ALE方法分析研究楔形體結(jié)構(gòu)入水問(wèn)題，并通過(guò)與試驗(yàn)數(shù)據(jù)的對(duì)比驗(yàn)證該方法的合理性。該方法可應(yīng)用于以慣性效應(yīng)為主的結(jié)構(gòu)入水流固耦合問(wèn)題?；谠摲椒ǎㄟ^(guò)系列數(shù)值分析，獲得以下建議。

（1） 采用罰函數(shù)模擬流固耦合時(shí)，因滑移能應(yīng)完全轉(zhuǎn)化為動(dòng)能，建議以最小滑移能來(lái)確定罰因子的取值。本文算例根據(jù)數(shù)值模擬結(jié)果建議罰函數(shù)罰因子取0.10。

（2） 模擬無(wú)限水域時(shí)，其截?cái)噙吔鐟?yīng)遠(yuǎn)離結(jié)構(gòu)，建議距離結(jié)構(gòu)3倍以上結(jié)構(gòu)尺寸為宜。

（3） 由于結(jié)構(gòu)表面水壓力振蕩較大，建議以距離結(jié)構(gòu)表面法向距離為1.4倍單元長(zhǎng)度網(wǎng)格對(duì)角線位置處的壓力為流固耦合面的表面壓力，評(píng)估結(jié)構(gòu)表面的壓力情況。

參考文獻(xiàn)：

[1] 鄒星， 李海濤， 宗智. 基于VOF 模型的結(jié)構(gòu)物出水過(guò)程數(shù)值模擬[J]. 武漢理工大學(xué)學(xué)報(bào)（信息與管理工程版）， 2012， 34（5）： 558561. DOI： 10.3963/j.issn.20953852.2012.5.007.

[2] 馬慶鵬， 魏英杰， 王聰， 等. 錐頭圓柱體高速入水空泡數(shù)值模擬[J]. 北京航空航天大學(xué)學(xué)報(bào)， 2014， 40（2）： 204209.

[3] JIANG C X， SHUAI Z J， ZHANG X Y， et al. Numerical study on transient behavior of waterentry supercavitating flow around a cylindrical projectile influenced by turbulent dragreducing additives[J]. Applied Thermal Engineering， 2016， 104： 450460. DOI： 10.1016/j.applthermaleng.2016.05.102.

[4] FACCI A L ， PORFIRI M ， UBERTINI S. Threedimensional water entry of a solid body： A computational study[J]. Journal of Fluids and Structures， 2016， 66： 3653. DOI： 10.1016/j.jfluidstructs.2016.07.015.

[5] ABRAHAM J P， GORMAN J M， FRANCO R， et al. Modeling and numerical simulation of forces acting on a sphere during earlywater entry[J]. Ocean Engineering， 2014， 76： 19. DOI： 10.1016/j.oceaneng.2013.11.015.

[6] OGER G， DORING M， ALESSANDRINI B， et al. Twodimensional SPH simulations of wedge water entries[J]. Journal of Computational Physics， 2006， 213（2）： 803822. DOI： 10.1016/j.jcp.2005.09.004.

[7] JI Z， XU F， TAKAHASHI A， et al. Large scale water entry simulation with smoothed particle hydrodynamics on single and multiGPU systems[J]. Computer Physics Communications， 2016， 209： 112. DOI： 10.1016/j.cpc.2016.05.016.

[8] GONG K， SHAO S D， LIU H， et al. Twophase SPH simulation of fluidstructure interactions[J]. Journal of Fluids and Structures， 2016， 65： 155179. DOI： 10.1016/j.jfluidstructs.2016.05.012.

[9] SHAO S D. Incompressible SPH simulation of water entry of a freefalling object[J]. International Journal for Numerical Methods in Fluids， 2009， 59（1）： 91115. DOI： 10.1002/fld.1813.

[10] 孫琦， 周軍， 林鵬. 基于LSDYNA的彈體撞水過(guò)程流固耦合動(dòng)力分析[J]. 系統(tǒng)仿真學(xué)報(bào)， 2010， 22（6）： 14981501.

[11] 鄭金偉， 宗智. 三維剛體橢圓頭結(jié)構(gòu)高速傾斜入水沖擊模擬[J]. 船海工程， 2012， 41（3）： 79. DOI： 10.3963/j.issn.16717953.2012.03.003.

[12] 張?jiān)狼啵?徐緋， 金思雅， 等. 飛船返回艙水上回收的沖擊響應(yīng)和入水姿態(tài)分析[J]. 振動(dòng)與沖擊， 2014， 33（18）： 204208. DOI： 110.13465/j.cnki.jvs.2014.18.033.

[13] WANG S， SOARES C G. Numerical study on the water impact of 3D bodies by an explicit finite element method[J]. Ocean Engineering， 2014， 78： 7388. DOI： 10.1016/j.oceaneng.2013.12.008.

[14] 田媛， 劉均， 汪浩， 等. 砰擊載荷下輕質(zhì)波紋夾芯夾層板動(dòng)力響應(yīng)特性分析[J]. 船舶力學(xué)， 2016， 20（10）： 12991307. DOI： 10.3969/j.issn.10077294.2016.10.010.

[15] 李艷臣， 熊偉， 方佩文. 拋落式玻璃鋼救生艇入水分析[J]. 船舶工程， 2016， 38（2）： 1417. DOI： 10.13788/j.cnki.cbgc.2016.02.014.

[16] YETTOU E， DESROCHERS A， CHAMPOUX Y. Experimental study on water impact of a symmetrical wedge[J]. Fluid Dynamics Research， 2006， 38（1）： 4766. DOI： 10.1016/j.fluiddyn.2005.09.003.

（編輯 武曉英）