辛建建, 方 田, 石伏龍

(1. 寧波大學(xué) 海運(yùn)學(xué)院， 浙江 寧波 315211； 2. 中國(guó)船舶科學(xué)研究中心， 江蘇 無(wú)錫 214082；3. 武漢理工大學(xué) 交通學(xué)院， 武漢 430063)

風(fēng)浪中航行的儲(chǔ)液船舶，如液化石油氣(LPG)、液化天然氣(LNG)船，在受到外部激勵(lì)時(shí)，液體會(huì)對(duì)容器壁產(chǎn)生強(qiáng)烈的砰擊作用，造成結(jié)構(gòu)的變形甚至破壞，對(duì)船舶穩(wěn)性和結(jié)構(gòu)強(qiáng)度造成嚴(yán)重威脅.液體發(fā)生晃蕩時(shí)會(huì)伴隨產(chǎn)生波面破碎、波對(duì)艙壁的高速砰擊、液滴形成及氣泡卷入等強(qiáng)非線(xiàn)性現(xiàn)象.而且，晃蕩引起的砰擊載荷會(huì)對(duì)艙壁產(chǎn)生強(qiáng)烈的局部沖擊，影響船舶的穩(wěn)定性，甚至破壞液艙的圍欄系統(tǒng)，導(dǎo)致沉船事故.因此，研究晃蕩現(xiàn)象和預(yù)報(bào)晃蕩砰擊載荷對(duì)液艙結(jié)構(gòu)設(shè)計(jì)和船舶的安全運(yùn)營(yíng)至關(guān)重要.

在過(guò)去，模型實(shí)驗(yàn)[1]和解析半解析方法[2]是船舶液艙晃蕩研究的常用手段.然而，模型實(shí)驗(yàn)費(fèi)用昂貴、條件苛刻(如低溫)，解析方法局限于簡(jiǎn)單幾何形狀液艙中的晃蕩問(wèn)題. 隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展，數(shù)值方法處理晃蕩問(wèn)題顯示出較大的優(yōu)越性，可以克服解析方法無(wú)法解決的困難，諸如液面大變形、搖蕩運(yùn)動(dòng)影響等.與費(fèi)用昂貴、時(shí)間較長(zhǎng)的模型實(shí)驗(yàn)相比，數(shù)值模擬方法的成本低廉、靈活方便.從自由表面跟蹤技術(shù)分，有流體體積分?jǐn)?shù)(VOF)法[3]、水平集(Level Set)法[4]和光滑粒子流體動(dòng)力學(xué)(SPH)法[5]等.根據(jù)計(jì)算網(wǎng)格可分為貼體網(wǎng)格法[6]、直角網(wǎng)格法[7]和無(wú)網(wǎng)格法[5].

在絕大多數(shù)的晃蕩流模擬中，對(duì)象為矩形液艙，邊界條件能直接施加，因?yàn)橛?jì)算網(wǎng)格與液艙邊界重合.然而在實(shí)際工程應(yīng)用中，液艙可能是任意形狀如菱形和橢圓形，另外液艙中也有可能存在內(nèi)部結(jié)構(gòu)如橫隔板.為了處理非矩形液艙或具有內(nèi)部結(jié)構(gòu)的情況，一個(gè)常用方法是貼體網(wǎng)格法[6].例如，Jiang等[8]基于 OpenFOAM模擬了菱形液艙的晃蕩問(wèn)題，Zhao等[9]基于重疊網(wǎng)格法模擬了部分充液 LNG 液艙的三維晃蕩流.然而，對(duì)于多自由度激勵(lì)、三維復(fù)雜形狀液艙和多個(gè)內(nèi)部結(jié)構(gòu)存在的情況，貼體網(wǎng)格方法面臨挑戰(zhàn).為了處理這些問(wèn)題，直角網(wǎng)格方法是一個(gè)具有吸引力的替代方法.在直角網(wǎng)格方法中，液艙邊界并不與背景網(wǎng)格線(xiàn)重合，通過(guò)引入力源項(xiàng)或重構(gòu)浸入邊界的插值模板以施加邊界條件.為了表示不規(guī)則邊界或內(nèi)部結(jié)構(gòu)物對(duì)晃蕩流的影響，Kim等[10]引入一個(gè)緩沖層以施加液艙斜邊界附近的自由表面邊界條件.Lee等[11]通過(guò)引入一個(gè)二次多項(xiàng)式內(nèi)插格式以施加液艙不規(guī)則邊界的無(wú)滑移邊界條件，根據(jù)不規(guī)則邊界與網(wǎng)格線(xiàn)之間的相對(duì)位置關(guān)系，內(nèi)插模板分為 13 種情況.

