(1.大連理工大學(xué) 船舶工程學(xué)院,遼寧 大連 116024;2.中國核動力研究設(shè)計(jì)院,成都 610041)
LNG儲罐的蒸發(fā)率預(yù)報(bào)對于儲罐設(shè)計(jì)乃至船舶設(shè)計(jì)十分重要。利用有限元方法對LNG液貨維護(hù)系統(tǒng)的溫度場數(shù)值模擬進(jìn)而得到蒸發(fā)率可以滿足工程設(shè)計(jì)精度需求,但是罐體設(shè)計(jì)方案修改后需要重構(gòu)有限元模型來完成蒸發(fā)率計(jì)算,建模成本較高。目前對于LNG蒸發(fā)率及其影響因子的研究逐漸增多,但是模型考慮的漏熱要素不全,因此適用范圍相對較小[1-6]。針對使用較廣泛的獨(dú)立雙C型罐體,分析儲罐溫度場和蒸發(fā)率的計(jì)算原理,根據(jù)實(shí)際尺寸數(shù)據(jù)確定40 500種組合方案,編寫ANSYS參數(shù)化建模腳本,利用MATLAB和ANSYS的協(xié)同工作實(shí)現(xiàn)參數(shù)化建模并同時(shí)計(jì)算出40 500種組合方案的漏熱量,基于均值分析法(analysis of means,ANOM)定量分析各因子的效應(yīng)并據(jù)此建立BP神經(jīng)網(wǎng)絡(luò)模型,實(shí)現(xiàn)雙C型獨(dú)立罐蒸發(fā)率的快速預(yù)測。
dxdydz+σdxdydz
(1)
儲罐系統(tǒng)需考慮以下初始條件和邊界條件。
1)初始條件。
(2)
(3)
當(dāng)t=0時(shí),有
θ(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ù)進(jìn)行估算[7]。
上述偏微分方程問題轉(zhuǎn)換為下面泛函數(shù)的極值問題[8]。
(8)
使用有限元法求解上面的函數(shù),可得各節(jié)點(diǎn)的溫度,基于溫度場積分可得到相應(yīng)的漏熱量Q。
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é)果影響不大,同時(shí)罐體厚度h是由內(nèi)部壓力決定的,因此將罐體厚度h和罐體半徑r取為定值。根據(jù)《散裝運(yùn)輸液化氣體船舶構(gòu)造與設(shè)備規(guī)范》[9]和實(shí)測海洋石油301的貨艙和貨物溫度,保溫層外部環(huán)境溫度取38.5 ℃,罐體內(nèi)壁溫度取-160 ℃作為計(jì)算依據(jù)。罐體長度l(在基準(zhǔn)長度基礎(chǔ)上以2.8 m為間隔取值,因?yàn)?.8 m是該儲罐對應(yīng)運(yùn)輸船強(qiáng)肋骨之間的距離,有利于分艙時(shí)艙壁位置的選取)。因子取值見表1。
表1 因子取值表
表1中共有40 500種數(shù)據(jù)組合,覆蓋了容積在5 000~10 000 m3之間的LNG液罐設(shè)計(jì)參數(shù),以此為基礎(chǔ)提出符合實(shí)際情況的罐體設(shè)計(jì)方案。
基于可變因子對儲罐進(jìn)行參數(shù)化建模。參數(shù)化有限元模型溫度場云圖見圖2。
圖2 儲罐溫度場云圖
該參數(shù)化模型能自動在1 min內(nèi)完成不同尺度方案的模型重建工作,對于后續(xù)的敏感性分析和神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)均有較大的輔助作用。
LNG罐漏熱量與蒸發(fā)率之間的關(guān)系為
(9)
式中:Q為罐體內(nèi)外溫度差異吸收的瞬時(shí)總熱量,W;γ為LNG汽化熱,520.77 kJ/kg;V為儲罐容積,m3;ρ為LNG密度,kg/m3。
LNG罐自然蒸發(fā)率計(jì)算常采用穩(wěn)態(tài)溫度場模型,內(nèi)外溫差使罐體內(nèi)部和外部產(chǎn)生熱量交換,使罐體外部溫度降低,由于空氣和罐體外部的對流作用,造成罐體外部溫度基本穩(wěn)定。根據(jù)前述貨物和貨艙實(shí)測溫度設(shè)置邊界條件,網(wǎng)格尺寸為0.5 m,數(shù)值模擬從外壁到內(nèi)壁的熱量輸入,即罐體漏熱量,并用于計(jì)算蒸發(fā)率。
以MATLAB作為主控程序,通過編寫ANSYS參數(shù)化建模腳本,實(shí)現(xiàn)MATLAB和ANSYS的協(xié)同工作[10],計(jì)算出40 500種組合方案的漏熱量。
采用均值分析法[11-13],因子效應(yīng)見表2。
表2中因子的水平表示該因子取對應(yīng)值時(shí)所有漏熱量總和的平均值,因子對應(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é)果。由此認(rèn)為,減小罐體絕熱層導(dǎo)熱率α1是減少漏熱量最有效的方式,而改變氣室的半徑和高度對漏熱量幾乎沒有影響。
選擇訓(xùn)練集和測試集時(shí)需要保證這兩個(gè)集合能代表整個(gè)數(shù)據(jù)的大致分布,否則模型在最終預(yù)測結(jié)果時(shí)誤差會變大[14-15]。采用反復(fù)迭代嘗試的方法來減小訓(xùn)練集和測試集的統(tǒng)計(jì)學(xué)差異[16]。按照實(shí)踐經(jīng)驗(yàn),選取25%的數(shù)據(jù)作為測試集,剩下的75%作為訓(xùn)練集,最終2個(gè)數(shù)據(jù)集的統(tǒng)計(jì)學(xué)性質(zhì)見表3。
表3 訓(xùn)練集與測試集數(shù)學(xué)統(tǒng)計(jì)性質(zhì)
單隱層BP神經(jīng)網(wǎng)絡(luò)可以表示任何閉區(qū)間的連續(xù)函數(shù)[17-18],因此本文模型隱層數(shù)為1。對于隱層節(jié)m,采用如下經(jīng)驗(yàn)公式。
(12)
式中:n為輸入層節(jié)點(diǎn)數(shù);k為輸出層節(jié)點(diǎn)數(shù);p為1~10之間的常數(shù)。
取n=8,k=1,于是m∈[4,13]。
BP神經(jīng)網(wǎng)絡(luò)的訓(xùn)練一般采用誤差反向傳播的梯度下降算法,但該方法收斂速度慢、易陷入局部極小點(diǎn),采用自適應(yīng)修改學(xué)習(xí)率算法traingdx。traingdx算法在在訓(xùn)練過程中不僅會適當(dāng)變化學(xué)習(xí)率,加快網(wǎng)絡(luò)的訓(xùn)練速度,同時(shí)還有機(jī)結(jié)合動量批梯度下降算法,有效避免局部極小點(diǎn)。對于傳遞函數(shù),常用的有tansig、logsig、purelin,由于數(shù)據(jù)的復(fù)雜性,考慮tansig-tansig、tansig-purelin、logsig-purelin、logsig-tansig 4種組合。
為確定最佳的隱節(jié)點(diǎn)數(shù)和傳遞函數(shù)組合,在每種傳遞函數(shù)組合下,隱節(jié)點(diǎn)數(shù)m∈[4,13]的范圍內(nèi)逐漸增加隱節(jié)點(diǎn)數(shù),用同一樣本訓(xùn)練,從中確定網(wǎng)絡(luò)誤差最小時(shí)對應(yīng)的傳遞函數(shù)組合與隱節(jié)點(diǎn)數(shù),結(jié)果見圖4。
圖4 各組合誤差分布
由圖4可知,當(dāng)傳遞函數(shù)組合為tansig-purelin,隱層節(jié)點(diǎn)數(shù)為5時(shí)模型表現(xiàn)最好。
模型訓(xùn)練參數(shù)設(shè)定值見表4,性能見表5。
表4 模型主要訓(xùn)練參數(shù)表
表5 模型性能表
由表5可見,無論在訓(xùn)練集還是測試集上,模型預(yù)測值的均方根誤差都在500 W左右,有少部分算例超過[-5%, 5%]的誤差區(qū)間,同時(shí)誤差值的范圍變化達(dá)到5 000 W的跨度。
為進(jìn)一步提高模型預(yù)測性能,將訓(xùn)練集分成不同的空間,在每個(gè)空間中單獨(dú)預(yù)測。由因子影響分析可知,罐體導(dǎo)熱率α1對漏熱量的影響最為顯著,于是按照α1將訓(xùn)練集分成3個(gè)子集并分別建立3個(gè)同結(jié)構(gòu)神經(jīng)網(wǎng)絡(luò)模型,在每個(gè)子集上面分別訓(xùn)練神經(jīng)網(wǎng)絡(luò)。最終模型中3個(gè)子模型的表現(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ù)測性能更加準(zhǔn)確。
為驗(yàn)證模型的可靠性,進(jìn)行隨機(jī)參數(shù)樣本測試,與有限元數(shù)值計(jì)算值進(jìn)行對比見表9,模型預(yù)測誤差在[-5%, 5%]內(nèi)的比率為100%。
表9 隨機(jī)參數(shù)樣本測試結(jié)果
某設(shè)定艙容為8 000 m3的獨(dú)立C型雙圓筒LNG罐,初始參數(shù)見表10,有限元數(shù)值計(jì)算得到其漏熱量為89 748 W。
對球心距j進(jìn)行優(yōu)化以減少漏熱量,并且要求j在5~10 m之間。分別在預(yù)測模型和有限元數(shù)值計(jì)算的基礎(chǔ)上使用粒子群算法(particle swarm optimization,PSO)進(jìn)行數(shù)值優(yōu)化[17-18]。結(jié)果見表11,其中T為消耗的計(jì)算時(shí)間。
表10 LNG儲罐初始參數(shù)
注:表中所有因子單位同表9。
表11 結(jié)果對比分析
由表11可見,2種方法得到的j值基本相等,但數(shù)值計(jì)算消耗的時(shí)間是預(yù)測模型的757倍。同時(shí)優(yōu)化前與優(yōu)化后相比漏熱量減少了10 549 W,對于航期較長的LNG船舶有巨大的經(jīng)濟(jì)效益。
1)漏熱因子按其影響程度的排序?yàn)榻^熱層導(dǎo)熱率>儲罐長度>罐體圓心距>墊木導(dǎo)熱率>墊木與罐體角度>墊木厚度>氣室半徑>氣室厚度。其中絕熱層導(dǎo)熱率影響最為顯著,是罐體長度的3.46倍,氣室半徑和氣室高度對蒸發(fā)率的影響很小,基本可以忽略。
2)按絕熱層導(dǎo)熱率α1將數(shù)值實(shí)驗(yàn)結(jié)果數(shù)據(jù)空間劃分為3個(gè)部分,建立3個(gè)同結(jié)構(gòu)的神經(jīng)網(wǎng)絡(luò)模型,整體的預(yù)測模型誤差范圍縮小60%,均方根誤差降低56%,數(shù)據(jù)集上誤差在 [-5%, 5%]的比率為100%。
3)對于體積在范圍在5 000~10 000 m3內(nèi)的獨(dú)立C型雙圓筒LNG罐,本文模型能比較準(zhǔn)確地預(yù)測其漏熱量。用建立的預(yù)測模型計(jì)算蒸發(fā)率的時(shí)間成本僅為基于參數(shù)化自動建模的有限元方法的1/700,適合早期方案優(yōu)化工作。