沖擊載荷作用下LNG燃料罐防波板參數(shù)多目標(biāo)優(yōu)化

馬寒陽(yáng)，董金善

（1.南京工業(yè)大學(xué)機(jī)械與動(dòng)力工程學(xué)院，南京 211816；2.華東理工大學(xué)機(jī)械與動(dòng)力工程學(xué)院，上海 200237）

0 引 言

船用LNG燃料罐作為一種移動(dòng)式低溫壓力容器，在運(yùn)行的過(guò)程中不僅要承受內(nèi)外壓載荷，還會(huì)受到由于慣性力或碰撞引起的液體沖擊載荷。沖擊載荷引發(fā)的液體晃動(dòng)現(xiàn)象會(huì)帶來(lái)不利的影響：一方面，對(duì)燃料罐結(jié)構(gòu)產(chǎn)生沖擊載荷；另一方面，液體的大幅晃動(dòng)將會(huì)改變載荷分布，產(chǎn)生附加的力和力矩，從而降低船的行駛穩(wěn)定性。因此，開(kāi)展船用LNG儲(chǔ)罐內(nèi)液體晃動(dòng)的數(shù)值方法研究及防波板的設(shè)計(jì)具有重要的學(xué)術(shù)價(jià)值和工程指導(dǎo)意義。目前，晃蕩問(wèn)題主要從理論、實(shí)驗(yàn)和模擬三個(gè)方向進(jìn)行研究。Li[1]研究了直接時(shí)域理論和間接時(shí)域理論對(duì)船舶運(yùn)動(dòng)的時(shí)延函數(shù)的影響，結(jié)果表明間接時(shí)域理論在處理無(wú)航速和較低航速工況下顯得更加高效實(shí)用；楊志勛[2]采用二維八邊形艙晃蕩模型進(jìn)行了實(shí)驗(yàn)，利用廣義極值分布、三參Weibull分布和兩參Weibull分布預(yù)測(cè)了晃蕩荷載峰值；Elahi[3]使用VOF方法研究了直線加速度和角加速度作用下二維燃料罐的晃蕩問(wèn)題，并與理論計(jì)算和SPH方法進(jìn)行了對(duì)比，證實(shí)了VOF方法的準(zhǔn)確性。

防波板的安裝能夠有效地減少燃料罐封頭承受的壓力，提高內(nèi)外罐間徑向支撐結(jié)構(gòu)環(huán)氧玻璃鋼的支撐能力。國(guó)內(nèi)外學(xué)者對(duì)防波板進(jìn)行了大量的研究：Kim[4]利用神經(jīng)網(wǎng)絡(luò)和遺傳算法對(duì)矩形容器在平移載荷作用下的防波板進(jìn)行分析，得到了防波板的最佳安裝高度；王鶴鵬[5]研究了防波板的數(shù)量、位置和幾何形狀對(duì)晃蕩沖擊力的影響；Ma[6]采用VOF模型和大渦模型（LES）模擬研究了不同類型垂直防波板對(duì)液體瞬態(tài)晃動(dòng)的有效性。結(jié)果表明兩個(gè)垂直擋板的適當(dāng)分布可以將峰值沖擊壓力降低三倍左右。不過(guò)，國(guó)內(nèi)外對(duì)防波板的研究主要集中在防波板對(duì)液體沖擊力的減少上，很少有學(xué)者針對(duì)防波板自身的強(qiáng)度進(jìn)行分析。目前，僅有Shimanovsky[7]和Kang[8]等對(duì)防波板的強(qiáng)度進(jìn)行了分析。且國(guó)產(chǎn)防波板易產(chǎn)生裂紋或松脫，據(jù)統(tǒng)計(jì)約占總數(shù)的90%以上[9]。因此，對(duì)防波板的防波效果和結(jié)構(gòu)強(qiáng)度兩個(gè)目標(biāo)進(jìn)行優(yōu)化分析，得到最佳的防波板形狀和位置具有重要的意義。

1 理論模型

1.1 控制方程

LNG船用燃料罐中流體部分的控制方程是Navier-Stokes方程和連續(xù)性方程。由于湍流過(guò)程非常復(fù)雜，根據(jù)統(tǒng)計(jì)學(xué)特征規(guī)律，采用流動(dòng)的時(shí)均特性信息（例如平均速度、平均壓力等）表示方程中任一瞬時(shí)的物理量。將Reynold時(shí)均方法用于湍流模型中，則Navier-Stokes方程形式為

