狄駿皓 王琪, 郭旭 丁徐強(qiáng) 楊鵬飛

1.江蘇科技大學(xué)機(jī)械工程學(xué)院 2.南通理工學(xué)院汽車(chē)工程學(xué)院 3.張家港中集圣達(dá)因工程有限公司

LNG是由天然氣冷凍至-162 ℃液化而成,其主要成分為甲烷,燃燒后產(chǎn)生的排放物對(duì)環(huán)境污染極小,能夠有效降低NOx、SO2、CO2和PM的排放,是全球最為清潔的化石能源[1]。2022年,我國(guó)天然氣消費(fèi)總量為3 663×108m3,占全國(guó)能源消耗總量的8.9%,遠(yuǎn)低于全球天然氣消耗占全球能源消耗的24.1%,預(yù)計(jì)到2040年我國(guó)天然氣消費(fèi)量將會(huì)達(dá)到6 155×108m3[2]。未來(lái)LNG在我國(guó)的使用量將會(huì)顯著提升,存在巨大的消費(fèi)市場(chǎng)。隨著全國(guó)最大的鹽城“綠能港”儲(chǔ)罐升頂成功,意味著我國(guó)的天然氣儲(chǔ)罐正在向超大型化、現(xiàn)代化方向發(fā)展。LNG儲(chǔ)罐的開(kāi)發(fā)研究對(duì)于開(kāi)發(fā)和擴(kuò)大LNG市場(chǎng)也具有較大的促進(jìn)作用,能有效緩解國(guó)內(nèi)其他化石能源的消耗壓力。

國(guó)內(nèi)外諸多學(xué)者從理論方面對(duì)儲(chǔ)罐的結(jié)構(gòu)優(yōu)化進(jìn)行了研究。謝劍等[3]利用ANSYS軟件進(jìn)行有限元分析并采用控制變量法對(duì)穹頂所受載荷進(jìn)行多次校核計(jì)算,通過(guò)分析實(shí)驗(yàn)數(shù)據(jù)得出影響穹頂應(yīng)力和位移的因素主要有外罐直徑、穹頂曲率半徑和穹頂厚度等;黃文健[4]對(duì)穹頂進(jìn)行了全面的受力分析,使用牛頓迭代法基本原理及JAVA程序求解穹頂最低預(yù)應(yīng)力來(lái)保證穹頂安全性;季冰乙[5]通過(guò)創(chuàng)建儲(chǔ)罐有限元模型,深入分析內(nèi)罐泄漏工況下,預(yù)應(yīng)力系統(tǒng)對(duì)罐壁的內(nèi)力分布、罐壁配筋和罐壁液密性影響。得出橫向預(yù)應(yīng)力和豎向預(yù)應(yīng)力對(duì)罐壁影響的大小關(guān)系,從而對(duì)罐壁進(jìn)行結(jié)構(gòu)優(yōu)化;Santhosh等[6]在儲(chǔ)罐屈曲約束狀態(tài)下對(duì)統(tǒng)一模型進(jìn)行了拓?fù)鋬?yōu)化并找出了最佳的穹頂形狀,并且相比原有模型在質(zhì)量上減少了6%,實(shí)現(xiàn)了儲(chǔ)罐的輕量化設(shè)計(jì)。

本研究以儲(chǔ)罐直徑、穹頂曲率半徑和罐壁厚度為變量,對(duì)儲(chǔ)罐結(jié)構(gòu)的輕量化設(shè)計(jì)進(jìn)行了優(yōu)化。通過(guò)建立薄膜罐數(shù)學(xué)模型和目標(biāo)函數(shù)對(duì)主容器進(jìn)行受力分析,并利用ANSYS軟件進(jìn)行參數(shù)化設(shè)計(jì)和應(yīng)力分析,利用粒子群算法求取最優(yōu)解,并對(duì)最終結(jié)果進(jìn)行了試驗(yàn)驗(yàn)證。

1 儲(chǔ)罐結(jié)構(gòu)設(shè)計(jì)

