徐國(guó)徽，顧學(xué)康

（中國(guó)船舶科學(xué)研究中心，江蘇 無(wú)錫 214082）

1 引 言

近年來(lái)隨著液化天然氣（LNG）需求的增加，液化天然氣船船體和液艙主尺度不斷大型化。大型化后的LNG船載貨量一般為25-35萬(wàn)立方，幾乎為傳統(tǒng)型LNG船的兩倍。由于吃水限制，液艙主要增加長(zhǎng)度和寬度。因而船體運(yùn)動(dòng)引起的液化天然氣的晃蕩運(yùn)動(dòng)及其對(duì)液艙結(jié)構(gòu)的沖擊問(wèn)題越來(lái)越突出，部分裝載的需求也越來(lái)越高[1]。液艙部分裝載時(shí)，船體大幅晃蕩引起的沖擊壓力對(duì)艙室結(jié)構(gòu)和圍護(hù)系統(tǒng)可能造成嚴(yán)重破壞，此問(wèn)題已經(jīng)引起了研究和設(shè)計(jì)者們的廣泛關(guān)注。另外，由于液艙邊界的彈性圍護(hù)系統(tǒng)的存在，大型LNG船艙內(nèi)液貨大幅晃蕩運(yùn)動(dòng)時(shí)作用在艙壁上的沖擊壓力將受到艙壁彈性的影響，導(dǎo)致晃蕩載荷和結(jié)構(gòu)響應(yīng)與剛性艙壁情況會(huì)有所不同。因此就液艙晃蕩中的流固耦合問(wèn)題開(kāi)展數(shù)值計(jì)算方法的研究具有明顯的工程實(shí)用意義，國(guó)內(nèi)外已有學(xué)者展開(kāi)了相關(guān)問(wèn)題的探索性研究。

朱仁慶[2]采用自編程序?qū)⒘鲌?chǎng)通度的概念、VOF法與水彈性力學(xué)理論相結(jié)合，進(jìn)行了二維粘性流體與彈性結(jié)構(gòu)耦合作用的理論與分析方法研究。Zhu等[3]應(yīng)用處理粘性流體晃蕩與彈性結(jié)構(gòu)相互耦合作用的理論及相應(yīng)數(shù)值計(jì)算方法，流體運(yùn)動(dòng)采用N-S方程描述，控制方程采用有限差分法離散，并由超松弛迭代法求解，液體自由表面通過(guò)流體體積法進(jìn)行重構(gòu)，考察二維液艙不同結(jié)構(gòu)剛度對(duì)液體晃蕩的影響。娜日薩[4]在DNV基于壓力的液艙晃蕩仿真及結(jié)構(gòu)強(qiáng)度評(píng)估方法的基礎(chǔ)上，引入時(shí)空放大、體積模擬減縮等技術(shù)提出了耦合的VLCC液艙晃蕩仿真及強(qiáng)度評(píng)估方法。Rognebakke和Faltinsen[5]考慮液體內(nèi)空泡和結(jié)構(gòu)的水彈性效應(yīng)，采用數(shù)值方法研究了部分充滿的矩形罐的誘捕空氣（entrapment of air）和水彈性問(wèn)題。Wang和Kim[6]采用聲固耦合分析法和簡(jiǎn)單的二維聲固耦合模型研究結(jié)構(gòu)響應(yīng)的水彈性效應(yīng)。Hisashi等[7]結(jié)合ALE方法做局部沖擊分析，計(jì)算貨艙系統(tǒng)的結(jié)構(gòu)響應(yīng)。目前液艙流固耦合問(wèn)題的研究還處于初步階段，學(xué)者們的方法基于不同的條件假設(shè)和理論，得出的結(jié)果和論斷也各有差異。

本文研究大型LNG船液艙晃蕩中的流固耦合效應(yīng)，采用基于顯式時(shí)間積分方法的有限元/有限體積法程序DYTRAN模擬流體和結(jié)構(gòu)的運(yùn)動(dòng)和變形及其相互作用，從二維剛性模型開(kāi)始，以計(jì)算流體力學(xué)方法結(jié)果和試驗(yàn)結(jié)果作為參照進(jìn)行數(shù)值方法研究和參數(shù)分析，進(jìn)而建立三維剛性和彈性液艙兩種模型，比較不同模型的數(shù)值計(jì)算結(jié)果，得到了對(duì)大型LNG船液艙晃蕩載荷和強(qiáng)度評(píng)估具有參考價(jià)值的結(jié)論。