目前，三維復(fù)雜液艙劇烈晃蕩的數(shù)值預(yù)報(bào)仍面臨巨大的挑戰(zhàn)，文獻(xiàn)中也鮮見(jiàn)對(duì)菱形液艙共振晃蕩頻率的研究.為此，本文采用一個(gè)基于直角網(wǎng)格的三維多相流模型模擬了菱形液艙劇烈晃蕩問(wèn)題，并預(yù)報(bào)了菱形液艙的共振頻率，重點(diǎn)專(zhuān)注于艙壁沖擊壓力和全局晃蕩流運(yùn)動(dòng)，忽略局部物理現(xiàn)象湍流和可壓縮的影響.在規(guī)則網(wǎng)格上求解不可壓縮黏性N-S方程，在長(zhǎng)方體四角布置4個(gè)楔形體以表示菱形斜邊，采用徑向基函數(shù)虛擬網(wǎng)格法[12]，通過(guò)在固體區(qū)域內(nèi)部布置合適的虛擬網(wǎng)格，以考慮不規(guī)則邊界對(duì)流場(chǎng)的影響.

另外，采用一個(gè)三維梯度增量Level Set(GALS)兩相流模型[13]捕捉強(qiáng)非線(xiàn)性自由表面.基于所開(kāi)發(fā)的模型，模擬了菱形液艙的大幅晃蕩現(xiàn)象，并進(jìn)行了數(shù)值收斂性驗(yàn)證.進(jìn)一步研究了不同充液水深下壁面沖擊壓力隨激勵(lì)頻率的變化規(guī)律，確定了菱形液艙的共振頻率，以期為L(zhǎng)NG船舶的液艙結(jié)構(gòu)設(shè)計(jì)提供理論和技術(shù)指導(dǎo).

1 數(shù)學(xué)模型

不可壓縮N-S方程由描述任意流場(chǎng)的質(zhì)量和動(dòng)量守恒方程組成.在三維直角坐標(biāo)系下，非定常黏性不可壓縮N-S方程表達(dá)如下：

(1)

(2)

式中：速度向量u=[uvw]，u、v、w分別是直角網(wǎng)格x、y和z方向上的速度分量；t為時(shí)間；ρ為密度；p為壓力；υ為流體的運(yùn)動(dòng)黏性系數(shù).對(duì)于外部激勵(lì)力fa，其在多數(shù)自由表面流動(dòng)中為重力加速度g.然而，對(duì)于晃蕩問(wèn)題，外部力包括重力和晃蕩運(yùn)動(dòng)引起的慣性力以考慮液艙壁面運(yùn)動(dòng)對(duì)艙內(nèi)流場(chǎng)的影響，其表達(dá)式為[7]

(3)

式中：V和Ω分別為非慣性坐標(biāo)系下的平移和旋轉(zhuǎn)速度；r和R分別為網(wǎng)格節(jié)點(diǎn)坐標(biāo)和旋轉(zhuǎn)中心.式(3)右邊5項(xiàng)分別是非慣性坐標(biāo)系的重力加速度、平移、旋轉(zhuǎn)、離心及科氏加速度.

對(duì)于自由表面的捕捉，通過(guò)求解GALS方程實(shí)時(shí)更新距離函數(shù)φ和梯度ψ以隱式捕捉兩相界面.GALS方程[14]表示為

(4)

(5)

在兩相流模式計(jì)算中，根據(jù)自由表面的位置變化，定義在計(jì)算網(wǎng)格單元中心密度ρ及黏性系數(shù)υ的計(jì)算公式為

ρ=Hε(φ)ρW+(1-Hε(φ))ρA

(6)

υ=Hε(φ)υW+(1-Hε(φ))υA

(7)

式中：Hε(φ)為Heaviside函數(shù)以區(qū)分流體和空氣界面的流體屬性；ρ及υ的下標(biāo)W和A分別表示水和空氣.為了減緩因密度(或黏性)介質(zhì)的跳躍對(duì)流動(dòng)的影響，使得過(guò)渡層內(nèi)物性參數(shù)光滑變化，抑制數(shù)值振蕩.

