王加夏,周天九,劉 昆,陳 月

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

船舶在航行過(guò)程中，會(huì)發(fā)生搖蕩運(yùn)動(dòng)，這時(shí)劇烈的波浪會(huì)對(duì)船體發(fā)生猛烈沖擊，這一現(xiàn)象被稱(chēng)為砰擊[1-2]．波浪砰擊過(guò)程涉及波浪強(qiáng)非線(xiàn)性、瞬時(shí)效應(yīng)的物理特性，也涉及水彈性效應(yīng)及流固耦合等復(fù)雜的分析技術(shù)，一直是海洋工程領(lǐng)域研究的熱點(diǎn)課題之一．嚴(yán)重的砰擊不僅可能會(huì)對(duì)沖擊區(qū)域產(chǎn)生巨大的沖擊壓力，造成局部結(jié)構(gòu)破壞;而且，頻繁的波浪砰擊也會(huì)誘導(dǎo)船體出現(xiàn)二節(jié)點(diǎn)的顫振現(xiàn)象，造成船體垂向彎矩極限值的增加，影響結(jié)構(gòu)的極限強(qiáng)度．

文獻(xiàn)[3]中首先對(duì)砰擊問(wèn)題進(jìn)行了理論研究，利用動(dòng)量定理得出了經(jīng)典的砰擊壓力的計(jì)算公式．在此基礎(chǔ)上，計(jì)入自由液面升高現(xiàn)象，將結(jié)構(gòu)入水砰擊過(guò)程簡(jiǎn)化為平板入水的模型，得到了沖擊理論模型[4]．后續(xù)學(xué)者在不同方面進(jìn)行了大量的修正．其中較為突出并且應(yīng)用較廣的分別為廣義Wagner模型[5]和修正的Logvinovich模型[6]．文獻(xiàn)[7]中對(duì)一艘豪華游輪進(jìn)行了尾部砰擊荷載和振動(dòng)響應(yīng)的模型試驗(yàn)研究，提出了一種用試驗(yàn)測(cè)量尾砰擊壓力場(chǎng)的方法．文獻(xiàn)[8]中在拖曳水池中對(duì)某艦船的分段模型進(jìn)行了尾砰擊響應(yīng)試驗(yàn)研究，指出在尾隨浪、零航速、波長(zhǎng)船長(zhǎng)比λ/L=1.0左右時(shí)觀察到了嚴(yán)重的尾砰擊現(xiàn)象．文獻(xiàn)[9]中針對(duì)靜水和波浪場(chǎng)景下的楔形體入水砰擊進(jìn)行試驗(yàn)研究，文獻(xiàn)[10]中針對(duì)平臺(tái)立柱的入水砰擊特性進(jìn)行研究．文獻(xiàn)[11-13]中研究了砰擊場(chǎng)景下多船體甲板結(jié)構(gòu)的動(dòng)態(tài)響應(yīng)問(wèn)題，通過(guò)水彈性梁模型來(lái)模擬甲板，并通過(guò)初始的結(jié)構(gòu)慣性階段和隨后的自由振動(dòng)階段模擬砰擊場(chǎng)景，結(jié)構(gòu)應(yīng)變和位移結(jié)果表明理論與試驗(yàn)結(jié)果吻合較好．文獻(xiàn)[14]中針對(duì)不同船級(jí)社關(guān)于滾裝船外飄砰擊校核公式進(jìn)行對(duì)比分析研究．文獻(xiàn)[15]中利用STAR-CCM+(STAR)對(duì)不同二維和三維剛性和彈性楔體的入水沖擊特性進(jìn)行了研究，認(rèn)為STAR可以較為準(zhǔn)確預(yù)測(cè)楔形體的位移和速度的時(shí)間歷程．

文中使用基于CFD求解器的STAR與有限元軟件Abaqus進(jìn)行交互耦合計(jì)算，對(duì)船體模型[16]在試驗(yàn)水池的砰擊現(xiàn)象進(jìn)行數(shù)值仿真，分析船體發(fā)生砰擊時(shí)的自由液面形態(tài)分布狀況、砰擊壓力、結(jié)構(gòu)變形以及應(yīng)力應(yīng)變分布狀況．

