王 玲，朱仁慶，李紅艷

(江蘇科技大學(xué) 船舶與海洋工程學(xué)院，江蘇 鎮(zhèn)江 212003)

入水砰擊問題涉及流體與固體相互耦合作用，是水彈性力學(xué)研究領(lǐng)域中的一個(gè)熱點(diǎn)。當(dāng)沖擊發(fā)生時(shí)，水彈性力學(xué)的相關(guān)特性主要取決于流體物理特性、結(jié)構(gòu)物幾何特性和物理特性以及入水的初速度及入水角度等因素。

由于實(shí)尺度試驗(yàn)成本太高、操作困難等原因，模型試驗(yàn)更受研究者的青睞[1–8]。張?jiān)狼嗟萚9]進(jìn)行了楔形體垂直入水沖擊實(shí)驗(yàn)，通過濾波消除高頻信號，得到了合理的加速度響應(yīng)值。在數(shù)值方法方面，早期入水問題的研究對象多為剛體，但是隨著理論知識的豐富，研究彈性體入水也逐漸增多[10–12]。張智[13]用 Fluent軟件模擬了楔形體垂直入水砰擊和轉(zhuǎn)角度入水砰擊的2種情況，并得到了加速度、速度和砰擊壓力等隨時(shí)間變化的曲線，同時(shí)，他還對砰擊影響進(jìn)行了敏感性分析。

本文數(shù)值模擬部分主要采用Ansys軟件中Fluent和Workbentch模塊，外加自主編程，結(jié)合動網(wǎng)格技術(shù)來模擬和計(jì)算剛性和彈性楔形體的入水過程。

1 入水問題的數(shù)學(xué)模型和數(shù)值模型

1.1 結(jié)構(gòu)入水水彈性力學(xué)的微分控制方程

考慮一個(gè)彈性體，其區(qū)域?yàn)閂s，入水的流體區(qū)域?yàn)閂L，彈性體表面為Sp，彈性體和流體的交界面為S（見圖1）。

設(shè)彈性體是各向同性且為線變形，彈性體的運(yùn)動微分方程為：

式中： ρs為彈性體的密度；下標(biāo)S為彈性體； ui為彈性體質(zhì)點(diǎn)沿 xi方向的位移； Fi為作用于彈性體上的體積力； σij為應(yīng)力，沿xi方向的應(yīng)力；t為時(shí)間。

應(yīng)變-位移關(guān)系為：

式中：eij為6個(gè)應(yīng)變分量，即e11，e22，e33，e12=e21，e23=e32，e31=e13。

各向同性彈性體的應(yīng)力-應(yīng)變關(guān)系為：

式中： λ和μ為拉梅常數(shù)。它們和楊氏模量E，泊松比ν的關(guān)系為

式中， ekk=e11+e22+e33體積膨脹應(yīng)變；k為啞標(biāo)。

對于流體部分的區(qū)域 VL；其運(yùn)動方程為N-S方程

式中： ρL為流體密度為流體質(zhì)點(diǎn)的加速度；p為壓強(qiáng)； Fi為體積力； pij為本構(gòu)關(guān)系，即

連續(xù)方程為

對于不可壓縮的液體，其狀態(tài)方程可取

若考慮壓縮性，對于純液體，其壓強(qiáng)與密度間的關(guān)系可由下列半經(jīng)驗(yàn)的公式給出：

式中：A，B為與液體性質(zhì)有關(guān)的常數(shù)，對于A=7，B=3 000大氣壓，下標(biāo)0表示未經(jīng)擾動時(shí)的量。

若考慮水的微量可壓縮性，則有

式中：K為體積彈性系數(shù)；對于水在 15 ℃ 時(shí)，有K=2.045×109Pa。

在彈性體和流體的交界面S上應(yīng)滿足壓強(qiáng)相等和法向速度相等的條件：

式中： ni為S面上單位法矢沿xi方向的分量，指向彈性體的外部為正。

1.2 入水問題的數(shù)值模型

1）控制方程的離散方法

本文流體域采用有限體積法，結(jié)構(gòu)域采用有限元法。

關(guān)于有限體積法對流域控制方程進(jìn)行離散的過程如下：

首先將連續(xù)方程和運(yùn)動方程統(tǒng)寫為

其中，左側(cè)第1項(xiàng)是瞬態(tài)項(xiàng)，第2項(xiàng)是對流項(xiàng)，右側(cè)依次代表擴(kuò)散項(xiàng)和源項(xiàng)。假定代表速度勢函數(shù)，且守恒方程為定常狀態(tài)，那么對任意控制單元，可以得到積分方程：