2 數(shù)值求解方法

2.1 N-S方程求解器

在本文研究中，基于自主開(kāi)發(fā)的時(shí)間半隱式有限差分法，在交錯(cuò)直角網(wǎng)格上求解非定常不可壓縮N-S方程.采用分步法解耦速度和壓力，TVD-RK2(Total Variation Diminishing Second-order Runge-Kutta)格式進(jìn)行時(shí)間推進(jìn)，對(duì)流項(xiàng)采用時(shí)間顯式處理，黏性項(xiàng)采用Crank-Nicolson格式半隱式處理.對(duì)于空間離散，對(duì)流項(xiàng)采用高階TVD-MUSCL(Monotonic Upstream-centered Scheme for Conservation Laws)格式離散，黏性項(xiàng)以中心差分格式離散.另外，對(duì)于壓力泊松方程離散形成的大型稀疏線(xiàn)性方程組，采用預(yù)處理的BICGSTAB(Bi-conjugate Gradients Stabilize)方法結(jié)合一種對(duì)角存儲(chǔ)格式求解.更多的細(xì)節(jié)請(qǐng)參考文獻(xiàn)[12].

2.2 GALS方法數(shù)值求解

GALS方程式(4)、(5)以特征形式表示為

(8)

求解過(guò)程分為兩步：① 采用Shu-Osher RK3方法從t+Δt到t進(jìn)行向后時(shí)間積分求解特征曲線(xiàn)方程?x(τ)/?τ=u(x(τ),τ)；② 然后沿著特征曲線(xiàn)采用顯式梯度格式進(jìn)行向前積分求解式(8)以更新φ和ψ.值得注意的是在任意位置的變量φ和ψ通過(guò)Hermite立方多項(xiàng)式插值得到以獲得三階精度，另外，通過(guò)對(duì)Hermite立方多項(xiàng)式格式進(jìn)行微分以解析方式計(jì)算φ和ψ的一階和二階導(dǎo)數(shù).詳細(xì)求解過(guò)程參考文獻(xiàn)[13].

2.3 不規(guī)則邊界處理

浸入邊界法的網(wǎng)格生成過(guò)程簡(jiǎn)單，但由于浸入邊界本身并未包含在計(jì)算網(wǎng)格中，精確表示背景網(wǎng)格中浸入邊界的存在是個(gè)難點(diǎn).本文采用基于徑向基函數(shù)(RBF)的隱式等值面表示方法[12]，光滑基函數(shù)擬合離散的物面控制點(diǎn)以構(gòu)造任意復(fù)雜物體的等值曲面，并且根據(jù)等值面距離函數(shù)的正負(fù)號(hào)識(shí)別場(chǎng)點(diǎn)屬性.下一步就可重構(gòu)虛擬網(wǎng)格的流體變量值以施加邊界條件.通過(guò)在浸入物體內(nèi)部布置有限數(shù)量的虛擬網(wǎng)格以表示靜止或運(yùn)動(dòng)邊界的存在，以適當(dāng)?shù)牟逯的K及技術(shù)(如雙線(xiàn)性插值方法)，施加浸入動(dòng)邊界條件.

虛擬網(wǎng)格的插值重構(gòu)過(guò)程可細(xì)分為3步：標(biāo)記鏡像點(diǎn)、插值重構(gòu)和滿(mǎn)足鏡像條件.首先，對(duì)于每個(gè)虛擬網(wǎng)格(GC)，沿著物體表面的外法向從虛擬網(wǎng)格單元中心延伸到流體域以確定唯一的一個(gè)鏡像點(diǎn)(MP)，如圖1所示.然后采用三線(xiàn)性插值格式計(jì)算每個(gè)鏡像點(diǎn)的流體變量.然而，存在一個(gè)或更多的模板點(diǎn)是虛擬網(wǎng)格本身的情況，為了避免繁瑣的內(nèi)插模板分類(lèi)過(guò)程，首先對(duì)封閉鏡像點(diǎn)(MP)的8個(gè)鄰近網(wǎng)格節(jié)點(diǎn)進(jìn)行判斷.如果該網(wǎng)格節(jié)點(diǎn)為流體網(wǎng)格點(diǎn)(FC)，則將它確定為模板點(diǎn)，1、2、3、4即為二維插值格式中識(shí)別出的模板點(diǎn)，如圖1(a)所示.如果為虛擬網(wǎng)格點(diǎn)(GC)，則相應(yīng)的邊界點(diǎn)(BP)作為模板點(diǎn)，而不是虛擬網(wǎng)格點(diǎn)本身，如圖1(b)所示.如果更多的鄰近網(wǎng)格點(diǎn)識(shí)別為虛擬網(wǎng)格點(diǎn)，同樣，相應(yīng)的邊界點(diǎn)BP作為模板點(diǎn)(見(jiàn)圖1(c)).在求得鏡像點(diǎn)的流體變量后，虛擬網(wǎng)格點(diǎn)的流體變量可通過(guò)滿(mǎn)足物面的鏡像條件直接得到.當(dāng)通過(guò)標(biāo)準(zhǔn)的七點(diǎn)中心差分格式在整個(gè)計(jì)算域求解N-S方程后，物面邊界條件就會(huì)得到隱式施加.更多的求解細(xì)節(jié)見(jiàn)文獻(xiàn)[12].