2 耦合理論介紹

依據(jù)耦合機(jī)理，如果耦合作用僅僅發(fā)生在兩相域的交界面上，在方程上的耦合是由兩相耦合面的平衡和協(xié)調(diào)關(guān)系引入的。對(duì)于這類(lèi)問(wèn)題按其兩相間相對(duì)運(yùn)動(dòng)的大小及相互作用的性質(zhì)又可分為三類(lèi)[8]:一是流體和固體結(jié)構(gòu)之間有較大的相對(duì)運(yùn)動(dòng)。如機(jī)翼顫振或懸橋震蕩中發(fā)生的氣動(dòng)耦合問(wèn)題，習(xí)慣上稱(chēng)為氣動(dòng)彈性力學(xué)問(wèn)題；二是具有流體有限位移的短期問(wèn)題。流固相互作用是在瞬間完成的，其特征是流體的密度發(fā)生急劇變化，如液體中的爆炸與沖擊問(wèn)題；三是液體域與固體域的動(dòng)力相互作用。其特點(diǎn)是流固相互作用時(shí)間較長(zhǎng)，相對(duì)位移有限，習(xí)慣上稱(chēng)為液動(dòng)彈性力學(xué)問(wèn)題，如近海結(jié)構(gòu)對(duì)波浪或地震的響應(yīng)、充液容器的晃動(dòng)問(wèn)題、壩水、船水耦合響應(yīng)以及水力機(jī)械在流體介質(zhì)中的振動(dòng)問(wèn)題；本文所研究的晃蕩問(wèn)題就屬于第三類(lèi)問(wèn)題。

流固耦合問(wèn)題按其研究方法可分為弱耦合法和強(qiáng)耦合法，或分別稱(chēng)之為分區(qū)迭代求解法和直接求解法。所謂弱耦合法是對(duì)流體域和固體域分別建立和求解運(yùn)動(dòng)方程，兩個(gè)運(yùn)動(dòng)方程間則通過(guò)插值函數(shù)進(jìn)行流體壓力、結(jié)構(gòu)位移等的數(shù)值傳遞，從而實(shí)現(xiàn)流體和結(jié)構(gòu)的相互耦合。采用弱耦合方法時(shí)，流體和結(jié)構(gòu)用各自的求解器在時(shí)域上積分，交錯(cuò)時(shí)間推進(jìn)，這種方法能充分利用現(xiàn)有的計(jì)算流體動(dòng)力學(xué)和計(jì)算結(jié)構(gòu)動(dòng)力學(xué)的方法和程序，只需作少量的修改，從而保持各程序的模塊化特性。它的主要缺點(diǎn)是由于積分是交錯(cuò)進(jìn)行的，流體和結(jié)構(gòu)的時(shí)間推進(jìn)積分總是存在時(shí)間滯后，耦合界面上的能量不能保持守恒，無(wú)法滿足動(dòng)平衡。強(qiáng)耦合法則是對(duì)流體和結(jié)構(gòu)分別建立運(yùn)動(dòng)方程，采用數(shù)值方法將兩個(gè)方程耦合起來(lái)形成同一控制方程即耦合方程，在單一時(shí)間步內(nèi)，同時(shí)求解耦合界面的壓力分布以及結(jié)構(gòu)對(duì)壓力做出的響應(yīng)，避免了弱耦合法中出現(xiàn)的時(shí)間滯后問(wèn)題。顯然，強(qiáng)耦合法能夠準(zhǔn)確地描述流體的運(yùn)動(dòng)，更真實(shí)地反映流體和結(jié)構(gòu)的相互作用，理論上更加合理，分析過(guò)程更加嚴(yán)密。但由于強(qiáng)耦合理論的復(fù)雜性以及計(jì)算機(jī)軟硬件資源的限制，該方法適合求解沒(méi)有接觸分析的小型到中型的瞬態(tài)耦合問(wèn)題?；问帉儆诘湫蛷?qiáng)耦合問(wèn)題，需要同時(shí)求解流體運(yùn)動(dòng)和結(jié)構(gòu)的響應(yīng)。

