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

許蘊(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ù)

引 言

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

砰擊問題是一個包含結(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軟件在流固耦合分析中的可靠性。同時，本文提出在楔形體入水問題中利用ANSYS前處理實(shí)現(xiàn)位移、速度曲線控制的方法，討論了楔形體入水時液面變化情況和噴射區(qū)的特點(diǎn)以及流場壓力變化情況，重點(diǎn)討論了斜升角為30°的剛性楔形體入水問題，與Zhao等人的模型實(shí)驗(yàn)結(jié)果和理論計(jì)算結(jié)果進(jìn)行了對比，得到了砰擊壓力系數(shù)。

1 ALE法

ALE法可以處理整個物體有空間的大位移并且本身有大變形等問題，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ū)時，采用拉格朗日和歐拉兩種算法，包含伴隨重映射的拉格朗日時間步。ALE法可以實(shí)現(xiàn)流體-固體耦合的動態(tài)分析，可以克服單元嚴(yán)重畸變引起的數(shù)值計(jì)算困難。ALE法先執(zhí)行一個或幾個拉格朗日時間步計(jì)算，此時材料的流動使單元網(wǎng)格產(chǎn)生變形，接著執(zhí)行ALE時間步計(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)整為指向流體。為了形成一個封閉的耦合面，在建模時選取對稱的楔形結(jié)構(gòu)的1/2部分，并將對稱面處的結(jié)構(gòu)單元定義為虛擬單元。楔形體只有一個自由度在z軸方向，模型只劃分一個網(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ù)耦合

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

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

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

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

3.3 壓力傳感器

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

3.4 時間步長尺度

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

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

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

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

4.2.1 壓力分布

圖4 -圖6顯示了入水砰擊過程中不同時刻的水面飛濺現(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 時流場壓力分布

圖8 0.117 s 時流場壓力分布

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

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

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

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

圖11 不同位置砰擊壓力時歷曲線

4.2.3 砰擊壓力系數(shù)

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

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

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

本文利用LS-DYNA軟件進(jìn)行仿真模擬，考慮到計(jì)算時間問題，選取4個點(diǎn)的結(jié)果進(jìn)行比較，因此本文計(jì)算給出的結(jié)果是獨(dú)立點(diǎn)，這些獨(dú)立點(diǎn)的計(jì)算結(jié)果與理論值吻合較好。其中文獻(xiàn)［1］ 給出的結(jié)果中，對縱坐標(biāo)做了無因次化處理，Z是楔形體液面距離楔形體底部的垂向高度，V表示楔形體在入水不同時刻的垂向速度。在砰擊過程中，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ù)，并與其理論推導(dǎo)結(jié)果進(jìn)行比較，與模型試驗(yàn)以及理論推導(dǎo)結(jié)果吻合較好。預(yù)報(bào)的砰擊壓力結(jié)果可正確描述入水砰擊的整個過程（包括入水初期、入水中期的砰擊過程階段）。LS-DYNA有效預(yù)報(bào)了砰擊壓力在楔形體表面的時間、空間分布，合理反映出砰擊過程的射流以及內(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ì)。