1 波浪砰擊的流固耦合數(shù)值模型

1.1 流體控制方程

STAR軟件基于FVM法將控制方程的積分形式進(jìn)行時(shí)間和空間維度的離散，從而得到可求解的代數(shù)方程組．首先計(jì)算域被細(xì)分為有限數(shù)量的鄰接控制體，這些控制體可以是任意多面體形狀．離散的控制方程在計(jì)算過(guò)程中還需要用到面積積分、體積積分以及時(shí)間和空間導(dǎo)數(shù)．

對(duì)于粘性的三維流動(dòng)，假定流動(dòng)由RANS方程控制，其中湍流效應(yīng)包括渦流模型和粘性模型．此時(shí)就需要求解連續(xù)性方程、動(dòng)量方程和兩個(gè)湍流特性方程．模型選取K-Epsilon湍流模型來(lái)模擬湍流對(duì)流體的影響，其中積分形式的質(zhì)量守恒和動(dòng)量守恒控制方程可以表示為[15]:

(1)

(2)

式中：ρ為流體密度；v為流體速度矢量；vb為控制體表面的速度；n為垂直于控制體表面的單位向量法線(xiàn)，其中S和V分別為控制體的面積和體積；T為指應(yīng)力張量(表示速度梯度和渦粘性)；p為壓力，I為單位張量．

由于VOF模型假定在控制體中各個(gè)不相融流體都具有相同的速度、壓力和溫度場(chǎng)．因此，在控制體中只需要求解單相流場(chǎng)的動(dòng)量、質(zhì)量和能量守恒的基本控制方程即可．描述體積分?jǐn)?shù)αi的輸運(yùn)的守恒方程為:

(3)

式中：sαi為第i相流場(chǎng)的源或匯；Dρi/Dt為相密度ρi的物質(zhì)導(dǎo)數(shù)．

1.2 結(jié)構(gòu)響應(yīng)方程

對(duì)于研究三維艦船在砰擊載荷下的動(dòng)態(tài)響應(yīng)時(shí)，需要利用有限單元法構(gòu)造考慮彈性結(jié)構(gòu)的響應(yīng)方程，將三維模型進(jìn)行網(wǎng)格離散，得到結(jié)構(gòu)的有限元模型，是進(jìn)行有限元計(jì)算的基本步驟．三維彈性動(dòng)力方程的基本方程[17]是在得到船體結(jié)構(gòu)的有限元模型的基礎(chǔ)上，進(jìn)行三維實(shí)體動(dòng)力分析的有限元求解步驟:

(4)

1.3 流固耦合