流固耦合問(wèn)題的計(jì)算方法應(yīng)用較多的有三類(lèi):

（1）固體和流體均采用修正拉格朗日法

此法優(yōu)點(diǎn)是編程方便，即用一種方法可對(duì)流固兩種介質(zhì)進(jìn)行分析。對(duì)不可壓流體相關(guān)的流固耦合問(wèn)題，此類(lèi)解法相當(dāng)簡(jiǎn)便。商業(yè)軟件多采用此方法處理簡(jiǎn)單的不可壓流體與固體的相互作用問(wèn)題。缺點(diǎn)是對(duì)可壓縮流場(chǎng)相關(guān)的流固耦合問(wèn)題，尤其是當(dāng)流場(chǎng)中出現(xiàn)旋渦或強(qiáng)間斷時(shí)，該方法將因網(wǎng)格變形過(guò)大或畸變，引起計(jì)算時(shí)間步長(zhǎng)銳減，甚至直接導(dǎo)致計(jì)算失敗。

（2）流體采用修正的歐拉方法處理，固體采用拉格朗日方法處理

此法屬于混合法，充分發(fā)揮了歐拉方法高精度捕捉流場(chǎng)與拉格朗日方法準(zhǔn)確描述結(jié)構(gòu)變形的優(yōu)點(diǎn)，但在流固界面（以下稱(chēng)為內(nèi)界面）的數(shù)值交換提出新的邊界條件，流場(chǎng)壁面邊界節(jié)點(diǎn)處的物理量要通過(guò)插值得到，對(duì)內(nèi)邊界條件處理的精度直接影響流場(chǎng)求解的精度。

（3）流體采用ALE方法處理，固體采用拉格朗日方法處理

此法屬于混合法，采用ALE網(wǎng)格，不需要進(jìn)行網(wǎng)格重分及物理量插值，通過(guò)流場(chǎng)網(wǎng)格移動(dòng)保證在內(nèi)邊界處流場(chǎng)與結(jié)構(gòu)網(wǎng)格的銜接。流場(chǎng)內(nèi)部的網(wǎng)格節(jié)點(diǎn)運(yùn)動(dòng)引起單元體積的變化滿足幾何守恒（Geometric Conservation Laws,GCL），內(nèi)邊界求解由ALE算法保證精度不降階。

在DYTRAN程序中同時(shí)提供Lagrange和Euler求解器，既模擬結(jié)構(gòu)的變形又模擬材料的流動(dòng)，Lagrange網(wǎng)格與Euler網(wǎng)格之間進(jìn)行耦合，從而分析流體和結(jié)構(gòu)之間的相互作用。DYTRAN程序提供后兩種耦合算法，一般耦合（general coupling）屬于第二類(lèi)，ALE耦合（Arbitrary Lagrange-Euler coupling）屬于第三類(lèi)。一般耦合方法對(duì)規(guī)則歐拉網(wǎng)格提出的快速算法有較高的計(jì)算速度由于它的耦合面不需要?dú)W拉節(jié)點(diǎn)與拉格朗日節(jié)點(diǎn)的一一對(duì)應(yīng)，所以對(duì)于結(jié)構(gòu)幾何形狀復(fù)雜的流固耦合問(wèn)題，一般耦合建模工作較為簡(jiǎn)單。但耦合交界面位置的計(jì)算及單元的合并增加了計(jì)算時(shí)間。ALE方法由于歐拉域的邊界是已經(jīng)確定的，這避免了像一般耦合那樣不斷計(jì)算耦合面位置的過(guò)程，也可以避免為實(shí)現(xiàn)耦合而建立過(guò)多的歐拉單元帶來(lái)的計(jì)算量的增加，大大提高了計(jì)算的效率。ALE方法不需要流體單元的不斷合并，數(shù)值計(jì)算結(jié)果中的噪聲較少。對(duì)于像液艙晃蕩這種結(jié)構(gòu)形式簡(jiǎn)單，固定內(nèi)部流體的流固耦合問(wèn)題，ALE方法是比較適合的。這也是本文研究液艙晃蕩問(wèn)題采用ALE法的原因。