通過高斯公式：

變換得到：

最后由耦合界面條件得到：

式中： Nfaces為控制體周圍單元面的數(shù)量； ?f為速度勢函數(shù)在面 f 上的對流； ρfuf·Af為流經(jīng)面 f的質(zhì)量通量； (??)n為 ??在面 f 上法線方向的大??； V為控制單元的體積。

2）流場的求解方法

本文選用了SIMPLEC算法。該方法考慮了SIMPLE算法中省略掉的速度修正方程中的項(xiàng)，改進(jìn)了速度修正方程中的“協(xié)調(diào)性”問題，因此得到了較為準(zhǔn)確的壓力修正值。3）流場與結(jié)構(gòu)耦合方程

式中：Ff和Fs分別是與Gf和Gs相應(yīng)的離散方程。注意到，耦合的流體和結(jié)構(gòu)方程可以分別表示為

2 楔形體入水砰擊試驗(yàn)

2.1 試驗(yàn)?zāi)康?/h3>

自主設(shè)計(jì)模型結(jié)構(gòu)，自主設(shè)計(jì)試驗(yàn)裝置。通過布置在楔形體結(jié)構(gòu)內(nèi)的傳感器，全程實(shí)時(shí)監(jiān)測楔形體結(jié)構(gòu)的下落的壓力變化，為研究楔形體結(jié)構(gòu)物入水砰擊問題的模型試驗(yàn)方法提供補(bǔ)充和改進(jìn)，也為該問題的數(shù)值模擬方法進(jìn)行驗(yàn)證。

2.2 試驗(yàn)?zāi)Ｐ?/h3>

為研究結(jié)構(gòu)入水砰擊過程中的砰擊壓力，結(jié)構(gòu)模型為鋁制模型。模型如圖2所示。

圖3是模型成品圖，在圖3中，為保證傳感器的引線能夠順利牽出，在楔形體頂部開一個(gè)孔。左右2個(gè)開孔，是為了保證楔形體重心仍位于結(jié)構(gòu)中心處。在圖5中，考慮到楔形體水密性，并為了保護(hù)傳感器，在楔形體左右兩端用有機(jī)玻璃進(jìn)行密封。模型的主要參數(shù)如下：

圖2 模型正視圖Fig.2 Model front view

圖3 模型成品圖Fig.3 Figure of finished mode

圖4 測點(diǎn)布置圖Fig.4 Measuring point layout

楔形體主尺度為:1 m×1 m×0.5 m，肋骨 2 道，縱骨3道（其中1道主龍骨）。

楔形體底升角：45°。

楔形體厚度：4 mm。

2.3 測點(diǎn)布置

該試驗(yàn)中的模型是軸對稱結(jié)構(gòu)，因此只布置一邊測點(diǎn)即可。如圖4所示，其中圖4（a）為測點(diǎn)示意圖，圖4（b）表示測點(diǎn)在模型中的真實(shí)位置。

圖5 模型入水裝置示意圖Fig.5 Schematic diagram of the model into water

該實(shí)驗(yàn)中總共布置了5個(gè)測試點(diǎn)，全部位于楔形體中垂面上，P1距頂部10 cm，隨后每2個(gè)點(diǎn)之間間距10 cm，考慮到中間加強(qiáng)筋結(jié)構(gòu)，所以P2和P3之間間隔 20 cm。

2.4 試驗(yàn)儀器

試驗(yàn)儀器包括高靈敏壓力傳感器、數(shù)據(jù)信號采集器、計(jì)算機(jī)、相關(guān)儀器電源和導(dǎo)線等。

壓力傳感器：若干個(gè)，量程為–10 kPa～2 MPa。該傳感器要滿足可以遇水的需求。

信號采集器：1臺，用來采集分析系統(tǒng)測量到的結(jié)構(gòu)入水瞬間的砰擊壓力值。

計(jì)算機(jī)：1臺，用來與信號采集器連接，接收與存儲信號采集器的數(shù)據(jù)。

該試驗(yàn)的采集工作由以上儀器配合完成，壓力傳感器測量得到結(jié)構(gòu)入水的砰擊壓力值，后經(jīng)數(shù)據(jù)信號采集器接收與儲存，最后通過計(jì)算機(jī)讀取。

實(shí)驗(yàn)開始前，先將壓力傳感器進(jìn)行標(biāo)定。

2.5 試驗(yàn)裝置和試驗(yàn)方案