本文研究的對(duì)象為立式圓筒平底鋼質(zhì)的LNG薄膜罐,薄膜罐的結(jié)構(gòu)簡(jiǎn)圖如圖1所示[7]。其由1個(gè)薄膜的鋼質(zhì)主容器、絕熱層以及1個(gè)混凝土罐共同組成。混凝土罐外部為工作平臺(tái)以及鋼質(zhì)走道,鋼質(zhì)穹頂下方為鋁質(zhì)吊頂,穹頂部分和吊頂部分用吊桿和吊頂加強(qiáng)筋連接,吊頂下方是玻璃纖維棉制的吊頂絕熱層。薄膜罐相較于其他儲(chǔ)罐而言,能將儲(chǔ)罐的絕緣性和氣密性作用明確分離,且罐容沒(méi)有極限。對(duì)LNG薄膜罐的穹頂結(jié)構(gòu)優(yōu)化研究更符合未來(lái)儲(chǔ)罐大型化、輕量化的發(fā)展趨勢(shì)。

LNG薄膜罐的主容器為立式圓筒鋼質(zhì)結(jié)構(gòu),本研究中的儲(chǔ)罐基于相關(guān)標(biāo)準(zhǔn)和規(guī)范設(shè)計(jì)。其各部分基本幾何參數(shù)見(jiàn)表1。

初始穹頂設(shè)計(jì)厚度的初始值按式(1)計(jì)算:

(1)

式中:S0為初始穹頂設(shè)計(jì)厚度,m;pd為主容器內(nèi)壓,MPa;Di為穹頂內(nèi)直徑,m;Sts為材料設(shè)計(jì)許用應(yīng)力,MPa;JE為焊接接頭系數(shù),無(wú)量綱;C為鋼板厚度負(fù)偏差,mm。初始設(shè)計(jì)穹頂材料厚度為6.5 mm,此時(shí)穹頂弧度為82°。

根據(jù)GB 50341-2014《立式圓筒形鋼制焊接油罐設(shè)計(jì)規(guī)范》中的定設(shè)計(jì)點(diǎn)法,初始罐壁設(shè)計(jì)厚度如表2所列。

2 優(yōu)化穹頂數(shù)學(xué)模型

2.1 建立數(shù)學(xué)模型

(2)

式中:h為矢高(BD長(zhǎng)度),m;A為穹頂面積,m2;R為曲率半徑,m。

2.2 目標(biāo)函數(shù)

穹頂?shù)慕Y(jié)構(gòu)尺寸設(shè)計(jì)取決于儲(chǔ)罐半徑、穹頂?shù)那拾霃揭约榜讽敽穸?。因?在儲(chǔ)罐半徑固定的前提下,選取最優(yōu)的曲率半徑R和鋼板厚度S,以尋求最輕穹頂質(zhì)量。S和R的數(shù)量關(guān)系如式(3)所示,取S的最小值。

(3)

式中:S為鋼板厚度,m;p1為儲(chǔ)罐設(shè)計(jì)內(nèi)壓,取0.025 MPa;φ為穹頂焊接接頭系數(shù),取0.35;[σ]t為設(shè)計(jì)溫度下許用應(yīng)力,取160 MPa;C1D為頂板鋼板負(fù)偏差,取0.3 mm;C2D為頂板腐蝕裕量,取0.3 mm。

在儲(chǔ)罐直徑和高度確定的情況下尋求罐頂?shù)淖顑?yōu)厚度以及罐壁的最優(yōu)厚度,從而得到最輕罐壁質(zhì)量[9]。目標(biāo)函數(shù)M如式(4)所示。

(4)

式中:M為主容器的最優(yōu)質(zhì)量,kg;ρ為鎳鋼材料密度,取7.20×103kg/m3;δ為罐壁厚度,m。

2.3 約束函數(shù)

根據(jù)儲(chǔ)罐設(shè)計(jì)標(biāo)準(zhǔn)GB 50341-2014規(guī)定,儲(chǔ)罐穹頂?shù)那拾霃綉?yīng)為儲(chǔ)罐直徑的0.8～1.2倍。且儲(chǔ)罐穹頂設(shè)計(jì)需要滿足式(5)所示的應(yīng)力要求。

(5)

式中:D為儲(chǔ)罐直徑,為35 m;σmax(i)為穹頂表面的最大應(yīng)力,MPa;[σi]為鎳鋼的材料許用應(yīng)力,MPa;σmax(k)為罐頂與罐壁連接處的最大應(yīng)力,MPa。

3 有限元分析參數(shù)化設(shè)計(jì)

