許蘊(yùn)蕾

（海軍駐上海地區(qū)艦艇設(shè)計(jì)研究軍代室 上海200011）

基于顯式有限元方法的二維楔形體砰擊壓力系數(shù)預(yù)報(bào)

許蘊(yùn)蕾

（海軍駐上海地區(qū)艦艇設(shè)計(jì)研究軍代室 上海200011）

利用LS-DYNA軟件對二維楔形體入水問題進(jìn)行研究。基于顯式有限元方法，選用任意拉格朗日-歐拉算法，建立了包含空氣、水和楔形體的完全耦合的二維有限元模型，研究流場的射流現(xiàn)象與壓力變化情況，預(yù)報(bào)了二維楔形體砰擊壓力系數(shù)，并與已公開發(fā)表的模型試驗(yàn)結(jié)果吻合較好，從而為后續(xù)砰擊載荷計(jì)算提供可靠的方法。

LS-DYNA軟件；二維楔形體；拉格朗日-歐拉算法；砰擊壓力系數(shù)

引 言

船舶在海上航行時(shí)，經(jīng)常會遇到惡劣海況，砰擊問題越來越受到學(xué)者們的關(guān)注。船舶遭受砰擊的瞬時(shí)，底部受到巨大的沖擊力，船體的垂向加速度會突然改變并隨即出現(xiàn)高頻振動，嚴(yán)重時(shí)會造成船體結(jié)構(gòu)的破壞。因此對于船舶設(shè)計(jì)工作者來說，必須準(zhǔn)確預(yù)報(bào)船體的艏部砰擊壓力。

砰擊問題是一個(gè)包含結(jié)構(gòu)、空氣和水三種介質(zhì)耦合的非線性定常問題。自從Von Karman和Wagner開展研究以來，針對砰擊問題，人們開展了大量的理論、數(shù)值和試驗(yàn)研究工作。Zhao等人［1］研究了邊界元方法，預(yù)報(bào)了二維楔形體和艏外飄入水砰擊問題，并開展了試驗(yàn)對比分析工作。駱寒冰［2］采用LS-DYNA軟件，對二維楔形體入水問題進(jìn)行數(shù)值預(yù)報(bào)，并與Zhao等人的模型試驗(yàn)結(jié)果進(jìn)行了比較，結(jié)果吻合較好。陳震［3］給出二維楔形體入水的基本理論，并采用有限元軟件MSC-Dytran對二維楔形體的入水過程進(jìn)行仿真模擬，獲得砰擊面上壓力分布以及入水角度對砰擊載荷的影響規(guī)律。

本文選用任意拉格朗日-歐拉算法，用拉格朗日網(wǎng)格描述結(jié)構(gòu)，多物質(zhì)歐拉網(wǎng)格描述空氣和水。采用顯式有限元方法、罰函數(shù)耦合算法，來模擬歐拉流體與拉格朗日固體接觸面之間的耦合效應(yīng)。二維楔形體模型，利用LS-DYNA軟件，對斜升角為10°的剛性楔形入水問題進(jìn)行仿真模擬，將數(shù)值結(jié)果與文獻(xiàn)［4］中計(jì)算結(jié)果進(jìn)行對比，驗(yàn)證LS-DYNA軟件在流固耦合分析中的可靠性。同時(shí)，本文提出在楔形體入水問題中利用ANSYS前處理實(shí)現(xiàn)位移、速度曲線控制的方法，討論了楔形體入水時(shí)液面變化情況和噴射區(qū)的特點(diǎn)以及流場壓力變化情況，重點(diǎn)討論了斜升角為30°的剛性楔形體入水問題，與Zhao等人的模型實(shí)驗(yàn)結(jié)果和理論計(jì)算結(jié)果進(jìn)行了對比，得到了砰擊壓力系數(shù)。

1 ALE法