3 數(shù)值模擬方法研究

3.1 二維剛性模型描述

為了便于比較分析，本文選取作者前期研究[8]的矩形液艙模型（圖1）及計(jì)算工況（表1）作為算例，艙室橫蕩激勵(lì)為簡(jiǎn)諧正弦函數(shù)，以艙室?guī)缀沃行臑樵c(diǎn)。采用任意拉格朗日（Lagrange）－歐拉（Euler）耦合（ALE）算法，流體域采用歐拉六面體單元?jiǎng)澐?，液艙結(jié)構(gòu)采用殼單元?jiǎng)澐?，兩者網(wǎng)格密度一致，流固耦合面上兩種網(wǎng)格節(jié)點(diǎn)空間重合，歐拉網(wǎng)格不再是空間固定而是隨著結(jié)構(gòu)的變形而移動(dòng)。耦合典型過(guò)程是:（a）歐拉材料流動(dòng)引起的壓力載荷作用到結(jié)構(gòu)網(wǎng)格上；（b）在壓力作用下對(duì)結(jié)構(gòu)響應(yīng)進(jìn)行積分；（c）把結(jié)構(gòu)的邊界位移和運(yùn)動(dòng)傳遞給流體；（d）更新流體動(dòng)態(tài)網(wǎng)格；（e）流體積分，計(jì)算新的流體壓力和應(yīng)力場(chǎng)。此耦合法避免了結(jié)構(gòu)流體兩種網(wǎng)格之間的交叉，減少了單元合并工作量，降低了結(jié)果的數(shù)值噪聲。

表1 計(jì)算工況Tab.1 Calculation conditions

圖1 艙室?guī)缀文Ｐ秃蛪毫y(cè)點(diǎn)位置(單位mm)Fig.1 Geometrical model of the tank and pressure measuring point arrangement(mm)

圖2所示的二維模型液艙周壁為剛性材料，引用MATRIG卡片可以把Lagrange網(wǎng)格定義為剛體；艙內(nèi)水為無(wú)粘、可壓縮線性流體本構(gòu)關(guān)系材料，引用EOSPOL卡片定義；通過(guò)定義剛體質(zhì)心處的強(qiáng)迫運(yùn)動(dòng)給以運(yùn)動(dòng)激勵(lì)，在流體的非主要運(yùn)動(dòng)方向設(shè)置一個(gè)單元，這樣流體僅在建模平面內(nèi)運(yùn)動(dòng)。

DYTRAN的材料模式中，包括了線彈性、彈塑性、剛性材料、橡膠材料、低密度泡沫材料、土壤材料、正交各向異性材料、層合復(fù)合材料、率相關(guān)材料以及各種屈服準(zhǔn)則、失效模式、狀態(tài)方程、多點(diǎn)爆炸燃燒模型等。

流體域采用線性流體來(lái)填充，多項(xiàng)式狀態(tài)方程EOSPOL定義水的多項(xiàng)式狀態(tài)方程，這里不計(jì)水的粘性，壓力是相對(duì)體積與比內(nèi)能的多項(xiàng)式函數(shù):

壓縮狀態(tài)（μ＞0）

圖2 二維剛性模型Fig.2 2D rigid model

拉伸狀態(tài)（μ≤0）

其中，μ=（ρ/ρ0）-1；ρ為水的密度，ρ0為初始狀態(tài)水的密度；e為單位質(zhì)量的比內(nèi)能。

忽略非線性項(xiàng)后可得:

其中，a1為水的體積模量。在計(jì)算的初始時(shí)刻，需定義歐拉單元的初始狀態(tài)。

對(duì)于彈性液艙的問(wèn)題，所用的結(jié)構(gòu)材料包括剛性材料和線彈性材料。剛體用RIGID材料，用MATRIG卡進(jìn)行定義，該卡將一部分Lagrange單元定義為剛體，這一部分單元為一個(gè)剛性整體，單元之間沒(méi)有相對(duì)移動(dòng)，單元本身也沒(méi)有彈性變形。該方法還可以通過(guò)定義任意幾何面為剛體，不計(jì)單元與幾何的形狀，其重量、重心和慣性矩也可以定義。

