(1.大連理工大學(xué) 船舶工程學(xué)院，遼寧 大連 116024；2.中國核動力研究設(shè)計院，成都 610041)

LNG儲罐的蒸發(fā)率預(yù)報對于儲罐設(shè)計乃至船舶設(shè)計十分重要。利用有限元方法對LNG液貨維護系統(tǒng)的溫度場數(shù)值模擬進而得到蒸發(fā)率可以滿足工程設(shè)計精度需求，但是罐體設(shè)計方案修改后需要重構(gòu)有限元模型來完成蒸發(fā)率計算，建模成本較高。目前對于LNG蒸發(fā)率及其影響因子的研究逐漸增多，但是模型考慮的漏熱要素不全，因此適用范圍相對較小[1-6]。針對使用較廣泛的獨立雙C型罐體，分析儲罐溫度場和蒸發(fā)率的計算原理，根據(jù)實際尺寸數(shù)據(jù)確定40 500種組合方案，編寫ANSYS參數(shù)化建模腳本，利用MATLAB和ANSYS的協(xié)同工作實現(xiàn)參數(shù)化建模并同時計算出40 500種組合方案的漏熱量，基于均值分析法(analysis of means，ANOM)定量分析各因子的效應(yīng)并據(jù)此建立BP神經(jīng)網(wǎng)絡(luò)模型，實現(xiàn)雙C型獨立罐蒸發(fā)率的快速預(yù)測。

1 LNG儲罐溫度場計算控制方程

dxdydz+σdxdydz

(1)

儲罐系統(tǒng)需考慮以下初始條件和邊界條件。

1)初始條件。

(2)

(3)

當(dāng)t=0時，有

θ(x,y,z,t)|t=0=θ0(x,y,z)

(4)

2)邊界條件。

θ(t)=q(t)

(5)

(6)

(7)

式中：θ1為儲罐周圍空氣的溫度；β為熱對流系數(shù)，主要由絕熱層邊界的粗糙度、空氣的導(dǎo)熱系數(shù)、黏滯系數(shù)、流速和流向等因子決定，通過相應(yīng)的量綱一的量數(shù)進行估算[7]。

上述偏微分方程問題轉(zhuǎn)換為下面泛函數(shù)的極值問題[8]。

(8)

使用有限元法求解上面的函數(shù)，可得各節(jié)點的溫度，基于溫度場積分可得到相應(yīng)的漏熱量Q。

2 參數(shù)化建模及蒸發(fā)率數(shù)值計算

2.1 蒸發(fā)率影響因子

T-5200型LNG儲罐凈質(zhì)量517 t，內(nèi)部可承受最大壓力為0.3 MPa，主要材料是X8Ni9，滿載載貨量為7 924.8 m3。其結(jié)構(gòu)與主要參數(shù)見圖1。

l-罐體長度;j-罐體同側(cè)半球體球心距;u-墊木與罐體夾角;h3-墊木厚度;r2-氣室半徑;h2-氣室高度;h-罐體厚度;r-罐體半徑圖1 儲罐結(jié)構(gòu)

罐體導(dǎo)熱率為α1，墊木導(dǎo)熱率為α2。由于罐體半徑r受到船的型寬和型深限制，而型寬和型深又受到港口和航道的限制，罐體厚度h的熱阻較小，對結(jié)果影響不大，同時罐體厚度h是由內(nèi)部壓力決定的，因此將罐體厚度h和罐體半徑r取為定值。根據(jù)《散裝運輸液化氣體船舶構(gòu)造與設(shè)備規(guī)范》[9]和實測海洋石油301的貨艙和貨物溫度，保溫層外部環(huán)境溫度取38.5 ℃，罐體內(nèi)壁溫度取-160 ℃作為計算依據(jù)。罐體長度l(在基準長度基礎(chǔ)上以2.8 m為間隔取值，因為2.8 m是該儲罐對應(yīng)運輸船強肋骨之間的距離，有利于分艙時艙壁位置的選取)。因子取值見表1。

表1 因子取值表

表1中共有40 500種數(shù)據(jù)組合，覆蓋了容積在5 000～10 000 m3之間的LNG液罐設(shè)計參數(shù),以此為基礎(chǔ)提出符合實際情況的罐體設(shè)計方案。

2.2 儲罐參數(shù)化建模

基于可變因子對儲罐進行參數(shù)化建模。參數(shù)化有限元模型溫度場云圖見圖2。

圖2 儲罐溫度場云圖

該參數(shù)化模型能自動在1 min內(nèi)完成不同尺度方案的模型重建工作，對于后續(xù)的敏感性分析和神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)均有較大的輔助作用。