式中：δij為Kronecker符號(hào)；μt為湍動(dòng)粘性系數(shù)，其取決于流體流動(dòng)狀態(tài)，而不是流體的物性參數(shù)。可建立μt與湍流動(dòng)能k和湍流動(dòng)能耗散率ε的關(guān)系式

式中，cμ是通過(guò)大量數(shù)據(jù)統(tǒng)計(jì)得到的適用于大部分湍流的常數(shù)值，一般取0.09。

1.2 VOF模型

1.3 Fluent環(huán)境設(shè)置

Fluent采用的有限體積法離散方程適用于結(jié)構(gòu)化網(wǎng)格和非結(jié)構(gòu)化網(wǎng)格。同時(shí)動(dòng)量修正方程和壓力修正方程分別使用二階迎風(fēng)格式和體積加權(quán)法。對(duì)于求解壓力速度耦合關(guān)系，則采用POSI算法[10]。對(duì)于沖擊載荷，通過(guò)UDF（User Defined Function）編寫動(dòng)量源接口程序來(lái)施加。因此，計(jì)算域中任一有限元網(wǎng)格單元在X、Y上的表達(dá)式為[11]

式中：SX和SY分別為計(jì)算域中任意有限單元在x和y上的動(dòng)量源項(xiàng)，單位為N/m3；αx和αy分別為計(jì)算域中任意有限單元在x和y上的加速度（本文運(yùn)動(dòng)方向加速度取19.62[12]），單位為m/s2；ρ1為氣相（空氣）的密度（1.225 kg/m3）；ρ2為液相（LNG）的密度（500 kg/m3）；Qx、Qy為計(jì)算域中任意有限單元?dú)庀?、液相的體積分?jǐn)?shù)。

2 防波板有限元計(jì)算

2.1 模型介紹

本文以船用低溫C型LNG燃料罐內(nèi)罐為研究對(duì)象，筒體內(nèi)徑為1 400 mm，兩端為標(biāo)準(zhǔn)橢圓形封頭，筒體長(zhǎng)度為3 100 mm，有效容積為5 m3。罐內(nèi)設(shè)置錐形防波板，防波板結(jié)構(gòu)參數(shù)如圖1所示。因模型對(duì)稱，為提高計(jì)算效率，取1/2模型進(jìn)行分析。防波板及筒體選用S30408，材料性能如表1所示。

圖1 防波板結(jié)構(gòu)參數(shù)Fig.1 Baffle structure parameters

表1 材料性能參數(shù)Tab.1 Material property parameters

2.2 流體有效性驗(yàn)證

本文借鑒Modaressi-Tehrani[13]論文中使用的模型和參數(shù)，用來(lái)驗(yàn)證流體算法的準(zhǔn)確性。Modaressi-Tehrani采用QS方法計(jì)算瞬態(tài)液體在縱向加速度為0.3g時(shí)對(duì)罐體壁面的沖擊力。本文通過(guò)Fluent有限差分法迭代計(jì)算充液率為0.4和0.6兩種情況下在一個(gè)振蕩周期內(nèi)的平均值。對(duì)比情況如圖2所示。

圖2 Fluent計(jì)算值和QS計(jì)算值對(duì)比Fig.2 Calculated value comparison of Fluent and quasi-static

可以看出，數(shù)值計(jì)算結(jié)果和理論計(jì)算結(jié)果吻合良好。但是，QS方法只能計(jì)算平均沖擊力的大小，而不能得到震蕩過(guò)程中的最大沖擊力。本文將采用最大沖擊力來(lái)計(jì)算晃蕩對(duì)防波板產(chǎn)生的影響。

2.3 不同充液率對(duì)罐體沖擊力的影響

充裝率是影響罐體壁面受力的主要因素之一，但燃料罐罐內(nèi)液體充裝率會(huì)經(jīng)歷一個(gè)從高到低變化的過(guò)程。因此，應(yīng)該以產(chǎn)生最大沖擊力時(shí)的充裝率作為研究對(duì)象。根據(jù)文獻(xiàn)[14-15]可知，在充裝率較小時(shí)不會(huì)產(chǎn)生最大沖擊力，因此，本文取0.5、0.55、0.6、0.65、0.7、0.75、0.8、0.85和0.9九種工況進(jìn)行研究，如圖3所示。

圖3 沖擊力與充液率的關(guān)系Fig.3 Correlation between impact force and filling ratio