彈性體本構(gòu)模型為彈性材料DMATEP，DMATEP卡只需要定義參考密度和4個(gè)彈性常數(shù)（彈性模量E、泊松比μ、體積模量K和剪切模量G）中的任意兩個(gè)就行，他們之間有如下關(guān)系:

3.2 網(wǎng)格更新方法

ALE耦合法的一個(gè)優(yōu)點(diǎn)就是流體網(wǎng)格獨(dú)立于材料變形而進(jìn)行空間運(yùn)動(dòng)，可以根據(jù)需要自由選擇其運(yùn)動(dòng)狀態(tài)，對(duì)數(shù)值分析物體的變形過(guò)程帶來(lái)了便利。不同于計(jì)算流體力學(xué)（CFD）里處理流體動(dòng)邊界的動(dòng)網(wǎng)格法和自由表面跟蹤技術(shù)VOF法，DYTRAN采用ALE方法描述帶自由液面的晃動(dòng)問(wèn)題時(shí)，處理計(jì)算區(qū)域內(nèi)動(dòng)態(tài)的網(wǎng)格與流體運(yùn)動(dòng)之間的域耦合，主要有以更新網(wǎng)格位置為目標(biāo)和以獲得網(wǎng)格節(jié)點(diǎn)運(yùn)動(dòng)速度為目標(biāo)的兩類(lèi)網(wǎng)格節(jié)點(diǎn)速度更新方法。

Hughes（1981）等[10]提出 L-E（Lagrangian-Eulerian）矩陣法，通過(guò)對(duì)單元節(jié)點(diǎn)定義不同的 L-E系數(shù)與相應(yīng)的流場(chǎng)質(zhì)點(diǎn)速度,共同確定網(wǎng)格內(nèi)部節(jié)點(diǎn)的運(yùn)動(dòng)；Donea（1982）等[11]提出了平均法，即用上一時(shí)刻相鄰網(wǎng)格點(diǎn)的速度平均值確定內(nèi)部網(wǎng)格節(jié)點(diǎn)速度；Huerta（1988）等[12]提出了混合法，即只在自由液面或運(yùn)動(dòng)邊界上進(jìn)行L-E網(wǎng)格更新，內(nèi)部網(wǎng)格節(jié)點(diǎn)速度用勢(shì)流方程或代數(shù)方法計(jì)算得到。目前ALE網(wǎng)格更新方法都是以上述方法為基礎(chǔ)，其中混合法應(yīng)用最廣。合理高效的網(wǎng)格算法可以在自由面附近獲得性態(tài)良好的計(jì)算單元，降低大幅度液體晃動(dòng)計(jì)算的網(wǎng)格重構(gòu)頻率，同時(shí)網(wǎng)格節(jié)點(diǎn)速度的數(shù)值還將影響到控制方程中對(duì)流項(xiàng)的量級(jí)，更重要的是改進(jìn)自由液面網(wǎng)格節(jié)點(diǎn)速度的求解方案，能使ALE方法在更多樣的幾何邊界和運(yùn)動(dòng)邊界問(wèn)題中得到廣泛應(yīng)用。DYTRAN程序中在ALEGRID卡片中有一些算法，對(duì)不同問(wèn)題可能需要不同的網(wǎng)格運(yùn)動(dòng)形式才適宜，否則更新過(guò)程中就會(huì)發(fā)生網(wǎng)格畸變。

3.3 數(shù)值噪聲控制

在動(dòng)力學(xué)數(shù)值分析中，數(shù)值噪聲對(duì)計(jì)算結(jié)果有很大的影響，控制數(shù)值噪聲非常重要。設(shè)置人工粘性是為了保持計(jì)算結(jié)果的平穩(wěn)性。體積粘性針對(duì)材料性質(zhì)，用于控制振蕩波的形成；沙漏阻尼針對(duì)單元特性，用于避免低精度單點(diǎn)積分體殼單元的無(wú)剛度變形模態(tài)，提高計(jì)算精度。

3.4 時(shí)間步長(zhǎng)控制

顯式求解方法要保證解的穩(wěn)定性，求解的時(shí)間步長(zhǎng)是隱式求解的時(shí)間步長(zhǎng)的1/1000～1/100，即時(shí)間步長(zhǎng)Δt必需小于應(yīng)力波跨越網(wǎng)格中最小單元的時(shí)間:Δt=min（）。臨界時(shí)間步長(zhǎng)Δt由單元的特cr征長(zhǎng)度和單元的材料特點(diǎn)決定，不同類(lèi)型的單元有其對(duì)應(yīng)的算法:

圖3 三種網(wǎng)格方案計(jì)算的P4點(diǎn)的壓力時(shí)間歷程（60%H）Fig.3 Simulated time histories of pressure P4 for three different grid schemes(60%H)

圖4 CFD計(jì)算的P4點(diǎn)壓力時(shí)間歷程（60%H）Fig.4 Time history of pressure P4 predicted by CFD(60%H)

上式中，Le是單元的特征長(zhǎng)度；c是聲音在材料中的傳播速度，E是彈性模量，ρ是質(zhì)量密度，σ是泊松系數(shù)，u是流體材料的速度，K是體積模量。在物理?xiàng)l件一定的情況下，時(shí)間步長(zhǎng)由最小單元決定，所以在網(wǎng)格中盡量避免很小的單元是至關(guān)重要的。 這里采用 30×15、60×30、120×60 三套網(wǎng)格方案，歐拉單元特征長(zhǎng)度分別為單位長(zhǎng)度0.04 m、0.02 m和0.01 m。圖3為幾種網(wǎng)格的計(jì)算結(jié)果比較，粗網(wǎng)格計(jì)算速度快，但誤差造成的數(shù)值振蕩比較明顯；細(xì)網(wǎng)格計(jì)算用時(shí)幾乎為粗網(wǎng)格的十倍，精度提高的優(yōu)勢(shì)也不明顯，因此中間網(wǎng)格（60×30）的計(jì)算方案是最有效的。與徐國(guó)徽文獻(xiàn)中的CFD計(jì)算結(jié)果（圖 4）及日本學(xué)者 Hinatsu（2001）模型試驗(yàn)結(jié)果（圖5）相比，變化趨勢(shì)和量值十分接近，表明本文數(shù)值計(jì)算方案可以接受。

3.5 體積模量減縮

圖5 模型試驗(yàn)測(cè)得P4點(diǎn)壓力時(shí)間歷程（60%H）Fig.5 Time history of pressure P4 measured in model tests(60%H)

為了進(jìn)一步開(kāi)展三維晃蕩問(wèn)題的數(shù)值模擬方法研究，必須提高計(jì)算效率縮短模擬時(shí)間。根據(jù)時(shí)間步長(zhǎng)公式，在保證計(jì)算精度歐拉網(wǎng)格尺寸一定的情況下，時(shí)間步長(zhǎng)與聲波在流體中傳播速度成反比，因此為降低波速可以減縮體積模量。水的真實(shí)體積模量K=2×103MPa，圖1的二維模型中采用1800個(gè)歐拉網(wǎng)格模擬10個(gè)周期的晃蕩運(yùn)動(dòng)耗時(shí)320分鐘（計(jì)算機(jī)內(nèi)存3.25 GB，Quad CPU 2.66 GHz），而當(dāng)體積模量減縮為K＝2 MPa時(shí)，同樣的模擬過(guò)程僅需23分鐘，大大提高了計(jì)算效率。圖6和圖7分別是真實(shí)體積模量和減縮體積模量情況下得到的壓力時(shí)間歷程，為了相互進(jìn)行比較，圖中對(duì)二者進(jìn)行了高頻濾波?？梢钥闯鰷p縮體積模量在相當(dāng)程度上減少了數(shù)值噪聲，并且壓力形式和時(shí)間特征并沒(méi)有受到影響。

圖6 P2點(diǎn)的壓力時(shí)間歷程（20%H）Fig.6 Time history of pressure P4(20%H)

圖7 P2點(diǎn)的壓力時(shí)間歷程（20%H，K減縮）Fig.7 Time history of pressure P4(20%H,reduced K)

4 三維剛性及三維彈性矩形液艙模型