圖1 虛擬網(wǎng)格法的二維插值格式Fig.1 Two-dimensional interpolation format of virtual grid method

3 三維菱形液艙晃蕩分析

為研究復(fù)雜液艙系統(tǒng)下的晃蕩特性，數(shù)值模擬了不同水深和激勵(lì)頻率下橫搖激勵(lì)的三維菱形液艙晃蕩問(wèn)題.根據(jù)Mikelis等[15]對(duì)不同充液水平下棱形液艙的晃蕩實(shí)驗(yàn)研究，液艙模型的長(zhǎng)、寬、高分別為L(zhǎng)=0.741 m、W=0.399 m、H=0.405 m，其幾何模型如圖2所示.液艙遭受強(qiáng)迫橫搖激勵(lì)，橫搖中心Y1、Y2在x軸中心位置距離液艙底部0.194 m.液艙橫搖激勵(lì)的運(yùn)動(dòng)描述為θy=θ0sinωt.θy為繞橫搖中心Y1、Y2的橫搖角；橫搖幅值θ0=0.1 rad；ω為激勵(lì)頻率，ω=2π/T，T為橫搖周期.在右壁面寬度方向中心位置設(shè)置7個(gè)壓力傳感器R1～R7，距離液艙底部的距離分別為0.075，0.113 5，0.205，0.296 5，0.341，0.390，0.405 m，R7距離右壁面0.02 m.在寬度中心距離左邊界 0.135 6 m位置設(shè)置浪高儀WR.在液艙所有邊界設(shè)置速度無(wú)滑移邊界條件，在上壁面設(shè)置壓力為0，其余邊界設(shè)置Neumann邊界條件.由于采用規(guī)則的矩形計(jì)算區(qū)域，為了考慮菱形液艙的斜邊界對(duì)流場(chǎng)的影響，在液艙橫剖面的四角布置4個(gè)楔形體，采用虛擬網(wǎng)格法施加楔形體表面的物面邊界條件.

圖2 菱形液艙幾何模型和傳感器布置(m)Fig.2 Geometrical model of prismatic tank and placement of sensors (m)

3.1 數(shù)值驗(yàn)證