ALE法可以處理整個(gè)物體有空間的大位移并且本身有大變形等問題，ALE法最早出現(xiàn)于數(shù)值模擬流體動力學(xué)問題的有限差分方法中。在內(nèi)部網(wǎng)格的劃分上，吸收了歐拉網(wǎng)格的長處，又不完全相同，ALE法不僅使內(nèi)部網(wǎng)格單元獨(dú)立于物質(zhì)實(shí)體而存在，還可以根據(jù)定義的參數(shù)在求解過程中適當(dāng)調(diào)整網(wǎng)格位置，避免出現(xiàn)嚴(yán)重的畸變。此法可使網(wǎng)格與網(wǎng)格之間的物質(zhì)流動。與歐拉網(wǎng)格一樣，ALE網(wǎng)格有兩層網(wǎng)格重疊在一起，空間網(wǎng)格可以在空間任意運(yùn)動，其余與歐拉方法描述一樣，在兩層網(wǎng)格中進(jìn)行物質(zhì)輸送。在結(jié)構(gòu)邊界運(yùn)動的處理上與拉格朗日方法的特點(diǎn)類似，ALE法能夠有效的跟蹤物質(zhì)結(jié)構(gòu)邊界的運(yùn)動。

ALE法執(zhí)行自動重分區(qū)時(shí)，采用拉格朗日和歐拉兩種算法，包含伴隨重映射的拉格朗日時(shí)間步。ALE法可以實(shí)現(xiàn)流體-固體耦合的動態(tài)分析，可以克服單元嚴(yán)重畸變引起的數(shù)值計(jì)算困難。ALE法先執(zhí)行一個(gè)或幾個(gè)拉格朗日時(shí)間步計(jì)算，此時(shí)材料的流動使單元網(wǎng)格產(chǎn)生變形，接著執(zhí)行ALE時(shí)間步計(jì)算：

（1）變形后的物體邊界條件保持不變，充分網(wǎng)格內(nèi)部單元，網(wǎng)格的拓?fù)潢P(guān)系保持不變；

（2）將變形網(wǎng)格中的單元變量和節(jié)點(diǎn)速度矢量輸運(yùn)到重分后的新網(wǎng)格中，稱為Advection Step。

2 二維楔形體入水問題仿真分析

2.1 建立仿真模型

圖1為楔形體的砰擊物理模型，楔形結(jié)構(gòu)為二維對稱結(jié)構(gòu)，只在z軸方向運(yùn)動。將靜水面處取為xoy平面，在y軸方向取單位長度。β為楔形體入水角，Ωfluid代表液體歐拉域，Ωair代表空氣歐拉域，ve為入水速度。

圖1 楔形體砰擊計(jì)算模型

楔形結(jié)構(gòu)用殼單元描述，空氣和水的流域用六面體歐拉單元表述。在歐拉域中，不同位置的網(wǎng)格大小有變化。越靠近邊界處的網(wǎng)格劃分越粗，而楔形體結(jié)構(gòu)位置附近的網(wǎng)格劃分則較細(xì)。選取楔形結(jié)構(gòu)的外側(cè)作為流固耦合作用面，楔形體的外表面法向量調(diào)整為指向流體。為了形成一個(gè)封閉的耦合面，在建模時(shí)選取對稱的楔形結(jié)構(gòu)的1/2部分，并將對稱面處的結(jié)構(gòu)單元定義為虛擬單元。楔形體只有一個(gè)自由度在z軸方向，模型只劃分一個(gè)網(wǎng)格沿y方向，并在該方向的單元表面設(shè)置剛性墻邊界條件模擬二維無限大流場。下頁圖2為楔形體部分網(wǎng)格模型。楔形體下落初速度為6.15 m/s。

圖2 楔形體部分網(wǎng)格模型

2.2 狀態(tài)方程

本文中流體的壓力采用了線性多項(xiàng)式狀態(tài)方程來描述，線性多項(xiàng)式狀態(tài)方程表示單位體積內(nèi)能的線性關(guān)系，壓力值由式（1）給定：

式中：C0、C1、C2、C3、C4、C5和C6為常數(shù)。如果式中μ＜0，則C2μ2和C6μ2兩項(xiàng)設(shè)置為0；其中μ=1/V-1，V表示相對體積。

線性多項(xiàng)式狀態(tài)方程可以用于模擬符合γ律狀態(tài)方程的氣體，各項(xiàng)系數(shù)可設(shè)置為C0=C1=C2=C3= C6=0，C4=C5=γ-1（γ為單位熱值律）。