三維剛性模型同二維的相似，但在流體的非主要運(yùn)動(dòng)方向增設(shè)單元，放開(kāi)流體僅在建模平面內(nèi)運(yùn)動(dòng)的限制，使其可在非主要運(yùn)動(dòng)方向也能運(yùn)動(dòng)。實(shí)際LNG船液艙的圍護(hù)系統(tǒng)是由兩層絕熱聚酯材料組成的。因此，我們考慮的三維彈性模型（圖9）在三維剛性（圖8）模型的基礎(chǔ)上，在流體域與剛性壁之間加了一層采用Lagrange體單元?jiǎng)澐值慕Y(jié)構(gòu)域模擬夾層泡沫，引用DMATEP卡片定義各向同性線彈性材料來(lái)模擬此層彈性結(jié)構(gòu)，參數(shù)設(shè)置依據(jù)真實(shí)的LNG船絕熱泡沫層物理量，體積模量K＝43.75 MPa，彈性模量E＝84 MPa。圖10和圖11是發(fā)生沖擊時(shí)液面圖。

圖8 三維剛性模型Fig.8 3D rigid model

圖9 三維彈性模型Fig.9 3D elastic model

圖10 14.3 s時(shí)刻自由液面及流體速度矢量圖（20%H）Fig.10 Free surface and fluid velocity vectors at 14.3 s

圖11 14.6 s時(shí)刻自由液面及流體速度矢量圖（20%H）Fig.11 Free surface and fluid velocity vectors at 14.6 s

5 壁面彈性對(duì)晃蕩沖擊壓力的影響分析

由于是三維模型，我們?nèi)∨c測(cè)點(diǎn)同樣水平高度的相鄰三單元做空間平均得到該處的壓力，通過(guò)剛性和彈性結(jié)構(gòu)模型處理計(jì)算得到的20%H裝載時(shí)液面附近P2點(diǎn)（圖12）和P3點(diǎn)（圖13）的壓力時(shí)間歷程，圖中我們可以觀察到壓力作用于彈性壁面的時(shí)間比剛性壁面的長(zhǎng)，這是由于彈性響應(yīng)形變的原因。

圖12 P2點(diǎn)的壓力時(shí)間歷程（20%H）Fig.12 Time history of pressure P2(20%H)

圖13 P3點(diǎn)的壓力時(shí)間歷程（20%H）Fig.13 Time history of pressure P3(20%H)

60%H裝載水平下泡沫材料的動(dòng)能時(shí)歷（圖14）反應(yīng)出彈性結(jié)構(gòu)在給定的簡(jiǎn)諧激勵(lì)下做簡(jiǎn)諧運(yùn)動(dòng)；同時(shí)通過(guò)泡沫材料的變形能時(shí)歷（圖15）我們可以看出彈性結(jié)構(gòu)在沖擊力作用下還發(fā)生了形變及反彈，開(kāi)始未發(fā)生沖擊時(shí)變形平緩，而到?jīng)_擊力作用時(shí)刻變形能值遠(yuǎn)高于其它時(shí)間；并且彈性結(jié)構(gòu)自身也存在反復(fù)振動(dòng)。

圖14 泡沫材料的動(dòng)能時(shí)間歷程（60%H）Fig.14 Time history of kinetic energy of foam(60%H)

圖15 泡沫材料的變形能時(shí)間歷程（60%H）Fig.15 Time history of deformation energy of foam(60%H)

局部結(jié)構(gòu)P6點(diǎn)處的大應(yīng)變發(fā)生時(shí)間與結(jié)構(gòu)大變形能（圖15）時(shí)刻基本一致，說(shuō)明整體結(jié)構(gòu)中該處的變形相對(duì)比較大，對(duì)變形能貢獻(xiàn)相對(duì)也最大；且應(yīng)變（圖16）主要沿X方向，即垂直于艙壁板格的方向也是晃蕩沖擊時(shí)刻流體相對(duì)艙壁板格運(yùn)動(dòng)方向，其它方向的應(yīng)變相對(duì)較小。對(duì)應(yīng)的等效應(yīng)力成分中（圖17），X方向的應(yīng)力占主要地位，其它方向應(yīng)力相對(duì)較小，剪切應(yīng)力與他們相比更小，放在同一幅圖中幾乎可以忽略。

圖16 P6點(diǎn)處單元應(yīng)變時(shí)間歷程（60%H）Fig.16 Time history of strain in element 55243 at P6(60%H)

圖17 P6點(diǎn)處單元應(yīng)力時(shí)間歷程（60%H）Fig.17 Time history of stress in element 55243 at P6(60%H)