3.1 主容器受力分析

在滿罐狀態(tài)下,儲(chǔ)罐罐壁主要受內(nèi)罐液體靜壓載荷以及地震作用,儲(chǔ)罐穹頂受到的力主要是受穹頂本身結(jié)構(gòu)的質(zhì)量產(chǎn)生的自重載荷、鋁質(zhì)吊頂以及吊頂拉桿等產(chǎn)生的附加靜載、內(nèi)載荷、雪載荷、風(fēng)載荷以及地震載荷等[10-12]。

3.1.1穹頂外載荷

穹頂及吊頂各部件質(zhì)量如表3所列。

表3 穹頂附件各部分質(zhì)量及材料kg部件名稱質(zhì)量材料鋁吊頂板1 180AL5083吊頂拉桿20 020S30408吊頂加強(qiáng)筋34 200S30408吊頂絕熱層16 800玻璃纖維棉

穹頂?shù)拿孑d荷集度指的是將穹頂所受的力以均布載荷的方式添加到整個(gè)穹頂上,可反映穹頂所受靜力的情況。穹頂所受均布載荷為1 200 N/m2。穹頂受外載荷的面載荷集度T按式(6)計(jì)算。

(6)

式中:T為面載荷集度,N/m2;g為重力加速度,取9.8 m/s2。

3.1.2風(fēng)載荷

對(duì)于儲(chǔ)罐風(fēng)載荷的計(jì)算參考GB 55001-2021《工程結(jié)構(gòu)通用規(guī)范》,儲(chǔ)罐穹頂風(fēng)載荷標(biāo)準(zhǔn)值可按式(7)計(jì)算。

(7)

式中:ωk為穹頂風(fēng)載荷標(biāo)準(zhǔn)值,N。

3.1.3地震作用

對(duì)于地震作用的受力分析方法,采用振型分解反應(yīng)譜法對(duì)其進(jìn)行分析,確定地震反應(yīng)譜后進(jìn)行儲(chǔ)罐的模型簡(jiǎn)化并初步模態(tài)分析得出最大加速度,接著計(jì)算不同狀態(tài)下需要施加的加速度大小,利用ANSYS軟件對(duì)滿罐狀態(tài)下的儲(chǔ)罐進(jìn)行豎向運(yùn)行基準(zhǔn)地震作用(operating basis earthquake, OBE)分析。滿罐狀態(tài)下的OBE最大響應(yīng)分析結(jié)果如表4所列。

表4 OBE的最大響應(yīng)分析結(jié)果地震作用質(zhì)點(diǎn)模態(tài)1模態(tài)2模態(tài)3最大加速度/(m·s-2)水平OBE(滿罐)內(nèi)罐+脈沖液體01.6881.2542.361縱向OBE(滿罐)內(nèi)罐+LNG液體1.1450.96901.500

根據(jù)GB 55001-2021選取儲(chǔ)罐罐頂受外載荷、內(nèi)載荷、風(fēng)載荷作用下的最不利載荷組合工況進(jìn)行應(yīng)力分析,穹頂所受載荷效應(yīng)組合值S如式(8)所示。

T=Gk+γQ1φC1Q1k+γQ2φC2Q2k

(8)

式中:T為穹頂載荷效應(yīng)組合值,MPa;Gk為永久載荷的標(biāo)準(zhǔn)值,N;Q1k為風(fēng)載荷的標(biāo)準(zhǔn)值,N;Q2k為地震作用下的可變載荷值,N;γ為可變載荷的分項(xiàng)系數(shù),無(wú)量綱;φ為可變載荷的載荷的組合值系數(shù),無(wú)量綱;Q1為風(fēng)載荷作用的標(biāo)準(zhǔn)值,MPa;Q2為地震作用的標(biāo)準(zhǔn)值,MPa;C1為風(fēng)載荷作用的可變值,MPa;C2為地震作用的可變值,MPa。

3.1.4罐壁靜液體靜壓載荷

不同液位高度下LNG液體對(duì)主容器側(cè)壁的壓力大小不同,兩者函數(shù)關(guān)系如式(9)所示。

p=ρg(20.80-Z)

(9)

式中:p為主容器側(cè)壁的壓力,kPa;Z為罐壁沿高度方向值,m;ρ為L(zhǎng)NG液體密度,kg/m3;g為重力加速度,N/kg。