首先模擬了充液比h/H=0.9 (h為水深)和周期T=0.97 s的工況以驗(yàn)證本文多相流方法的精度和可靠性.分別采用80×50×60、100×60×80和120×80×100的3套網(wǎng)格離散計(jì)算區(qū)域以研究網(wǎng)格收斂性，根據(jù)CFL(Courant-Friedreichs-Lewy)條件決定變時(shí)間步長(zhǎng)，其中時(shí)間松弛系數(shù)設(shè)置為ωCFL=0.3.本文方法3套網(wǎng)格下計(jì)算的在R3點(diǎn)的壓力和WR位置的波浪爬高與實(shí)驗(yàn)結(jié)果進(jìn)行了比較，如圖3所示.圖中η為波浪爬高.在3套網(wǎng)格下的壓力和波浪幅值與實(shí)驗(yàn)數(shù)據(jù)吻合較好，盡管觀察到峰值壓力有略微差別.然而，隨著網(wǎng)格細(xì)化，峰值壓力逐漸趨近實(shí)驗(yàn)數(shù)據(jù)，如圖3(a)所示，驗(yàn)證了本文方法良好的網(wǎng)格收斂性.另外，波峰在η=0.4 m的位置發(fā)生截?cái)?，因?yàn)椴ɡ伺郎搅艘号擁敳?為了研究時(shí)間步長(zhǎng)的收斂性，選擇3個(gè)不同的時(shí)間松弛系數(shù)ωCFL=0.1、ωCFL=0.3和ωCFL=0.6，分別在x、y和z方向選用100×60×80的計(jì)算網(wǎng)格.圖4展示了3個(gè)時(shí)間松弛系數(shù)下計(jì)算的壓力和波浪爬升與實(shí)驗(yàn)數(shù)據(jù)的比較結(jié)果.總體而言，在3個(gè)時(shí)間松弛系數(shù)下本文計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)吻合.隨著時(shí)間步長(zhǎng)的縮小，壓力峰值會(huì)更加尖銳，更趨向于實(shí)驗(yàn)結(jié)果，ωCFL=0.1和ωCFL=0.6時(shí)的壓力峰值相比實(shí)驗(yàn)數(shù)據(jù)分別略小4.8%和7.1%，偏差在合理范圍內(nèi)，再一次驗(yàn)證了本文方法的精度和不同時(shí)間步長(zhǎng)下的可靠性.為了在求解質(zhì)量和計(jì)算成本間取得平衡， 在其他工況計(jì)算中時(shí)間松弛系數(shù)選為ωCFL=0.3，網(wǎng)格設(shè)置為100×60×80.

圖3 h/H=0.9時(shí)三套網(wǎng)格下本文方法和實(shí)驗(yàn)得到的壓力和波浪爬高時(shí)間歷程Fig.3 Time history of pressure and wave elevation by present method on three grids and experiment at h/H=0.9

圖4 h/H=0.9時(shí)三個(gè)時(shí)間松弛系數(shù)下本文方法和實(shí)驗(yàn)得到的壓力和波浪爬高時(shí)間歷程Fig.4 Time history of pressure and wave elevation by present method on three time coefficients and experimental data at h/H=0.9

分別模擬了h/H=0.3、T=1.517 s和h/H=0.61、T=1.112 s的兩個(gè)工況以進(jìn)一步驗(yàn)證本文方法的可靠性.圖5、6分別為兩個(gè)工況下波高和3個(gè)測(cè)量點(diǎn)壓力與Mikelis等實(shí)驗(yàn)數(shù)據(jù)的比較結(jié)果.當(dāng)h/H=0.3時(shí)，在垂向壁面R2點(diǎn)的壓力脈沖值略高于實(shí)驗(yàn)結(jié)果，自由表面曲線(xiàn)在波谷與實(shí)驗(yàn)值吻合，但波峰略高于實(shí)驗(yàn)值(見(jiàn)圖5).當(dāng)水深繼續(xù)增加到h/H=0.61時(shí)，壓力變化趨勢(shì)和形態(tài)與實(shí)驗(yàn)結(jié)果吻合，自由表面曲線(xiàn)的波高略大于實(shí)驗(yàn)值(見(jiàn)圖6).本文結(jié)果與實(shí)驗(yàn)值的偏差可能是由于湍流或空氣可壓縮的影響所致，考慮到三維晃蕩的復(fù)雜性、非線(xiàn)性和瞬態(tài)性，偏差在合理范圍內(nèi).從圖5可看出，晃蕩波周期性地砰擊壁面產(chǎn)生脈沖形壓力曲線(xiàn).同時(shí)，在每個(gè)周期內(nèi)出現(xiàn)連續(xù)的雙波峰現(xiàn)象，其由強(qiáng)烈的非線(xiàn)性波浪砰擊產(chǎn)生.當(dāng)液艙運(yùn)動(dòng)到一側(cè)時(shí)，流體被高速運(yùn)動(dòng)液艙驅(qū)動(dòng)到同一側(cè)，快速運(yùn)動(dòng)的水體砰擊垂向壁面，導(dǎo)致第1個(gè)壓力峰值產(chǎn)生，盡管液艙達(dá)到最大橫搖角，但水體在慣性作用下繼續(xù)向前運(yùn)動(dòng)，導(dǎo)致第2個(gè)壓力峰值產(chǎn)生.當(dāng)水深增加到h/H=0.61時(shí)，由于此時(shí)晃蕩流的非線(xiàn)性程度降低，壓力曲線(xiàn)呈現(xiàn)周期性的單壓力峰值現(xiàn)象，自由表面曲線(xiàn)呈現(xiàn)較為規(guī)則的正弦曲線(xiàn)變化趨勢(shì)，R2點(diǎn)的壓力呈現(xiàn)周期性的簡(jiǎn)諧振蕩特性，均衡壓力在 1 250 Pa左右，該點(diǎn)在靜水面以下，砰擊壓力由靜水壓力和動(dòng)壓共同主導(dǎo).