圖18中是20%H裝載水平下P2和P3測(cè)點(diǎn)和60%H裝載水平下P4和P6測(cè)點(diǎn)的壓力統(tǒng)計(jì)值，液面下P2和P4點(diǎn)差異不大，而液面上的P3和P6點(diǎn)彈性模型計(jì)算值小于同工況下剛性模型計(jì)算值。

6 結(jié) 論

本文建立了基于顯式時(shí)間積分方法的有限元/有限體積法的流固耦合數(shù)值模擬技術(shù)，基于動(dòng)力學(xué)仿真軟件DYTRAN對(duì)液艙晃蕩數(shù)值模擬方法進(jìn)行了研究，并與CFD計(jì)算結(jié)果及試驗(yàn)結(jié)果進(jìn)行了對(duì)比分析。比較結(jié)果表明，本文提出的液艙晃蕩數(shù)值模擬方法是可行的。在此基礎(chǔ)上，進(jìn)一步對(duì)大型LNG船液艙圍護(hù)系統(tǒng)的彈性壁面帶來(lái)的流固耦合問(wèn)題開(kāi)展了初步研究，發(fā)現(xiàn)艙壁彈性對(duì)晃蕩壓力幅值和持續(xù)時(shí)間有一定影響，需要引起研究和設(shè)計(jì)者的重視。

圖18 剛性模型和彈性模型各測(cè)點(diǎn)的壓力均值Fig.18 Mean pressures at different measuring points in rigid and elastic models

[1]Bergheim P,Thiagarajan K P.The air-water sloshing problem:Parametric studies on excitation magnitude and frequency[C]//Proceedings of 27th International Conference on Offshore Mechanics and Arctic Engineering.Estoril,Portugal,2008.

[2]朱仁慶.液體晃蕩及其與結(jié)構(gòu)的相互作用[D].無(wú)錫:中國(guó)船舶科學(xué)研究中心,2002.

[3]Zhu Renqing,Fang Zhiyong,Wu Yousheng.Numerical simulation of viscous liquid sloshing coupled with elastic structures[J].Journal of Ship Mechanics,2006,10(3):61-70.

[4]娜日薩.VLCC液艙晃蕩仿真及其結(jié)構(gòu)強(qiáng)度評(píng)估方法研究[D].哈爾濱:哈爾濱工程大學(xué),2006.

[5]Rognebakke O,Faltinsen O M.Hydroelastic sloshing induced impact with entrapped air[J].Hydroelasticity in Marine Technology,2006,1:169-180.

[6]Wang B,Kim J W.Strength evaluation of LNG containment system considering fluid-structure interaction under sloshing impact pressure[C]//Proceedings of 26th International Conference on Offshore Mechanics and Arctic Engineering.San Diego,California,USA,2007.

[7]Hisashi lto.A direct assessment approach for structural strength evaluation of cargo containment system under sloshing inside LNGC tanks based on fluid structure interaction[C]//Proceedings of 27th International Conference on Offshore Mechanics and Arctic Engineering.Estoril,Portugal,2008.

[8]邢景棠,周 盛,崔爾杰.流固耦合力學(xué)概述[J].力學(xué)進(jìn)展,1997,27(1):19-38.

[9]徐國(guó)徽,顧學(xué)康.液艙晃蕩砰擊壓力的數(shù)值計(jì)算方法研究[C]//第二十二屆全國(guó)水動(dòng)力學(xué)研討會(huì)暨第九屆全國(guó)水動(dòng)力學(xué)學(xué)術(shù)會(huì)議.成都,2009.

[10]Hughes T J R,Liu W K,Thomas K Z.Lagrangian-Eulerian finite element formulation for incompressible viscous flows[J].Computer Methods in Applied Mechanics and Engineering,1981,29:329-349.

[11]Donea J,Giuliani S.An arbitrary Lagrangian-Eulerian finite element method for trainsient dynamic fluid-structure interaction[J].Computer Methods in Applied Mechanics and Engineering,1982,33:689-723.

[12]Huerta A,Liu W K.Viscous flow with large free surface motion[J].Computer Methods in Applied Mechanics and Engineering,1988,69(3):277-324.