2.3 LNG罐蒸發(fā)率計算原理

LNG罐漏熱量與蒸發(fā)率之間的關(guān)系為

(9)

式中：Q為罐體內(nèi)外溫度差異吸收的瞬時總熱量，W；γ為LNG汽化熱，520.77 kJ/kg；V為儲罐容積，m3；ρ為LNG密度，kg/m3。

LNG罐自然蒸發(fā)率計算常采用穩(wěn)態(tài)溫度場模型，內(nèi)外溫差使罐體內(nèi)部和外部產(chǎn)生熱量交換，使罐體外部溫度降低，由于空氣和罐體外部的對流作用，造成罐體外部溫度基本穩(wěn)定。根據(jù)前述貨物和貨艙實測溫度設(shè)置邊界條件，網(wǎng)格尺寸為0.5 m，數(shù)值模擬從外壁到內(nèi)壁的熱量輸入，即罐體漏熱量，并用于計算蒸發(fā)率。

3 數(shù)據(jù)獲取與分析

3.1 數(shù)據(jù)獲取及有效性說明

以MATLAB作為主控程序，通過編寫ANSYS參數(shù)化建模腳本，實現(xiàn)MATLAB和ANSYS的協(xié)同工作[10]，計算出40 500種組合方案的漏熱量。

3.2 基于均值分析法的因子效應(yīng)分析

采用均值分析法[11-13]，因子效應(yīng)見表2。

表2中因子的水平表示該因子取對應(yīng)值時所有漏熱量總和的平均值，因子對應(yīng)的效應(yīng)為該因子最大水平值與最小水平值之差。

表2 ANOM分析

(10)

effectfactor=max(levelfactor)-min(levelfactor)

(11)

式中:level為水平，W；factor為因子；val為因子取值；n為總數(shù)據(jù)量；effect為效應(yīng)，W；Cal為漏熱量，W。

因子效應(yīng)排名見圖3。

圖3 因子效應(yīng)排名

α1對漏熱量的影響非常明顯，約為l的3.46倍，u和h3影響較小，r2和h2對漏熱量的影響最小，分別只有α1的1.1%和0.73%，符合之前的顯著性分析結(jié)果。由此認為，減小罐體絕熱層導(dǎo)熱率α1是減少漏熱量最有效的方式，而改變氣室的半徑和高度對漏熱量幾乎沒有影響。

4 建立人工神經(jīng)網(wǎng)絡(luò)(ANN)模型

4.1 訓(xùn)練集和測試集的選擇

選擇訓(xùn)練集和測試集時需要保證這兩個集合能代表整個數(shù)據(jù)的大致分布，否則模型在最終預(yù)測結(jié)果時誤差會變大[14-15]。采用反復(fù)迭代嘗試的方法來減小訓(xùn)練集和測試集的統(tǒng)計學(xué)差異[16]。按照實踐經(jīng)驗，選取25%的數(shù)據(jù)作為測試集，剩下的75%作為訓(xùn)練集，最終2個數(shù)據(jù)集的統(tǒng)計學(xué)性質(zhì)見表3。

表3 訓(xùn)練集與測試集數(shù)學(xué)統(tǒng)計性質(zhì)

4.2 BP神經(jīng)網(wǎng)絡(luò)模型的構(gòu)造

單隱層BP神經(jīng)網(wǎng)絡(luò)可以表示任何閉區(qū)間的連續(xù)函數(shù)[17-18]，因此本文模型隱層數(shù)為1。對于隱層節(jié)m，采用如下經(jīng)驗公式。

(12)

式中：n為輸入層節(jié)點數(shù)；k為輸出層節(jié)點數(shù)；p為1～10之間的常數(shù)。

取n=8，k=1，于是m∈[4,13]。

BP神經(jīng)網(wǎng)絡(luò)的訓(xùn)練一般采用誤差反向傳播的梯度下降算法，但該方法收斂速度慢、易陷入局部極小點，采用自適應(yīng)修改學(xué)習(xí)率算法traingdx。traingdx算法在在訓(xùn)練過程中不僅會適當(dāng)變化學(xué)習(xí)率，加快網(wǎng)絡(luò)的訓(xùn)練速度，同時還有機結(jié)合動量批梯度下降算法，有效避免局部極小點。對于傳遞函數(shù)，常用的有tansig、logsig、purelin，由于數(shù)據(jù)的復(fù)雜性，考慮tansig-tansig、tansig-purelin、logsig-purelin、logsig-tansig 4種組合。