試驗(yàn)水池的尺寸為 5 m×2.5 m×3 m。砰擊試驗(yàn)塔架的高度為5m，塔架的橫梁中端固定一小型起重機(jī)。試驗(yàn)中采用電動葫蘆和手動脫鉤相結(jié)合的方式，并在橫梁上潤滑油，目的是為了盡量減少由起吊裝置造成的速度損失。為了減少砰擊的離散性，將人為因素降到最低，減小偶然誤差對測量結(jié)果的影響，每組試驗(yàn)都進(jìn)行5次重復(fù)試驗(yàn)，2次試驗(yàn)之間大約間隔25 min。與此同時(shí)，等待水面重新平靜，模型充分干燥以及準(zhǔn)備下一次的入水試驗(yàn)。圖5為模型入水裝置示意圖。

試驗(yàn)開始時(shí)，用電動葫蘆將模型拉升至指定高度，然后剪斷吊繩使其自由下落。

2.6 實(shí)驗(yàn)過程及結(jié)果分析

本次試驗(yàn)主要是模擬結(jié)構(gòu)物自由落體的入水砰擊試驗(yàn)。實(shí)驗(yàn)開始時(shí)，首先將模型吊于空中懸臂梁上，然后讓結(jié)構(gòu)在距離自由液面0.5 m，0.8 m和1 m的高度自由下落。

表1為模型入水試驗(yàn)工況表。每次試驗(yàn)開始時(shí)，先將各個(gè)測試儀器歸零，移動起重機(jī)，將模型提升到所在工況的高度，然后用剪刀剪斷吊繩，讓模型自由落下。模型入水過程中，底部首先與水面發(fā)生碰撞，濺射出來的水花沿著中部向兩邊擴(kuò)散，此時(shí)，數(shù)據(jù)采集儀器會自動將壓力傳感器的電壓信號歷程記錄下來。每一次試驗(yàn)結(jié)束后，將模型重新提升至下一個(gè)目標(biāo)高度，等待水面平靜后進(jìn)行下一次試驗(yàn)。

表1 模型入水試驗(yàn)工況表Tab.1 Model test conditions into the water

由圖6的模型落水過程圖可以發(fā)現(xiàn)，模型落入水中后，在其兩側(cè)激起水花。隨著下落高度的增加，水花成噴射狀。

表2表示不同高度下有壓力傳感器測得的在模型中心處各點(diǎn)的砰擊壓力峰值。圖7和圖8分別表示0.5 m和1 m高度各點(diǎn)落水砰擊壓力歷程圖。由圖表可知，隨著下落高度的增加，即砰擊入水速度的增大，砰擊壓力峰值逐漸增大。在楔形體下落過程中，越靠近楔形體底部的點(diǎn)越早接觸水面，得到的砰擊壓力值也越大。由于測點(diǎn)布置時(shí)，未充分考慮浮力與重力的關(guān)系，模型落水瞬間，水流未來得及沖擊到P5點(diǎn)時(shí)，浮力已經(jīng)大于重力，因而P5點(diǎn)未有監(jiān)測值。

綜上所述，隨著下落高度的增加，入水砰擊速度的增大，砰擊壓力峰值逐漸增大；砰擊發(fā)生時(shí)，即底部剛接觸水面的時(shí)候，模型的砰擊壓力迅速增大到峰值，因而越靠近底部受到的砰擊壓力也越大，隨著入水深度地增加，砰擊壓力逐漸減小，直至水面平靜。

3 三維柔性楔形體入水

3.1 幾何建模及網(wǎng)格劃分

本節(jié)研究的三維楔形體結(jié)構(gòu)如圖9所示，圖9（a）表示計(jì)算域，圖9（b）表示三維對稱楔形體平面圖。其中對稱楔形的底邊長為 1 m，高度 1 m，α=β=45°。計(jì)算域尺寸為：高5 m，長15 m，寬度4.4 m，基本滿足不反射邊界條件。此時(shí)楔形體結(jié)構(gòu)完全按照試驗(yàn)?zāi)Ｐ蛠砟M，即與上述剛性楔形體模型不同的是，楔形體主體有2道肋骨，3道縱骨。

圖6 模型落水過程Fig.6 Model drowning process

表2 模型中心處各點(diǎn)的砰擊壓力峰值試驗(yàn)值Tab.2 The test slamming peak of each pointat the center of the model

圖7 0.5 m 高度各點(diǎn)落水砰擊壓力時(shí)間歷程圖Fig.7 The pressure time history diagram of each pointfrom 0.5 m height

圖8 1 m 高度各點(diǎn)落水砰擊壓力時(shí)間歷程圖Fig.8 The pressure time history diagram of each pointfrom 1 m height