表1 線性多項(xiàng)式狀態(tài)方程中空氣和水的參數(shù)值

2.3 砰擊壓力系數(shù)

基于沖量砰擊理論，楔形面處的砰擊力峰值可以用式（2）表示：

因此，無因次壓力系數(shù)k表達(dá)式為：

式中： ρ為水的質(zhì)量密度；V為楔形面在法向的速度分量；k為無因次化的壓力峰值系數(shù)。

3 DYNA流固耦合分析要點(diǎn)

3.1 罰函數(shù)耦合

為檢查每一個(gè)節(jié)點(diǎn)對主表面的穿透，用罰函數(shù)耦合系數(shù)追蹤結(jié)構(gòu)從物質(zhì)-拉格朗日節(jié)點(diǎn)和歐拉流體間的相對位移d。罰函數(shù)的物理意義相當(dāng)于為了限制從節(jié)點(diǎn)和主面之間的穿透，在節(jié)點(diǎn)和被穿透面之間放置一個(gè)法向彈簧，如果不出現(xiàn)穿透，就不進(jìn)行任何操作。如果有穿透發(fā)生，即從節(jié)點(diǎn)對主物質(zhì)表面的穿透，則從節(jié)點(diǎn)與被穿透表面之間引入一個(gè)大小與穿透深度成正比的較大界面接觸力。罰函數(shù)數(shù)值的大小不僅影響穿透現(xiàn)象的發(fā)生，如果罰函數(shù)數(shù)值太大，還會使物質(zhì)產(chǎn)生鋼化現(xiàn)象，物質(zhì)表面過分剛硬會影響程序計(jì)算的穩(wěn)定。

3.2 剛性體模型速度位移曲線

剛性體模型在顯示動力學(xué)分析中具有非常重要的意義。當(dāng)定義一個(gè)剛體后，剛體內(nèi)所有節(jié)點(diǎn)的自由度都耦合到剛性體的質(zhì)量中心上，每個(gè)時(shí)間步的節(jié)點(diǎn)力和力矩合成為作用在剛性體上的力和力矩。然后計(jì)算剛性體的運(yùn)動，再轉(zhuǎn)換到節(jié)點(diǎn)位移。由于無論定義多少個(gè)節(jié)點(diǎn)，剛體僅有6個(gè)自由度，因此將有限元模型中剛硬部分用剛性體模型定義，可以大大縮減顯示分析的計(jì)算時(shí)間。

定義剛體模型可約束平移參數(shù)和轉(zhuǎn)動約束參數(shù)，所有剛體的加載必須在前處理ANSYS中定義的part號上而不能使用在組件上。在剛體上施加位移條件，相當(dāng)于在DYNA中定義位移曲線。在前處理ANSYS中定義數(shù)組參數(shù)位移曲線是將時(shí)間間隔和對應(yīng)的位移值集合在一起。它可以分成兩部分，一部分為時(shí)間間隔，另一部分為與時(shí)間間隔對應(yīng)的位移值，定義位移曲線之前，必須首先定義位移曲線的縱坐標(biāo)和橫坐標(biāo)值所對應(yīng)的數(shù)組，且其維數(shù)必須相同。因此，在計(jì)算楔形體以恒定速度入水時(shí)，可根據(jù)此方法控制楔形體速度保持恒定。

3.3 壓力傳感器

在DYNA中對剛性體模型的單元進(jìn)行壓力的提取，需要定義壓力傳感器，在輸入文件中定義壓力傳感器位置，這些位置同樣是拉格朗日算符單元所在的位置，壓力傳感器壓力的提取是從流體單元中得到的壓力。

3.4 時(shí)間步長尺度

在LS-DYNA中自動選擇的時(shí)間步長是基于以下條件的：在一個(gè)單位時(shí)間步內(nèi)聲波不允許穿過最小的單元，時(shí)間步長也受到全系統(tǒng)剛度的限制，然后與時(shí)間步長相關(guān)的過大的接觸剛度可能會產(chǎn)生非物理的能量的接觸，時(shí)間步長影響了震蕩壓力信號，其作用效果與較高的接觸剛度的作用效果一樣，時(shí)間步長的改變對結(jié)果沒有明顯的影響。