3.2 ANSYS參數(shù)化設(shè)置

在DesignModeler模塊中繪制二維草圖尺寸以及繪制平面時(shí)對(duì)穹頂高度和頂板厚度進(jìn)行參數(shù)化設(shè)置,罐壁厚度按表2中數(shù)據(jù)進(jìn)行建模,且為減少在使用ANSYS軟件仿真時(shí)的計(jì)算時(shí)間和計(jì)算量,將儲(chǔ)罐的三維模型簡(jiǎn)化成整體的1/4進(jìn)行分析[13],在DM中進(jìn)行對(duì)稱處理。為防止在后續(xù)穹頂優(yōu)化研究過(guò)程中網(wǎng)格劃分對(duì)實(shí)驗(yàn)結(jié)果造成的影響,需要先對(duì)穹頂應(yīng)力最大處進(jìn)行網(wǎng)格無(wú)關(guān)性檢驗(yàn)[14],根據(jù)表5網(wǎng)格無(wú)關(guān)性檢驗(yàn)數(shù)據(jù)可知,當(dāng)網(wǎng)格數(shù)超過(guò)1 336個(gè)時(shí),最大等效應(yīng)力之間的誤差在2%以內(nèi),因此,本研究選用2 783個(gè)網(wǎng)格單元進(jìn)行分析。

表5 網(wǎng)格無(wú)關(guān)性檢驗(yàn)數(shù)據(jù)網(wǎng)格數(shù)量/個(gè)最大等效應(yīng)力/MPa19849.321 33657.362 78358.485 41658.2211 86457.7430 05059.10

在Workbench模塊中將靜態(tài)結(jié)構(gòu)分析、模態(tài)響應(yīng)譜分析以及風(fēng)載荷的分析3種案例進(jìn)行SRSS組合分析。接著生成dxdb文件并提取薄膜罐的最大應(yīng)力σmax。

在對(duì)罐壁進(jìn)行優(yōu)化時(shí),隨機(jī)給定一組設(shè)計(jì)變量作為罐壁厚度,此變量數(shù)據(jù)放置于可供ANSYS軟件調(diào)用的文件中,此時(shí)粒子坐標(biāo)為(δ集,R),δ集=(δ1,……δ8),R=17.5為定量。在ANSYS中進(jìn)行循環(huán)運(yùn)算,選取滿足應(yīng)力條件的5組數(shù)據(jù)繪制成表6。

表6 罐壁優(yōu)化數(shù)據(jù)(選取)mm序號(hào)δ1δ2δ3δ4δ5δ6δ7δ8113.1612.5511.009.958.056.545.103.68213.0512.5210.829.587.966.425.043.55312.8312.3010.709.427.926.384.833.40412.7312.1410.569.387.786.334.783.25512.5411.9510.549.147.756.254.703.20最優(yōu)解名義厚度13.0013.0011.0010.009.007.006.006.00

4 優(yōu)化方法設(shè)計(jì)

4.1 粒子群算法設(shè)計(jì)

粒子群算法屬于進(jìn)化算法的一種[15],是受鳥(niǎo)群覓食集群活動(dòng)的規(guī)律性啟發(fā), 鳥(niǎo)群在搜尋食物的過(guò)程中通過(guò)相互信息交換,從而判斷自身是否處于最佳位置。此算法可在限定范圍內(nèi)共享群體和個(gè)體的運(yùn)動(dòng)信息,使復(fù)雜無(wú)序問(wèn)題變得有序,并通過(guò)適應(yīng)度來(lái)確定解的優(yōu)劣性,最終通過(guò)多次迭代而獲得最優(yōu)解,具有易實(shí)現(xiàn)、精度高和收斂快等優(yōu)點(diǎn)[16]。對(duì)于確定罐體直徑的儲(chǔ)罐穹頂?shù)膬?yōu)化問(wèn)題,以穹頂曲率半徑及板材厚度為變量,以面載荷集度為目標(biāo)函數(shù),使用粒子群算法能快速有效地求得最優(yōu)解。