圖9 計(jì)算域和結(jié)構(gòu)平面圖Fig.9 The sketch diagram of computational domain and structure plan

與模擬三維剛性楔形體入水近似，對液面附近網(wǎng)格進(jìn)行加密。網(wǎng)格單元總數(shù)為102 453。采用六面體網(wǎng)格能夠減少網(wǎng)格數(shù)量，并且不需要進(jìn)行網(wǎng)格重構(gòu)，極大地減少了計(jì)算量。楔形體結(jié)構(gòu)全部采用四面體網(wǎng)格劃分。網(wǎng)格劃分如圖10所示。

圖10 網(wǎng)格劃分Fig.10 Meshing

3.2 計(jì)算求解

本節(jié)選用3D壓力基求解器，流動為非定常流，考慮重力加速度大小為9.81 m/s2。物理求解模型選擇VOF多相流模型，同時(shí)打開Implicit Body Force，以提高解的穩(wěn)定性。選擇標(biāo)準(zhǔn)壁面函數(shù)（Standard Wall Functions）。邊界條件作如下設(shè)置，將整個(gè)區(qū)域的最上面設(shè)置成壓力出口，其他邊界均設(shè)置為無滑移壁面邊界。在這里采用彈簧光滑模型和網(wǎng)格重構(gòu)模型，同時(shí)開啟6DOF模型。在運(yùn)動屬性中，將楔形體定義為剛體（rigid body），楔形體的質(zhì)量及轉(zhuǎn)動慣量在udf中設(shè)置。

3.3 彈性體入水的響應(yīng)分析

本節(jié)研究三維彈性對稱楔形體距自由液面1 m高度下落的情況。

圖11是非對稱楔形體從1 m高度自由落下不同時(shí)刻楔形體表面壓力分布云圖。從壓力分布云圖可以看出，在入水過程中，楔形體的高壓區(qū)主要分布在楔形體頭部，隨著楔形體下落深度的增加，高壓波逐漸傳播出去，楔形體表面的壓力趨于平穩(wěn)。

圖11 楔形體不同時(shí)刻的壓力分布云圖（h=1 m）Fig.11 Contours of surface pressure of wedge in different time (h=1 m)

從圖12楔形體不同時(shí)刻的等效應(yīng)力云圖中可以看出，應(yīng)力最大區(qū)域依然是楔形體頭部區(qū)域，在t=0.1 s時(shí)刻，達(dá)到了15 520 Pa。隨著楔形體下落深度的增加，應(yīng)力在楔形體表面逐漸趨于穩(wěn)定。

圖13是不同時(shí)刻流體的速度云圖。從圖中可以看出，在楔形體下落過程中，流體的速度增大，在t=0.5 s時(shí)速度達(dá)到了10.06 m/s，最大速度在楔形體底端兩側(cè)頂角處。隨著楔形體下落深度的增加，流體速度逐漸減小。在t=0.2 s時(shí)，速度又增大，此時(shí)是因?yàn)樾ㄐ误w的下落造成了流體飛濺，飛濺的速度較楔形體下落時(shí)流體的速度大，最大速度處仍在楔形體底端兩側(cè)頂角處。

圖12 楔形體不同時(shí)刻的等效應(yīng)力云圖（h=1 m）Fig.12 Contours of von mises stress of wedge in different time(h=1 m)

圖13 流體不同時(shí)刻的速度云圖（h=1 m）Fig.13 Contours of velocity of fluid in different time (h=1 m)

為更全面驗(yàn)證上述試驗(yàn)的可行性與有效性，本文在柔性楔形體入水的計(jì)算中加入了5個(gè)測點(diǎn)，這5個(gè)測點(diǎn)位置與試驗(yàn)?zāi)Ｐ偷臏y點(diǎn)位置相匹配，如圖7（a）所示。

圖14是測點(diǎn)的砰擊壓力圖。從圖中可以看出，最大的砰擊壓力為P1點(diǎn)，依次減小，P5因位置原因，未得到水流的沖擊，因而沒有數(shù)據(jù)，與試驗(yàn)中測點(diǎn)的趨勢相符。

圖14 測點(diǎn)的砰擊壓力圖Fig.14 Slamming pressure peak of test points

表3是數(shù)值模擬測點(diǎn)試驗(yàn)值和模型試驗(yàn)測點(diǎn)砰擊壓力峰值對比。以上4組對比數(shù)據(jù)中，誤差范圍在10%以內(nèi)，排除空氣阻力等，本次數(shù)值模擬結(jié)果基本符合實(shí)際。