圖6 h/H=0.61時(shí)的的壓力和自由表面時(shí)間歷程Fig.6 Time history of pressure and free surface at h/H=0.61

圖7、8分別為兩個(gè)工況下在兩個(gè)瞬時(shí)時(shí)刻的三維波面圖，其中4個(gè)藍(lán)色楔形體表示菱形液艙的幾何外形，用來(lái)考慮不規(guī)則邊界與流場(chǎng)的相互作用.h/H=0.3時(shí)，整個(gè)水體隨著液艙搖蕩向一側(cè)產(chǎn)生劇烈運(yùn)動(dòng)，水體上層流體速度大于底部流體，形成大的涌浪，如圖7(a)所示.當(dāng)液艙達(dá)到最大橫搖角時(shí)，涌浪砰擊壁面，隨后，液艙向反方向運(yùn)動(dòng)，驅(qū)動(dòng)水體向反方向運(yùn)動(dòng)，并形成相同的涌浪，如圖7(b)所示.h/H=0.61時(shí)，水體跟隨著液艙運(yùn)動(dòng)并一直爬升到液艙頂部斜邊，生成陡峭的、與上斜邊貼合的非線(xiàn)性波浪.當(dāng)波浪爬升到斜邊最大位置時(shí)，其在重力作用下下落，形成波浪翻卷，同時(shí)，后方的波浪繼續(xù)向前運(yùn)動(dòng)，形成皺褶形的波浪紋，如圖8所示.

圖7 h/H=0.3時(shí)兩時(shí)刻的瞬時(shí)自由表面Fig.7 Instantaneous free surface at two moments at h/H=0.3

圖8 h/H=0.61時(shí)兩時(shí)刻的瞬時(shí)自由表面Fig.8 Instantaneous free surface at two moments at h/H=0.61

3.2 共振頻率預(yù)報(bào)

在載液船舶液艙設(shè)計(jì)中，液艙共振頻率(共振頻率≈自然頻率)預(yù)報(bào)對(duì)于避免危險(xiǎn)工況的發(fā)生、保障航行安全具有重要的意義.對(duì)于矩形液艙，液艙自然頻率可解析計(jì)算為

(9)

m,n=0, 1, 2, …

圖9 4個(gè)充液水深下艙壁壓力幅值與激勵(lì)頻率之間的關(guān)系Fig.9 Pressure amplitude at tank walls versus excitation frequency at four filling water levels

然而，對(duì)于不規(guī)則形狀液艙如菱形液艙，其自然頻率無(wú)法解析計(jì)算得到.鑒于此，本文數(shù)值研究了不同水深下液艙激勵(lì)頻率與艙壁局部位置壓力幅值之間的關(guān)系，進(jìn)一步確定了不同水深下菱形液艙的共振頻率.圖9展示了4個(gè)充液水深下艙壁壓力測(cè)量點(diǎn)上的壓力幅值隨液艙激勵(lì)頻率的變化趨勢(shì)，壓力幅值取穩(wěn)定狀態(tài)時(shí)5個(gè)周期壓力峰值的平均值，以避免初始擾動(dòng)和瞬時(shí)波動(dòng)的影響.可以看出，隨著激勵(lì)頻率的增加，艙壁各點(diǎn)處的壓力幅值迅速增加到最大值，當(dāng)激勵(lì)頻率繼續(xù)增加時(shí)，壓力幅值迅速減小，而且壓力幅值減小的速度大于增加的速度.壓力幅值最大時(shí)對(duì)應(yīng)的激勵(lì)頻率即為共振頻率，也就是自然頻率.由此可得到不同水深下菱形液艙的自然頻率，如表1所示.為了提供參考，計(jì)算了相同主尺度下矩形液艙的最低階自然頻率，此時(shí)m=1,n=0.由表1可看出，相同水深下菱形液艙的最低階自然頻率小于相同主尺度下矩形液艙的最低階自然頻率，原因可能是液艙與晃蕩流體強(qiáng)烈的非線(xiàn)性相互作用.應(yīng)該指出的是，本文預(yù)報(bào)激勵(lì)頻率的精度為小數(shù)點(diǎn)后一位，因?yàn)椴捎玫淖钚〖?lì)頻率間隔是0.1 rad/s，為了提高預(yù)報(bào)精度，需要取更小的激勵(lì)頻率步長(zhǎng).