從圖中可知，液體晃動(dòng)產(chǎn)生的沖擊力先隨著充裝率增大而增大，充裝率大于0.85后，晃動(dòng)產(chǎn)生的沖擊力隨充裝率的增加而減小，這與文獻(xiàn)[14]中研究得到的結(jié)果相同。發(fā)生這種現(xiàn)象的原因主要是因?yàn)楫?dāng)充裝率較大時(shí)，液體的質(zhì)量也比較大，因此，液體的晃動(dòng)幅度就會(huì)比較小，所產(chǎn)生的沖擊力相應(yīng)減小。

2.4 數(shù)值模擬結(jié)果與分析

根據(jù)JB4732《鋼制壓力容器—分析設(shè)計(jì)標(biāo)準(zhǔn)》[16]和《天然氣燃料動(dòng)力船舶規(guī)范》[12]，對(duì)于防波板的應(yīng)力評(píng)定是基于第三強(qiáng)度理論，并對(duì)其進(jìn)行線性化處理。線性化處理是選擇危險(xiǎn)截面并把各應(yīng)力分量沿一條應(yīng)力處理線進(jìn)行應(yīng)力計(jì)算和分類。根據(jù)計(jì)算結(jié)果，得到的最大應(yīng)力點(diǎn)位于防波板與筒體上角鋼連接處，從最大應(yīng)力點(diǎn)沿網(wǎng)格厚度方向?qū)?yīng)力進(jìn)行線性化。為保證線性化結(jié)果的可靠性，在防波板厚度方向至少有三層以上網(wǎng)格密度。應(yīng)力云圖如圖4所示。

圖4 LNG燃料罐應(yīng)力云圖Fig.4 Stress contour of LNG fuel tank

根據(jù)線性化結(jié)果，薄膜應(yīng)力為155.02 MPa≤1.5Sm，彎曲應(yīng)力強(qiáng)度為236.24 MPa≤3Sm。由此可知，C型LNG燃料罐防波板滿足強(qiáng)度要求，并以此為目標(biāo)進(jìn)行優(yōu)化。

3 防波板優(yōu)化設(shè)計(jì)

3.1 正交試驗(yàn)設(shè)計(jì)

防波板的結(jié)構(gòu)復(fù)雜且可變參數(shù)較多。參考之前學(xué)者進(jìn)行的大量研究，選取防波板距封頭安裝距離、防波板底邊距中心線距離、防波板頂邊距中心線距離和防波板形狀四個(gè)參數(shù)進(jìn)行研究，每個(gè)參數(shù)在限定條件內(nèi)取4個(gè)水平。利用有限元軟件的流固耦合分別獲取封頭處的最大沖擊力（y1）、防波板薄膜應(yīng)力（y2）和防波板彎曲應(yīng)力（y3）計(jì)算結(jié)果。如果把每一種工況都進(jìn)行計(jì)算則需要進(jìn)行256次計(jì)算。為了節(jié)省計(jì)算成本，采用正交設(shè)計(jì)法設(shè)計(jì)試驗(yàn)。設(shè)計(jì)正交表時(shí)留有一個(gè)空列，以便進(jìn)行殘差分析。因此本文采用L16（45）正交表。

列出本文采用的試驗(yàn)方案和仿真結(jié)果，如表2所示。其中，A代表防波板距封頭安裝距離，B代表防波板底邊距中心線距離，主要是為了使防波板上部弓形面積小于罐體橫斷面積的20%；C代表防波板頂邊距中心線距離，主要是使防波板的有效面積大于罐體橫斷面積的40%；D代表防波板形狀為平板形、錐形、橢圓形和波紋板形，平板型防波板即為厚度為5 mm的平板；錐形防波板的半頂角α=70°；橢圓形防波板是由長(zhǎng)短軸之比為四比一的橢圓構(gòu)成的；波紋板形防波板的中間區(qū)域是由振幅為25 mm，周期為400 mm的正弦函數(shù)構(gòu)成，并隨著距離筒體越近振幅逐漸減少直至變成平板，為方便之后進(jìn)行參數(shù)化研究，令它們的形狀系數(shù)[17]值為1、2、3、4。

表2 試驗(yàn)方案和仿真結(jié)果Tab.2 Experimental schemes and simulation results

3.2 方差分析

方差分析可以了解各因素對(duì)模擬結(jié)果的影響，還可以把由試驗(yàn)條件引起的數(shù)據(jù)波動(dòng)和由試驗(yàn)誤差引起的數(shù)據(jù)波動(dòng)區(qū)分開(kāi)來(lái)。對(duì)于Ln(mk)正交表，把實(shí)驗(yàn)數(shù)據(jù)按表3填寫。