船體結(jié)構(gòu)波浪砰擊問(wèn)題屬于強(qiáng)非線(xiàn)性的流固耦合問(wèn)題，文中通過(guò)STAR與Abaqus之間的協(xié)同耦合功能，自動(dòng)在流固交界面處進(jìn)行數(shù)據(jù)的交互傳遞實(shí)現(xiàn)雙向的流固耦合計(jì)算．其中STAR求解流體方程，Abaqus求解結(jié)構(gòu)方程，屬于分區(qū)式流固耦合．采用STAR和Abaqus之間通過(guò)隱式的耦合方法進(jìn)行耦合計(jì)算，主要是因?yàn)殡[式耦合方法適用于結(jié)構(gòu)和流體之間的強(qiáng)耦合應(yīng)用，它允許STAR和Abaqus在每一時(shí)間步內(nèi)交換不止一次數(shù)據(jù)，不停迭代直到滿(mǎn)足某些收斂準(zhǔn)則．這種耦合方案更穩(wěn)定，二階精度更高．

由于波浪砰擊現(xiàn)象屬于強(qiáng)耦合現(xiàn)象，因此在數(shù)值模型中，選擇隱式強(qiáng)耦合算法．文中依單個(gè)耦合步驟為例，針對(duì)隱式耦合時(shí)間步的流程進(jìn)行介紹，具體如圖1．在STAR中設(shè)置運(yùn)動(dòng)規(guī)范為“變形”，在Abaqus設(shè)置模型的彈性模量、板厚、初始速度、重力、時(shí)間步長(zhǎng)和耦合面等參量．具體計(jì)算流程如下:首先STAR開(kāi)始流場(chǎng)載荷計(jì)算，然后將壓力和剪切力傳遞給FEM求解器，之后Abaqus根據(jù)耦合界面的載荷進(jìn)行結(jié)構(gòu)響應(yīng)分析并將得到的位移和變形量傳遞給STAR，STAR再根據(jù)傳遞的位移量通過(guò)“變形”功能實(shí)現(xiàn)網(wǎng)格的移動(dòng)，這樣完成一次協(xié)同耦合進(jìn)行計(jì)算更新，重復(fù)直至達(dá)到最大物理時(shí)間，tn為計(jì)算時(shí)間步．

圖1 隱式算法的單個(gè)耦合時(shí)間步的流程Fig.1 Procedure for a single coupling step with implicit solvers

2 模型設(shè)置

2.1 模型尺寸

文中建立的波浪砰擊載荷下船體結(jié)構(gòu)動(dòng)態(tài)響應(yīng)的數(shù)值仿真模型和文獻(xiàn)[16]中的鋁制模型船試驗(yàn)作為對(duì)比，文獻(xiàn)模型的主尺度參數(shù)如表1，船體的有限元模型如圖2．

表1 主尺度參數(shù)Table 1 Principal dimension

圖2 鋁船的有限元模型Fig.2 FEM model of the aluminum ship

2.2 計(jì)算域設(shè)置及網(wǎng)格劃分情況

坐標(biāo)系原點(diǎn)在船尾垂線(xiàn)與船底平面的交點(diǎn)處，X軸沿船長(zhǎng)方向指向船首，Z軸垂直向上，并滿(mǎn)足右手直角坐標(biāo)系．計(jì)算域如圖3，其中背景區(qū)域(background)坐標(biāo)為[-8，-4，-4]，[24， 4， 4]，重疊區(qū)域(overset)坐標(biāo)為[-1.5，-1.2，-0.5]，[5， 1.2， 0.8]，網(wǎng)格的總數(shù)為1 178 788個(gè)．模擬收斂精度為98.779%．計(jì)算域背景區(qū)域[18]和重疊區(qū)域的邊界條件設(shè)置可參見(jiàn)表2．

圖3 計(jì)算域Fig.3 Computational domain

表2 邊界條件Table 2 Boundary condition

整體計(jì)算域在Y=0處剖面網(wǎng)格劃分如圖4，最大網(wǎng)格尺寸為0.2 m，加密區(qū)域網(wǎng)格尺寸為0.1 m，切割體網(wǎng)格單元為各項(xiàng)同性．另外，為捕捉自由液面處波浪形態(tài)和存在的波浪翻卷及破碎現(xiàn)象，對(duì)船體表面鄰近區(qū)域的網(wǎng)格進(jìn)行加密，如圖4，水面附近網(wǎng)格Z方向加密，Z方向尺寸0.025 m，重疊區(qū)域網(wǎng)格尺寸為0.1 m，船表面尺寸為0.05 m，棱柱層厚度為0.025 m．

圖4 液面處網(wǎng)格加密區(qū)域Fig.4 Section meshing at the free surface domain

文中在Abaqus中設(shè)置試驗(yàn)?zāi)Ｐ偷陌搴?、彈性模量等材料屬性，模型使? mm厚度的鋁板，并對(duì)結(jié)構(gòu)施加初始速度和重 力，導(dǎo)出inp文件．在STAR的“協(xié)同仿真”中確定導(dǎo)入位移和速度，導(dǎo)出壓力和壁面剪切應(yīng)力，讀取Abaqus中模型的inp文件，實(shí)現(xiàn)STAR與Abaqus之間的數(shù)據(jù)交換．設(shè)置區(qū)域中重疊區(qū)域(overset)的運(yùn)動(dòng)規(guī)范為“變形”．