表1 4個(gè)充液水深下矩形液艙和菱形液艙的自然頻率Tab.1 Natural frequency of rectangular and prismatic tanks at four filling water levels

隨著水深的增加，自然頻率增加，與式(9)一致.而且，隨著水深的增加，晃蕩波浪產(chǎn)生的沖擊壓力幅值逐漸減小，從h/H=0.46時(shí)P1點(diǎn)的 2 300 Pa到h/H=0.46時(shí)P1點(diǎn)的 1 200 Pa.從3.1節(jié)也可看出，隨著充液比的增加，壁面沖擊壓力的非線(xiàn)性變化特征逐漸減弱，從h/H=0.3時(shí)的雙壓力峰值現(xiàn)象過(guò)渡到h/H=0.61時(shí)的單壓力峰值現(xiàn)象，最后在h/H=0.9時(shí)呈現(xiàn)規(guī)則的正弦曲線(xiàn)變化特征.因?yàn)殡S著充液比的增加，晃蕩液體的運(yùn)動(dòng)幅值會(huì)更多受到周?chē)后w黏性阻尼作用的約束，而且在高充液比如h/H=0.75或0.9時(shí)，從圖9(c)、9(d)中點(diǎn)P5的壓力幅值可看出，晃蕩波浪砰擊到液艙頂部，頂部砰擊的阻尼作用也會(huì)對(duì)波浪的非線(xiàn)性產(chǎn)生抑制作用.另外，在充液比h/H=0.46或0.61時(shí)，底部斜邊點(diǎn)P1和壁面下部P2點(diǎn)遭受巨大的沖擊壓力，超過(guò) 2 000 Pa，壁面中上部的沖擊壓力也較大，接近 1 500 Pa，如圖9(a)、9(b)，說(shuō)明此水深下整個(gè)垂向壁面都遭受較大的砰擊載荷，液艙壁的整體強(qiáng)度需要加強(qiáng).在充液比增加到h/H=0.75或0.90時(shí)，艙壁的整體沖擊壓力顯著降低，但艙壁上斜邊和頂部會(huì)遭受顯著的沖擊壓力.由此可知，載液船舶在海上航行時(shí)，理想的充液水深是接近滿(mǎn)載時(shí).當(dāng)充液水深近似為液艙高度一半時(shí)，液艙處于非常不利的工作環(huán)境.

4 結(jié)論

采用基于直角網(wǎng)格法的三維多相流模型模擬了橫搖激勵(lì)下菱形液艙大幅晃蕩問(wèn)題，研究了不同水深下艙壁局部位置壓力幅值與激勵(lì)頻率之間的關(guān)系，得出菱形液艙的共振頻率，主要結(jié)論有：

(1) 高充液比下晃蕩研究表明本文計(jì)算模型具有良好的網(wǎng)格和時(shí)間收斂性，不同充液比下的沖擊壓力和波浪爬高也與實(shí)驗(yàn)數(shù)據(jù)吻合，驗(yàn)證了本文方法的精度.

(2) 隨著充液比增加，壁面沖擊壓力的非線(xiàn)性變化特征逐漸減弱，從充液比h/H=0.3時(shí)的雙壓力峰值現(xiàn)象過(guò)渡到h/H=0.61時(shí)的單壓力峰值現(xiàn)象，最后在h/H=0.9時(shí)呈現(xiàn)規(guī)則的正弦曲線(xiàn)變化特征.

(3) 同一充液比下菱形液艙自然頻率低于同主尺度矩形液艙的自然頻率，隨著充液比增加，最大的壓力幅值逐漸減小，因?yàn)橛懈嗟酿ば宰枘嶙饔?，而且高充液比時(shí)，也存在頂部砰擊的阻尼作用.

(4) 當(dāng)充液水深近似為液艙高度一半時(shí)，整個(gè)液艙垂向艙壁都遭受較大的砰擊載荷，液艙處于非常不利的工作環(huán)境.