為確定最佳的隱節(jié)點數(shù)和傳遞函數(shù)組合，在每種傳遞函數(shù)組合下，隱節(jié)點數(shù)m∈[4,13]的范圍內(nèi)逐漸增加隱節(jié)點數(shù)，用同一樣本訓(xùn)練，從中確定網(wǎng)絡(luò)誤差最小時對應(yīng)的傳遞函數(shù)組合與隱節(jié)點數(shù)，結(jié)果見圖4。

圖4 各組合誤差分布

由圖4可知，當(dāng)傳遞函數(shù)組合為tansig-purelin，隱層節(jié)點數(shù)為5時模型表現(xiàn)最好。

4.3 神經(jīng)網(wǎng)絡(luò)訓(xùn)練及仿真

模型訓(xùn)練參數(shù)設(shè)定值見表4,性能見表5。

表4 模型主要訓(xùn)練參數(shù)表

表5 模型性能表

由表5可見，無論在訓(xùn)練集還是測試集上，模型預(yù)測值的均方根誤差都在500 W左右，有少部分算例超過[-5%, 5%]的誤差區(qū)間，同時誤差值的范圍變化達到5 000 W的跨度。

4.4 模型改進

為進一步提高模型預(yù)測性能，將訓(xùn)練集分成不同的空間，在每個空間中單獨預(yù)測。由因子影響分析可知，罐體導(dǎo)熱率α1對漏熱量的影響最為顯著，于是按照α1將訓(xùn)練集分成3個子集并分別建立3個同結(jié)構(gòu)神經(jīng)網(wǎng)絡(luò)模型，在每個子集上面分別訓(xùn)練神經(jīng)網(wǎng)絡(luò)。最終模型中3個子模型的表現(xiàn)性能分別見表6、7、8。

表6 ANN1性能

表7 ANN2性能

表8 ANN3性能

由表6、7、8可見，模型誤差范圍由5 000 W降到2 000 W，均方根誤差由540 W降到240 W，誤差分布在[-5%，5%]以內(nèi)的比率為100%，相比與原模型，新模型的預(yù)測性能更加準確。

為驗證模型的可靠性，進行隨機參數(shù)樣本測試，與有限元數(shù)值計算值進行對比見表9，模型預(yù)測誤差在[-5%, 5%]內(nèi)的比率為100%。

表9 隨機參數(shù)樣本測試結(jié)果

4.5 基于PSO的方案優(yōu)化實例

某設(shè)定艙容為8 000 m3的獨立C型雙圓筒LNG罐，初始參數(shù)見表10，有限元數(shù)值計算得到其漏熱量為89 748 W。

對球心距j進行優(yōu)化以減少漏熱量，并且要求j在5～10 m之間。分別在預(yù)測模型和有限元數(shù)值計算的基礎(chǔ)上使用粒子群算法(particle swarm optimization，PSO)進行數(shù)值優(yōu)化[17-18]。結(jié)果見表11，其中T為消耗的計算時間。

表10 LNG儲罐初始參數(shù)

注：表中所有因子單位同表9。

表11 結(jié)果對比分析

由表11可見，2種方法得到的j值基本相等，但數(shù)值計算消耗的時間是預(yù)測模型的757倍。同時優(yōu)化前與優(yōu)化后相比漏熱量減少了10 549 W，對于航期較長的LNG船舶有巨大的經(jīng)濟效益。

5 結(jié)論

1)漏熱因子按其影響程度的排序為絕熱層導(dǎo)熱率>儲罐長度>罐體圓心距>墊木導(dǎo)熱率>墊木與罐體角度>墊木厚度>氣室半徑>氣室厚度。其中絕熱層導(dǎo)熱率影響最為顯著，是罐體長度的3.46倍，氣室半徑和氣室高度對蒸發(fā)率的影響很小，基本可以忽略。

2)按絕熱層導(dǎo)熱率α1將數(shù)值實驗結(jié)果數(shù)據(jù)空間劃分為3個部分，建立3個同結(jié)構(gòu)的神經(jīng)網(wǎng)絡(luò)模型，整體的預(yù)測模型誤差范圍縮小60%，均方根誤差降低56%，數(shù)據(jù)集上誤差在 [-5%, 5%]的比率為100%。

3)對于體積在范圍在5 000～10 000 m3內(nèi)的獨立C型雙圓筒LNG罐，本文模型能比較準確地預(yù)測其漏熱量。用建立的預(yù)測模型計算蒸發(fā)率的時間成本僅為基于參數(shù)化自動建模的有限元方法的1/700，適合早期方案優(yōu)化工作。