文中基于建立的波浪砰擊耦合數(shù)值模型，研究了規(guī)則波迎浪情況下鋁制模型船的砰擊場(chǎng)景及其船體結(jié)構(gòu)響應(yīng)．選取船體無(wú)航速(工況1)和有航速(工況2)進(jìn)行數(shù)值仿真，分析結(jié)構(gòu)運(yùn)動(dòng)及砰擊響應(yīng)特性．選取的兩組工況如表3．

表3 仿真選取的兩種工況Table 3 Two working conditions selected by simulation

為保證船體波浪砰擊耦合分析過(guò)程中數(shù)值模型的正確性，將數(shù)值結(jié)果和文獻(xiàn)結(jié)果在砰擊載荷進(jìn)行對(duì)比，具體如圖5，t為試驗(yàn)及仿真計(jì)算時(shí)刻，圖中給出了數(shù)值和試驗(yàn)關(guān)于支桿與船體連接處的垂向載荷的對(duì)比，從圖中可以看出，仿真結(jié)果與文獻(xiàn)值整體曲線(xiàn)趨勢(shì)一致，吻合度較高，仿真載荷幅值稍小于文獻(xiàn)值，這可能是由于波浪試驗(yàn)是在水域有限的水池中進(jìn)行的，壁面對(duì)波的反射可能會(huì)造成試驗(yàn)結(jié)果的偏大，而且工況1的整體誤差稍大于工況2的情況，這可能是由于工況1的波高比工況2更大，導(dǎo)致文獻(xiàn)中試驗(yàn)結(jié)果與仿真之間的差值也越大，另一方面，數(shù)值計(jì)算誤差也對(duì)精度有一定的影響．

圖5 垂向力的試驗(yàn)值和仿真結(jié)果的對(duì)比Fig.5 Comparison of experimental and numerical values of heave force

根據(jù)圖5的載荷對(duì)比曲線(xiàn)可知，文中所提出的波浪砰擊下船體耦合響應(yīng)數(shù)值仿真方法是有效的．

3 模型船體砰擊壓力

3.1 砰擊時(shí)刻自由液面形態(tài)和模型壓力云圖

圖6為鋁制模型船在工況1時(shí)發(fā)生砰擊的前后時(shí)刻情況，當(dāng)船舶漂浮在波浪上時(shí)，由于波浪的周期性起伏運(yùn)動(dòng)，船體也隨之出現(xiàn)升沉和縱搖現(xiàn)象，圖中主要給出了不同時(shí)刻模型自由液面形態(tài)及船底表面砰擊壓力云圖．圖6(a)為t=1.48 s時(shí)自由液面形態(tài)和模型表面壓力云圖．從圖中可知，該時(shí)刻波峰位于模型中部偏船尾位置，模型首部位于波谷處，首部表面完全出水，壓力集中于船底偏尾部區(qū)域，首部壓力最小．由圖6(b)可知，t=2.06 s時(shí)模型尾部即將進(jìn)入波谷，船首位于波浪波峰處，船首底部?jī)A斜處發(fā)生波浪沖擊，模型船首表面與波浪水面完全接觸，最大壓力集中在首部．該工況時(shí)波浪表面未發(fā)生大的表面變形甚至波浪破碎現(xiàn)象．

圖6 工況1船體各時(shí)刻自由液面形態(tài)和表面壓力云圖Fig.6 Free surface elevation and the pressure distribution on the ship model in case 1

圖7為鋁制模型船在工況2時(shí)發(fā)生砰擊前后時(shí)刻的模型自由液面形態(tài)及船底表面砰擊壓力云圖．

圖7 工況2船體各時(shí)刻自由液面形態(tài)和表面壓力云圖Fig.7 Free surface elevation and the pressure distribution on the ship model in case 2

由圖7(a)可知，在t=1.25 s時(shí)波峰處于模型后方近船尾處，模型首部位于波谷處，首部表面完全出水，壓力主要集中在船底偏尾部，而船首壓力較?。蓤D7(b)可知，在t=1.56 s時(shí)刻波峰運(yùn)動(dòng)到船首處，船尾處離開(kāi)上一個(gè)波峰位置即將進(jìn)入波谷，模型首部?jī)A斜處發(fā)生相對(duì)波浪沖擊，波浪表面沿模型表面出現(xiàn)射流攀升現(xiàn)象，最大壓力出現(xiàn)在模型首部?jī)A斜位置．