4 計(jì)算結(jié)果

4.1 10°斜升角楔形體計(jì)算結(jié)果比較

文獻(xiàn)［4］中二維剛性楔形體斜升角為10°，該模型采用剛性材料以5.425 m/s的速度等速入水。本文利用LS-DYNA軟件進(jìn)行仿真，利用前處理ANSYS中施加位移曲線的方式實(shí)現(xiàn)等速入水，得到如表2所示，砰擊壓力的理論值與數(shù)值解和對應(yīng)位置本文的仿真值之間的比較結(jié)果。本文計(jì)算結(jié)果誤差較小，計(jì)算結(jié)果較為可靠。圖3為砰擊壓力時(shí)歷曲線。

表2 文獻(xiàn)值與本文仿真值之間的比較

圖3 砰擊壓力時(shí)域結(jié)果

4.2 30°斜升角楔形體仿真模型計(jì)算結(jié)果

4.2.1 壓力分布

圖4 -圖6顯示了入水砰擊過程中不同時(shí)刻的水面飛濺現(xiàn)象。圖7 -圖9顯示流場中壓力分布情況。在該計(jì)算模型中并沒有觀察到液體的滲漏現(xiàn)象，表明耦合模擬合理?？梢钥闯?，在射流處，砰擊壓力急劇減少，可以忽略不計(jì)；在外域水中，離液面升高處越遠(yuǎn)，壓力越小。液面分離之前，砰擊壓力的最大值出現(xiàn)在液面升高處。在液面分離之前的入水過程，砰擊壓力的最大值明顯降低，砰擊壓力的最大值出現(xiàn)的位置，將從楔形體底部向上端折角處移動。這是由于仿真楔形體質(zhì)量較輕，速度下降較快，導(dǎo)致砰擊壓力逐漸減少，也與壓力傳感器得到的壓力相對應(yīng)，位置越靠上的壓力傳感器測得的壓力最大值越小。

圖4 入水初期

圖5 入水中期

圖6 入水后期

圖7 0.051 s 時(shí)流場壓力分布

圖8 0.117 s 時(shí)流場壓力分布

4.2.2 時(shí)域內(nèi)砰擊壓力曲線

圖10顯示在對應(yīng)于Zhao等［1］模型試驗(yàn)中，壓力傳感器位置處的仿真模型壓力傳感器的壓力時(shí)歷曲線。最大砰擊壓力的逐漸降低，同時(shí)也對應(yīng)于圖7 -圖9顯示的砰擊壓力最大值的降低。

圖10 楔形體入水實(shí)驗(yàn)?zāi)Ｐ?/p>

圖11所示為5個(gè)壓力傳感器測得的砰擊壓力時(shí)域預(yù)報(bào)結(jié)果。橫坐標(biāo)是入水時(shí)間，5個(gè)壓力傳感器位置對應(yīng)于圖10所示試驗(yàn)壓力測點(diǎn)位置。

圖11 不同位置砰擊壓力時(shí)歷曲線

4.2.3 砰擊壓力系數(shù)

文獻(xiàn)［1］中給出不同斜升角楔形體的不同理論算法算得的砰擊壓力系數(shù)，對于30°斜升角的砰擊壓力系數(shù)曲線如圖12所示。

圖12 砰擊壓力空間分布比較

本文計(jì)算結(jié)果與文獻(xiàn)［1］ 中模型試驗(yàn)結(jié)果吻合較好，圖中所示4個(gè)點(diǎn)的數(shù)據(jù)來源于前4個(gè)壓力傳感器。其中，文獻(xiàn)［1］ 給出的理論計(jì)算結(jié)果是樣條曲線的形式，是由于理論計(jì)算中選取計(jì)算點(diǎn)數(shù)目較多，將一系列計(jì)算點(diǎn)結(jié)果進(jìn)行曲線擬合，便得到圖12中的計(jì)算結(jié)果。