使用粒子群算法進(jìn)行優(yōu)化時(shí),利用無(wú)質(zhì)量粒子進(jìn)行運(yùn)算,該粒子只具有速度和位置兩個(gè)屬性,在進(jìn)行穹頂結(jié)構(gòu)優(yōu)化時(shí)粒子位置坐標(biāo)相對(duì)應(yīng)于穹頂?shù)陌宀暮穸群颓拾霃竭@兩個(gè)自變量。設(shè)計(jì)算法運(yùn)算流程如下[17]:

(1) 參數(shù)設(shè)置。設(shè)置算法的慣性權(quán)重ω=0.9,粒子維度Q=2,個(gè)體學(xué)習(xí)因子c1=2、群體學(xué)習(xí)因子c2=2,粒子群規(guī)模n=20,迭代次數(shù)K=20,最大迭代步數(shù)Smax=2。

搜索邊界Plimit=(xmin,xmax)如式(10)所示:

(10)

式中:S為板材厚度,mm;R為穹頂曲率徑,mm。最大速度vmax=(Simax,Rimax)(v為粒子速度,mm/s;Simax為第i個(gè)粒子的最大板材厚度,mm;Rimax為第i個(gè)粒子的最大曲率半徑,mm)。

(2) 初始化。初始化種群中每個(gè)粒子的位置xi=(6,28)及速度vi。

(3) 適應(yīng)度計(jì)算。計(jì)算每個(gè)粒子的適應(yīng)度值f,適應(yīng)度函數(shù)為穹頂最大應(yīng)力值。根據(jù)粒子初始位置以曲率半徑和板材厚度為依據(jù)建立參數(shù)化模型并進(jìn)行ANSYS有限元分析計(jì)算,利用ANSYS生成dxdb文件,提取最大等效應(yīng)力值σmax。

(4) 數(shù)據(jù)記錄。記錄并更新最佳曲率半徑和罐壁材料厚度組合即粒子的最佳位置xbest,種群最佳尺寸組合Pg及其適應(yīng)度值fbest。

(5) 更新粒子位置與速度如式(11)所示:

(11)

式中:rand1()、rand2()為0～1之間的隨機(jī)數(shù);z為隨機(jī)粒子,無(wú)量綱。

(6) 迭代規(guī)則。返回第3步進(jìn)行迭代,當(dāng)?shù)罄肕ATLAB讀取生成的dxdb文件后提取的最大應(yīng)力值若大于120 MPa時(shí),則跳過(guò)步驟(4)直接找出最佳尺寸組合(Sbest,Rbest)和群體最優(yōu)解。

(7) 當(dāng)達(dá)到最大迭代次數(shù)時(shí)停止迭代。輸出全局最優(yōu)解和穹頂最小等效應(yīng)力σmin,此時(shí)粒子的位置則為全局最優(yōu)尺寸組合。

在對(duì)罐壁進(jìn)行優(yōu)化時(shí),程序隨機(jī)給定一組設(shè)計(jì)變量作為罐壁厚度,此變量數(shù)據(jù)放置于可供ANSYS軟件調(diào)用的文件中,此時(shí)粒子坐標(biāo)為(δ集,R),δ集=(δ1,……δ8),R=17.5為定量。

根據(jù)上述算法步驟設(shè)計(jì),繪制如圖3所示的粒子群算法流程圖。

4.2 罐頂優(yōu)化結(jié)果

根據(jù)4.1節(jié)所示算法步驟并利用MATLAB設(shè)計(jì)算法程序[18],得出目標(biāo)函數(shù)關(guān)于穹頂曲率半徑及板材厚度的三維圖像曲線,算法迭代完成后得到粒子的最優(yōu)位置及最優(yōu)解[19-20]。

經(jīng)過(guò)多次優(yōu)化,最終得到穹頂曲率半徑為28 000 mm、穹頂板材厚度均為6 mm時(shí)為最優(yōu)解且滿足設(shè)計(jì)要求。連續(xù)5次優(yōu)化結(jié)果對(duì)比如表7所列。

表7 連續(xù)5次優(yōu)化結(jié)果對(duì)比優(yōu)化次數(shù)曲率半徑/mm板材厚度/mm表面最大法向應(yīng)力/MPa128 0006.0258.32228 0006.0058.33328 0006.0058.33428 0006.0158.32528 0006.0058.33

5 優(yōu)化后結(jié)果

5.1 優(yōu)化前后結(jié)果對(duì)比