對(duì)比圖6、7可知，在工況2時(shí)，由于模型具有向前運(yùn)動(dòng)的航速，使得船體與波浪之間存在較大的相對(duì)速度，造成模型在較小波高的情況下，船首傾斜處與波浪水面交界處出現(xiàn)了液面的猛烈沖擊現(xiàn)象，甚至出現(xiàn)波浪破碎射流分離的強(qiáng)非線(xiàn)性現(xiàn)象，造成船體首部出現(xiàn)較大的砰擊載荷．

3.2 測(cè)點(diǎn)砰擊壓力和速度對(duì)比曲線(xiàn)

如圖8，在船底布置一系列測(cè)點(diǎn)以探測(cè)船體在波浪中受到的砰擊壓力及結(jié)構(gòu)變形位移．P2位于船首底部轉(zhuǎn)折處，坐標(biāo)為(3.4， 0， 0);P4點(diǎn)位于船尾底部轉(zhuǎn)折處，坐標(biāo)為(0.4， 0， 0)．

圖8 測(cè)點(diǎn)位置示意Fig.8 Schematic map of measuring point position

圖9為工況1時(shí)不同測(cè)點(diǎn)的壓力和速度曲線(xiàn)歷程曲線(xiàn)，其中實(shí)線(xiàn)為砰擊壓力曲線(xiàn)，虛線(xiàn)為速度曲線(xiàn)，左側(cè)Y軸為壓力值，右側(cè)Y軸為速度值．

圖9 工況1測(cè)點(diǎn)壓力和速度對(duì)比曲線(xiàn)Fig.9 Correlation curve of pressure and velocity at measuring point in case 1

從圖中可知，壓力p和速度v曲線(xiàn)均呈現(xiàn)周期性變化，反映了砰擊載荷的頻繁出現(xiàn)的特性．并且圖中可知P2測(cè)點(diǎn)(船首底部)的壓力幅值最大為1.64 kPa，P4測(cè)點(diǎn)(船尾底部)幅值為0.90 kPa，船體首尾處均相繼出現(xiàn)砰擊現(xiàn)象．圖中壓力曲線(xiàn)峰值稍有一定的衰減，可能是由于湍流模型存在粘度以及阻尼造成波高衰減造成的，而且P2測(cè)點(diǎn)壓力曲線(xiàn)呈現(xiàn)正弦函數(shù)形式，與規(guī)則波的輸入相符，此時(shí)砰擊壓力曲線(xiàn)與速度曲線(xiàn)峰值向?qū)?yīng)，表明了砰擊載荷與砰擊速度的相關(guān)性，其余測(cè)點(diǎn)的壓力幅值與速度幅值存在一定的相位差，且測(cè)點(diǎn)的壓力曲線(xiàn)均較為規(guī)則．

圖10為工況2時(shí)測(cè)點(diǎn)的壓力和速度曲線(xiàn)對(duì)比，對(duì)比圖10、9，可以發(fā)現(xiàn)兩種工況下相應(yīng)各測(cè)點(diǎn)整體曲線(xiàn)特征相似，工況2時(shí)船首區(qū)域測(cè)點(diǎn)壓力峰值為3.294 kPa，大于工況1中對(duì)應(yīng)測(cè)點(diǎn)峰值，這主要是因?yàn)楣r2中增加了船舶航速，使得砰擊現(xiàn)象較工況1更加明顯．工況2中P4測(cè)點(diǎn)均出現(xiàn)負(fù)壓載荷，這是由測(cè)點(diǎn)所處位置的流場(chǎng)載荷的往復(fù)變化以及結(jié)構(gòu)的變形效應(yīng)共同作用下的結(jié)果．P2測(cè)點(diǎn)位于船艏底部，周?chē)鲌?chǎng)呈規(guī)律性變化，而P4測(cè)點(diǎn)位于船艉，且周?chē)鲌?chǎng)由于受到波浪破碎的原因，導(dǎo)致壓力分布特性變化較大．