表3 正交表Tab.3 Orthogonal table

方差分析的關(guān)鍵是偏差平方和的分解。

總偏差平方和與總自由度為

各列偏差平方和與自由度為

若將因素A安排在正交表的第j列(j=1,2,...,k)上，則有SSA=SSj，這里的r=n/m。

誤差平方和與自由度為

于是得方差分析表4。

表4 方差分析Tab.4 Analysis of variance

表5為各因素對(duì)封頭沖擊力和防波板應(yīng)力值的方差分析結(jié)果。

表5 y1方差分析Tab.5 Analysis of y1 variance

表5中，偏差平方和表示因素對(duì)試驗(yàn)結(jié)果產(chǎn)生的影響，e表示試驗(yàn)誤差，利用因素偏差平方和相對(duì)于誤差平方和來(lái)衡量因素的效應(yīng)，當(dāng)因素偏差平方和較小時(shí)，則可以認(rèn)為該因素可以忽略不計(jì)[18]。取檢驗(yàn)水平α為0.05和0.01，查F分布臨界值表可知，F(xiàn)0.05(3,3)=9.28，F(xiàn)0.01(3,3)=29.46。當(dāng)F大于F0.01(3,3)=29.46時(shí)，稱為因素高度顯著，記為**；當(dāng)F大于F0.05(3,3)=9.28時(shí)，稱為因素顯著，記為*；當(dāng)F小于F0.05(3,3)=9.28時(shí)，稱為因素不顯著。

綜上所述，A、B、C和D四個(gè)因素對(duì)封頭沖擊力有顯著性影響，其中B、C、D三個(gè)為高度顯著。根據(jù)表6和表7的分析，防波板薄膜應(yīng)力和彎曲應(yīng)力產(chǎn)生顯著性影響的是A、B和D。選擇顯著性影響因素作為篩選因子評(píng)判標(biāo)準(zhǔn)可以有效地減少工作量，為后續(xù)的多目標(biāo)優(yōu)化縮小范圍。

然而站上演講臺(tái)，一切都不一樣了。3個(gè)多小時(shí)的論壇上，David侃侃而談，詳述新西蘭葡萄酒的趨勢(shì)與變化；兩個(gè)半小時(shí)的大師班，David悉心講解新西蘭每個(gè)產(chǎn)區(qū)、每款葡萄酒的特色與亮點(diǎn)，大師班的開(kāi)始高呼三聲毛利語(yǔ)，激起了整個(gè)課堂的氣氛；1個(gè)多小時(shí)站在品鑒會(huì)的入口背景板處與葡萄酒愛(ài)好者合照、交流，聚光燈的照射下汗水一次次浸濕手帕；4個(gè)小時(shí)的晚宴上，致開(kāi)幕詞不忘感謝活動(dòng)的每一位組織者……晚上11點(diǎn)，第十屆金樽獎(jiǎng)?lì)C獎(jiǎng)典禮落下帷幕，David穿過(guò)走廊，慢慢地走回房間。留給他與廣州相處的時(shí)間不多，第二天的中午他又要登上18個(gè)小時(shí)的飛機(jī)，回到新西蘭，開(kāi)始新一周的工作。

表6 y2方差分析Tab.6 Analysis of y2 variance

表7 y3方差分析Tab.7 Analysis of y3 variance

3.3 構(gòu)建回歸方程

為了使回歸方程可以精確高效地預(yù)測(cè)非試驗(yàn)點(diǎn)的響應(yīng)值，回歸方程的基本形式采用:

式中：n為設(shè)計(jì)變量數(shù)目；ai、aj和aij為多項(xiàng)式系數(shù)。

構(gòu)造封頭沖擊力（y1）、防波板薄膜應(yīng)力（y2）和防波板彎曲應(yīng)力（y3）的回歸方程表達(dá)式如下（保留2位有效數(shù)字）：

由于回歸方程的擬合是一種近似求解，試驗(yàn)數(shù)據(jù)和響應(yīng)面回歸方程的選擇都會(huì)給擬合的結(jié)果造成影響。因此，選用復(fù)相關(guān)系數(shù)評(píng)價(jià)回歸方程的擬合精度，復(fù)相關(guān)系數(shù)越接近1則說(shuō)明擬合誤差越小?；貧w方程y1、y2和y3的復(fù)相關(guān)系數(shù)分別為0.990、0.969和0.981。可以看出，回歸方程的擬合滿足精度要求。