罐壁部分優(yōu)化后實(shí)際所需厚度發(fā)生了改變,罐底部分名義厚度減小了1 mm。罐頂部分優(yōu)化后的穹頂材料厚度減小了0.5 mm,此時(shí)的穹頂半徑為28 000 mm,弧度由初始的82°變?yōu)?8°,此時(shí)穹頂?shù)馁|(zhì)量減輕了103 t,實(shí)現(xiàn)了罐體的輕量化目標(biāo)。優(yōu)化前后結(jié)果對(duì)比如表8所列,優(yōu)化后的模型法向應(yīng)力分析圖如圖4所示,滿足應(yīng)力要求。

表8 優(yōu)化前后結(jié)果對(duì)比優(yōu)化前后穹頂半徑/mm穹頂厚度/mm穹頂弧度/(°)罐底最厚部分/mm優(yōu)化前30 0906.58213.16優(yōu)化后28 0006.07812.54

5.2 優(yōu)化后實(shí)驗(yàn)結(jié)果

基于上文對(duì)罐頂結(jié)構(gòu)的優(yōu)化分析,得出此罐穹頂?shù)淖顑?yōu)曲率半徑為儲(chǔ)罐直徑的0.8倍,穹頂厚度為6 mm。此時(shí)主容器表面所受的最大應(yīng)力滿足要求,對(duì)優(yōu)化后結(jié)果進(jìn)行氣壓實(shí)驗(yàn)。氣壓試驗(yàn)中,環(huán)境溫度為20 ℃,操作介質(zhì)溫度為20 ℃的干燥空氣,使用精度等級(jí)為1.6級(jí)、量程為0～40 kPa的壓力表對(duì)儲(chǔ)罐進(jìn)行耐壓試驗(yàn),試驗(yàn)壓力為31.25 kPa,保壓時(shí)間為60 min,儲(chǔ)罐無(wú)滲漏,無(wú)可見(jiàn)異常變形,無(wú)異常響聲。進(jìn)行嚴(yán)密性試驗(yàn)時(shí),試驗(yàn)壓力為25 kPa,保壓時(shí)間為24 h,儲(chǔ)罐無(wú)滲漏,無(wú)明顯異響。符合設(shè)計(jì)要求。優(yōu)化后結(jié)果被應(yīng)用于中集圣達(dá)因山東臨沂兩萬(wàn)方薄膜罐項(xiàng)目。

6 結(jié)論

本文以LNG薄膜罐為研究對(duì)象,為得到主容器最輕質(zhì)量,對(duì)儲(chǔ)罐結(jié)構(gòu)進(jìn)行了優(yōu)化設(shè)計(jì)。設(shè)計(jì)了適合的粒子群算法并利用ANSYS和MATLAB軟件進(jìn)行分析優(yōu)化。所完成的主要工作如下:

(1) 根據(jù)相關(guān)標(biāo)準(zhǔn)對(duì)相關(guān)參數(shù)進(jìn)行初步設(shè)計(jì)。利用ANSYS對(duì)LNG儲(chǔ)罐罐頂進(jìn)行受力分析,根據(jù)薄膜罐的幾何結(jié)構(gòu)即其中的數(shù)量關(guān)系得出了目標(biāo)函數(shù),并確定了影響因子。

(2) 利用參數(shù)化、模塊化設(shè)置解決工程結(jié)構(gòu)應(yīng)力分析問(wèn)題。分析過(guò)程中在ANSYS不同模塊中進(jìn)行參數(shù)化設(shè)置并通過(guò)多組數(shù)據(jù)對(duì)罐壁厚度進(jìn)行了優(yōu)化。

(3) 將工程問(wèn)題和算法相結(jié)合。設(shè)計(jì)合適的粒子群算法,并利用MATLAB進(jìn)行編程和ANSYS提取其計(jì)算過(guò)程中的最大等效應(yīng)力,得到了罐頂最優(yōu)曲率半徑以及罐頂和罐壁的最優(yōu)材料厚度。

(4) 對(duì)優(yōu)化后的薄膜罐主容器進(jìn)行氣壓試驗(yàn)檢驗(yàn),優(yōu)化后的薄膜罐實(shí)現(xiàn)了輕量化的目標(biāo),且更具經(jīng)濟(jì)性。