圖10 工況2測(cè)點(diǎn)壓力和速度對(duì)比曲線(xiàn)Fig.10 Correlation curve of pressure and velocity at measuring point in case 2

3.3 船體結(jié)構(gòu)響應(yīng)

圖11為工況1時(shí)不同時(shí)刻的模型結(jié)構(gòu)應(yīng)力(單位:MPa)和應(yīng)變?cè)茍D.圖12為工況2時(shí)不同時(shí)刻的模型結(jié)構(gòu)應(yīng)力和應(yīng)變?cè)茍D．

圖11 工況1模型應(yīng)力和應(yīng)變?cè)茍DFig.11 Stress and strain nephogram of case 1

圖11為工況1下t=1.48 s時(shí)模型應(yīng)力和應(yīng)變?cè)茍D，由圖6可知，t=1.48 s時(shí)刻船體中后部處于波峰位置，船首處完全出水，船體出現(xiàn)中拱現(xiàn)象，對(duì)應(yīng)圖11中船體底部及舷側(cè)出現(xiàn)大的應(yīng)力和應(yīng)變集中，而船首處應(yīng)力、應(yīng)變較?。藭r(shí)的對(duì)應(yīng)的最大應(yīng)力值為2.2 Mpa，最大應(yīng)變值為3.053×10-5．

圖12為工況2時(shí)船體受波浪砰擊前后時(shí)刻應(yīng)力、應(yīng)變?cè)茍D，對(duì)應(yīng)圖12、7可知，工況2中t=1.56 s時(shí)船首處發(fā)生波浪砰擊現(xiàn)象，此時(shí)船首底部與波浪發(fā)生猛烈砰擊，出現(xiàn)波浪破碎的情況，因此造成船首應(yīng)力、應(yīng)變?cè)龃?，出現(xiàn)相應(yīng)的應(yīng)力、應(yīng)變集中區(qū)，船體尾部區(qū)域的應(yīng)力明顯較?。藭r(shí)的對(duì)應(yīng)的最大應(yīng)力值為0.683 MPa，最大應(yīng)變值為7.34×10-6．

圖12 工況2模型應(yīng)力和應(yīng)變?cè)茍DFig.12 Stress and strain nephogram of case 2

4 結(jié)論

基于文獻(xiàn)中的模型試驗(yàn)進(jìn)行了仿真分析，通過(guò)對(duì)比試驗(yàn)測(cè)量值和數(shù)值仿真結(jié)果，檢驗(yàn)文中建立的波浪砰擊載荷下船體結(jié)構(gòu)動(dòng)態(tài)響應(yīng)的數(shù)值仿真模型的可行性．分析船體的發(fā)生砰擊時(shí)的自由液面形態(tài)分布狀況、砰擊壓力、結(jié)構(gòu)變形以及應(yīng)力應(yīng)變分布狀況．

(1) 所建立的計(jì)及結(jié)構(gòu)變形效應(yīng)的波浪砰擊流固耦合仿真方法與文獻(xiàn)中試驗(yàn)結(jié)果整體趨勢(shì)吻合較好，只是幅值相對(duì)試驗(yàn)稍小，這可能是由于數(shù)值模擬與實(shí)船試驗(yàn)場(chǎng)景的差異如忽略了水池實(shí)驗(yàn)的壁面效應(yīng)以及數(shù)值誤差的影響造成的．因此文中建立的數(shù)值模型對(duì)于模擬船舶在波浪中的運(yùn)動(dòng)及砰擊現(xiàn)象具有較大的適用性和可行性．

(2) 當(dāng)船體在波浪中運(yùn)動(dòng)時(shí)，船首區(qū)域由于波浪的沖擊會(huì)發(fā)生砰擊現(xiàn)象，且當(dāng)船體有航速時(shí)自由液面出現(xiàn)了波浪表面破碎等強(qiáng)非線(xiàn)性的砰擊現(xiàn)象，并分析了在不同航行階段中船體結(jié)構(gòu)應(yīng)力和應(yīng)變的變化特征，觀察到船底應(yīng)力變化，表明船首底部是波浪砰擊問(wèn)題中需要重要關(guān)注的區(qū)域．