從表3中4組對比數(shù)據(jù)中，可以看出考慮彈性體的數(shù)值模擬的砰擊壓力值均大于模型實(shí)驗(yàn)得到的砰擊壓力值，這是由于在考慮彈性體的數(shù)值模擬中，默認(rèn)選擇了先進(jìn)行耦合計(jì)算，然后考慮固體的彈性形變，但是在實(shí)際試驗(yàn)中，耦合與形變同時(shí)發(fā)生。固體變形，相當(dāng)于一個(gè)卸力的效果，使壓力有效降低，考慮彈性體計(jì)算短暫的忽略了這個(gè)變形，所以數(shù)值稍微高一點(diǎn)。

表3 數(shù)值模擬測點(diǎn)試驗(yàn)值和模型試驗(yàn)測點(diǎn)砰擊壓力峰值對比Tab.3 Contrast of slamming pressure peak between numerical simulation and model test points

4 結(jié) 語

本文自行設(shè)計(jì)了楔形體入水試驗(yàn)，以便驗(yàn)證改數(shù)值模擬方法的正確性。自行設(shè)計(jì)了試驗(yàn)?zāi)Ｐ秃驮囼?yàn)裝置，完成對試驗(yàn)?zāi)Ｐ偷募庸ぶ圃欤{(diào)試完整個(gè)水池試驗(yàn)系統(tǒng)后進(jìn)行了不同高度的入水砰擊試驗(yàn)。試驗(yàn)完整記錄了模型入水砰擊的過程，得到了模型在不同高度下入水砰擊過程中的砰擊壓力。最后，進(jìn)行三維彈性體入水的數(shù)值模擬，數(shù)值模擬模型完全參照試驗(yàn)?zāi)Ｐ?，得到了砰擊壓力云圖、等效應(yīng)力云圖以及流體的速度云圖，對比分析了彈性體測點(diǎn)與試驗(yàn)中測點(diǎn)的砰擊壓力，驗(yàn)證了該數(shù)值模擬方法的可行性與有效性。

[1]OCHI M K, MOTTER L E.Prediction of slamming characteristics and hull responses for ship design[J].Transactions SNAME, 1973, Vol.81, 144–177.

[2]PUKHNACHOV V V.Linear approximation in the problem on a blunt body entry in water [J].Din.Sploshnoi Sredy.1979, 38:143–150.

[3]TAKAGI K, 1994, Influence of e1asticity on hydrodynamic impact problem[J], J.Kansai Soc.N.A., Japan, Sep.1994, No.222, pp.97–106.

[4]A.EL MALKI ALAOUI, A.NêME, A.TASSIN, et al.Experimental study of coefficients during vertical water entry of axisymmetric rigid shapes at constant speeds[J].Applied Ocean Research.2012, 37: 183–197.

[5]SUN S Y, SUN S L, WU G X.Oblique water entry of a wedge into waves with gravity effect[J].Jounal of Fluids and structures, 2015, 52: 49–64 (IF, 2.021).

[6]張偉, 郭子濤, 肖新科, 等.彈體高速入水特性實(shí)驗(yàn)研究[J].爆炸與沖擊.2011, 31(6): 579–584.

[7]何春濤, 王聰, 魏英杰, 等.圓柱體垂直入水空泡形態(tài)試驗(yàn)[J].北京航空航天大學(xué)學(xué)報(bào), 2012(11): 1542–1546.

[8]楊衡, 張阿漫, 龔小超, 等.不同頭型彈體低速入水空泡試驗(yàn)研究[J].哈爾濱工程大學(xué)學(xué)報(bào), 2014, 09: 1060–1066.

[9]張?jiān)狼? 徐緋, 金思雅, 等.楔形體入水沖擊響應(yīng)的試驗(yàn)研究及應(yīng)用[J].機(jī)械強(qiáng)度, 2015, 37(2): 226–231.

[10]楊衡, 孫龍泉, 龔小超, 等.彈性結(jié)構(gòu)入水砰擊載荷特性三維數(shù)值模擬研究[J].振動與沖擊, 2014, 33 (19): 28–33.

[11]張偉偉, 金先龍, 劉濤.正交異性復(fù)合材料儲液容器入水過程數(shù)值分析[J].計(jì)算力學(xué)學(xué)報(bào), 2015, 32(15): 633–638, 698.

[12]SUN P, MING F, ZHANG A.Numerical simulation of interactions between free surface and rigid body using a robust SPH method[J].Ocean Engineering, 2015, 98: 32–49.

[13]張智.基于CFD技術(shù)的三維楔形體入水砰擊試驗(yàn)和數(shù)值預(yù)報(bào)研究[D].天津: 天津大學(xué).2013.