3.4 多目標(biāo)優(yōu)化及帕累托（Pareto）最優(yōu)解

多目標(biāo)優(yōu)化問(wèn)題不存在最優(yōu)解，而是在相互矛盾的目標(biāo)函數(shù)和限定條件下組成的一系列解的集合，往往一個(gè)目標(biāo)函數(shù)的增加是以另一個(gè)目標(biāo)函數(shù)的降低為代價(jià)，則稱這些解為帕累托（Pareto）最優(yōu)解。因此，在求得帕累托最優(yōu)解后，要根據(jù)工程實(shí)際問(wèn)題從眾多最優(yōu)解中挑選出一個(gè)或幾個(gè)最終解。

本文采用MATLAB優(yōu)化工具箱中的基于遺傳算法的gamultiobj函數(shù)進(jìn)行優(yōu)化，gamultiobj函數(shù)所采用的算法是基于NSGA-Ⅱ改進(jìn)的一種多目標(biāo)優(yōu)化算法。本文以求得的3個(gè)回歸方程為目標(biāo)函數(shù)，以4個(gè)因素的水平上下限為約束條件，求這3個(gè)目標(biāo)函數(shù)同時(shí)較小時(shí)因素的數(shù)值。由于防波板線性化結(jié)果薄膜應(yīng)力應(yīng)小于1.5倍許用應(yīng)力，彎曲應(yīng)力應(yīng)小于3倍許用應(yīng)力，以此為限定條件，對(duì)Pareto解集進(jìn)行篩選，優(yōu)選出的Pareto非劣解集如表8所示。

表8 多目標(biāo)優(yōu)化Pareto非劣解集Tab.8 Pareto noninferiority solutions

2個(gè)Pareto非劣解中，第一組解集具有較小的封頭沖擊力，相對(duì)而言y2和y3的值則偏大。第二組解集的y2和y3較小，但是封頭沖擊力較大。由于封頭采用的是形狀系數(shù)，需要在平板封頭、錐形封頭和橢圓形封頭三種情況下進(jìn)行驗(yàn)證，驗(yàn)證結(jié)果如表9所示。

表9 多目標(biāo)優(yōu)化驗(yàn)證結(jié)果Tab.9 Verification results

經(jīng)過(guò)對(duì)比考慮選定1號(hào)解為最優(yōu)解，防波效果最好且能夠滿足強(qiáng)度要求。1號(hào)解的封頭沖擊力為32 328.1 N，相對(duì)于未優(yōu)化之前降低了19.1%。圖5為優(yōu)化后線性化應(yīng)力沿厚度方向分布圖?？梢钥闯?，經(jīng)過(guò)優(yōu)化防波板薄膜應(yīng)力和彎曲應(yīng)力均得到降低，優(yōu)化度分別為12.7%和33.4%。優(yōu)化得到了良好的效果，對(duì)以后工程設(shè)計(jì)具有借鑒意義。

圖5 線性化應(yīng)力Fig.5 Linearized stress

4 結(jié) 論

本文綜合運(yùn)用流固耦合方法、正交試驗(yàn)設(shè)計(jì)、響應(yīng)面回歸方程和多目標(biāo)優(yōu)化對(duì)防波板布置、防波板形狀和防波板橫截面等4個(gè)因素進(jìn)行了分析，得到的主要結(jié)論如下：

（1）在2g加速度作用下的燃料罐壁面沖擊力先隨充裝率的增加而增加，當(dāng)充裝率大于0.85后，隨充裝率的增加而減小，即當(dāng)充裝率為0.85時(shí)，燃料罐壁面受到的沖擊力最大。

（2）通過(guò)方差分析得到各因素對(duì)防波效果和結(jié)構(gòu)強(qiáng)度的敏感程度。其中，四個(gè)因素對(duì)封頭沖擊力有顯著性影響，防波板距封頭安裝距離、防波板底邊距中心線距離和防波板形狀三個(gè)因素會(huì)對(duì)防波板薄膜應(yīng)力和彎曲應(yīng)力產(chǎn)生顯著性影響。

（3）通過(guò)多目標(biāo)優(yōu)化得到最優(yōu)結(jié)構(gòu)參數(shù)組合，與原有防波板進(jìn)行對(duì)比，封頭沖擊力、防波板薄膜應(yīng)力和彎曲應(yīng)力分別減少了19.1%、12.7%和33.4%，具有良好的優(yōu)化效果。