本文利用LS-DYNA軟件進(jìn)行仿真模擬，考慮到計(jì)算時(shí)間問題，選取4個(gè)點(diǎn)的結(jié)果進(jìn)行比較，因此本文計(jì)算給出的結(jié)果是獨(dú)立點(diǎn)，這些獨(dú)立點(diǎn)的計(jì)算結(jié)果與理論值吻合較好。其中文獻(xiàn)［1］ 給出的結(jié)果中，對縱坐標(biāo)做了無因次化處理，Z是楔形體液面距離楔形體底部的垂向高度，V表示楔形體在入水不同時(shí)刻的垂向速度。在砰擊過程中，LS-DYNA預(yù)報(bào)的砰擊壓力系數(shù)與試驗(yàn)值吻合非常好，反映顯式有限元技術(shù)在預(yù)報(bào)二維剛性體砰擊問題上的可靠性以及本文所采用仿真方法的合理性。

5 結(jié) 論

本文采用LS-DYNA軟件，基于顯式有限元方法，采用任意拉格朗日-歐拉算法和罰函數(shù)耦合算法，針對文獻(xiàn)［1］ 中楔形體入水砰擊試驗(yàn)結(jié)果以及理論計(jì)算結(jié)果，建立二維剛性楔形體模型，展開砰擊壓力系數(shù)預(yù)報(bào)工作。預(yù)報(bào)的10°楔形體砰擊壓力值與文獻(xiàn)［4］中吻合較好，驗(yàn)證了LS-DYNA在二維楔形體入水仿真中比較可靠。

文中重點(diǎn)分析與Zhao等［1］楔形體入水砰擊實(shí)驗(yàn)相同30°斜升角的楔形體入水模型，得到不同測點(diǎn)的砰擊壓力時(shí)域分布以及求得砰擊壓力系數(shù)，并與其理論推導(dǎo)結(jié)果進(jìn)行比較，與模型試驗(yàn)以及理論推導(dǎo)結(jié)果吻合較好。預(yù)報(bào)的砰擊壓力結(jié)果可正確描述入水砰擊的整個(gè)過程（包括入水初期、入水中期的砰擊過程階段）。LS-DYNA有效預(yù)報(bào)了砰擊壓力在楔形體表面的時(shí)間、空間分布，合理反映出砰擊過程的射流以及內(nèi)域和外域中的壓力分布，為后續(xù)計(jì)算砰擊壓力提供了砰擊壓力系數(shù)求解的方法。

［1］ ZHAO R，F(xiàn)ALTINSEN O M.Water entry of twodimensional bodies［J］. Journal of Fluid Mechanics，1993，246：593-612.

［2］ 駱寒冰，吳景健，王珊，等.基于顯式有限元方法的二維楔形剛體入水砰擊載荷并行計(jì)算預(yù)報(bào)［J］.船舶力學(xué)，2012，16（8）：907-904.

［3］ 陳震，肖熙.二維楔形體入水砰擊仿真研究［J］.上海交通大學(xué)學(xué)報(bào)，2007，41（9）：1425-1428.

［4］ AQUELET N，SOULI M，OLOVSSON L. Euler-Lagrange coupling with damping effects：Application to slamming problems［J］.Computer methods in applied mechanics and engineering，2006，195（1-3）：110-132.

Slamming pressure coef fi cient prediction for water entry of 2D wedge body with explicit fi nite element method

XU Yun-lei
(Representative Of fi ce of Naval Warship Design & Research, Shanghai 200011, China)

The water entry of a two-dimension wedge body has been studied by LS-DYNA software. Based on the explicit finite element and Arbitrary Lagrangian-Eulerian (ALE) algorithm, it builds up a two-dimension coupling fi nite element model including air, water and wedge, studies jet fl ow and pressure change, and predicts the slamming pressure coef fi cient of the two-dimension wedge body. The results are consistent with the published model test results, which can provide the reliable method for calculating of the slamming pressure coef fi cient.

LS-DYNA software; 2D wedge body; Arbitrary Lagrangian-Eulerian (ALE) algorithm; slamming pressure coef fi cient

U661.44

A

1001-9855（2014）03-0042-06

2014-02-18 ；

2014-04-16

許蘊(yùn)蕾（1981-），女，碩士，工程師，研究方向：船舶總體設(shè)計(jì)。