蒸氣云爆炸作用下鋼儲(chǔ)罐動(dòng)力響應(yīng)的雙向流固耦合分析

胡 可， 趙 陽(yáng)， 王 鋼， 郭艷軍

(1. 浙江大學(xué) 空間結(jié)構(gòu)研究中心， 杭州 310058； 2. 四川電力設(shè)計(jì)咨詢有限責(zé)任公司， 成都 610094)

鋼儲(chǔ)罐內(nèi)部爆炸通常是由罐內(nèi)所儲(chǔ)存的石化物資揮發(fā)產(chǎn)生的可燃?xì)怏w與空氣混合后形成的蒸氣云遇明火或電火花引發(fā)的爆炸事故，其過(guò)程主要包括：蒸氣云爆炸、爆炸沖擊波傳播、沖擊波與儲(chǔ)罐結(jié)構(gòu)相互作用等。相對(duì)應(yīng)的研究包括爆炸沖擊波參數(shù)計(jì)算、儲(chǔ)罐內(nèi)表面爆炸荷載確定以及儲(chǔ)罐結(jié)構(gòu)動(dòng)力響應(yīng)分析等[1]。其中，結(jié)構(gòu)動(dòng)力響應(yīng)分析是鋼儲(chǔ)罐安全設(shè)計(jì)與抗爆優(yōu)化的重要基礎(chǔ)。目前，密閉空間內(nèi)部蒸氣云爆炸的CFD (Computational Fluid Dynamics)分析大多僅限于流場(chǎng)的模擬計(jì)算，而將密閉容器假定為剛體，未考慮流場(chǎng)與結(jié)構(gòu)的相互作用(即流固耦合效應(yīng))[2]。然而，內(nèi)部爆炸流場(chǎng)的作用會(huì)使儲(chǔ)罐結(jié)構(gòu)產(chǎn)生運(yùn)動(dòng)與變形，這些運(yùn)動(dòng)和變形又會(huì)反過(guò)來(lái)作用于內(nèi)部爆炸流場(chǎng)，進(jìn)而影響儲(chǔ)罐內(nèi)表面上爆炸荷載的大小與分布[3]。因此，求解鋼儲(chǔ)罐在內(nèi)部爆炸作用下的結(jié)構(gòu)動(dòng)力響應(yīng)實(shí)際上是一個(gè)復(fù)雜的雙向流固耦合問(wèn)題。

從對(duì)控制方程的解法上來(lái)區(qū)分，流固耦合問(wèn)題可分為強(qiáng)耦合與弱耦合。強(qiáng)耦合方法把流場(chǎng)、固體結(jié)構(gòu)及耦合作用全部構(gòu)造在一個(gè)控制方程里，在單一時(shí)間步內(nèi)對(duì)所有變量同時(shí)進(jìn)行求解，具有物理概念清晰、計(jì)算精度高等優(yōu)點(diǎn)，但由于流固控制方程性質(zhì)不同，求解難度較大。而弱耦合方法是在每一個(gè)時(shí)間步內(nèi)分別對(duì)固體結(jié)構(gòu)方程和流體方程進(jìn)行迭代求解，固體域和流體域的計(jì)算結(jié)果需要通過(guò)搭建中間數(shù)據(jù)交換平臺(tái)進(jìn)行彼此信息的映射來(lái)實(shí)現(xiàn)耦合求解。弱耦合方法能夠使用現(xiàn)有的商業(yè)軟件或程序，并可以根據(jù)實(shí)際情況的需要進(jìn)行固體域和流體域信息的傳遞，比強(qiáng)耦合方法更加容易實(shí)現(xiàn)，近年來(lái)得到了廣泛的使用[4]。Baum等[5]將弱耦合方法應(yīng)用到了武器爆炸領(lǐng)域，對(duì)爆炸以及殘片的運(yùn)動(dòng)過(guò)程進(jìn)行了模擬，過(guò)程中在流固耦合界面采用了快速插值方法和守恒技術(shù)進(jìn)行數(shù)據(jù)傳遞。Chun[6]使用ANSYS和FLUENT模擬了爆炸過(guò)程中沖擊波對(duì)柔性防護(hù)結(jié)構(gòu)的影響，研究結(jié)果表明，弱耦合方法是分析爆炸荷載下剛體或者柔性結(jié)構(gòu)響應(yīng)的很好工具。Haupt等[7]采用弱耦合方法對(duì)有熱效應(yīng)的流固耦合問(wèn)題進(jìn)行了計(jì)算，采用的CFD求解器是軟件TAU，而CSD求解器則是軟件ANSYS，數(shù)據(jù)之間的傳遞則采用軟件MpCCI。王彬等[8]采用MpCCI平臺(tái)將有限元軟件ANSYS與流體計(jì)算軟件STAR-CD結(jié)合起來(lái)，采用弱耦合方法對(duì)風(fēng)與膜結(jié)構(gòu)的耦合作用進(jìn)行了數(shù)值模擬。武岳等[9]采用弱耦合方法針對(duì)大跨度張拉結(jié)構(gòu)發(fā)展了風(fēng)振分析的數(shù)值風(fēng)洞方法，在研究中考慮了幾何非線性因素，編制了相應(yīng)的有限元分析程序MFS。薛祖杰[10]采用ANSYS 14.0的Workbench數(shù)值仿真平臺(tái)對(duì)超高層建筑的流固耦合風(fēng)振響應(yīng)進(jìn)行了數(shù)值模擬，得到了建筑表面風(fēng)壓分布與動(dòng)向響應(yīng)。

本文采用弱耦合方法，分別通過(guò)軟件FLUENT與ANSYS對(duì)內(nèi)部蒸氣云爆炸流場(chǎng)和鋼儲(chǔ)罐結(jié)構(gòu)進(jìn)行獨(dú)立求解，并通過(guò)ANSYS Workbench仿真平臺(tái)提供的System Coupling模塊進(jìn)行耦合界面的數(shù)據(jù)傳遞，實(shí)現(xiàn)蒸氣云爆炸作用下鋼儲(chǔ)罐結(jié)構(gòu)動(dòng)力響應(yīng)的雙向流固耦合分析。同時(shí)進(jìn)行解耦計(jì)算，并將耦合與解耦兩種方法進(jìn)行對(duì)比分析。

1 雙向流固耦合的基本理論與求解方法

1.1 流固耦合的基本控制方程

流固耦合問(wèn)題需要求解的耦合方程的定義域包含流體域與固體域。對(duì)于鋼儲(chǔ)罐內(nèi)部蒸氣云爆炸流場(chǎng)來(lái)說(shuō)，其需要遵循質(zhì)量守恒、動(dòng)量守恒、能量守恒和化學(xué)組分平衡方程。質(zhì)量守恒即連續(xù)性方程，其微分形式為[11]

(1)

式中：ux，uy，uz分別為x，y，z三個(gè)方向的速度分量；t為時(shí)間；ρ為流體質(zhì)量密度。

動(dòng)量守恒的本質(zhì)是滿足牛頓第二定律，式(2)～式(4)為x，y，z三個(gè)方向的動(dòng)量方程

(2)

(3)

(4)

能量守恒的本質(zhì)是熱力學(xué)第一定律，表達(dá)式為

(5)

組分質(zhì)量守恒可表示為系統(tǒng)內(nèi)某化學(xué)組分質(zhì)量對(duì)時(shí)間的變化率等于通過(guò)系統(tǒng)界面凈擴(kuò)散量與化學(xué)反應(yīng)產(chǎn)生的該組分生產(chǎn)率之和，組分方程為[12]

(6)

式中：cm為組分m的體積濃度；ρcm為該組分的質(zhì)量濃度；Dm為該組分的擴(kuò)散系數(shù)；Sm為生成率。

流體誘發(fā)固體振動(dòng)、位移的控制方程為

(7)

式中：Ms為固體質(zhì)量矩陣；Cs為阻尼矩陣；Ks為剛度矩陣；r為固體位移；τs為固體受到的應(yīng)力。

耦合控制方程為

(8)

式中：q為熱流量；T為溫度；下標(biāo)f為流體的特性；下標(biāo)s為固體的特性。

1.2 湍流燃燒模型與數(shù)值方法

湍流模型采用在關(guān)于湍動(dòng)能k方程的基礎(chǔ)上引入一個(gè)關(guān)于湍流耗散率ε方程而形成的k-ε兩方程模型[13]。該湍流模型的準(zhǔn)確性與收斂性都十分符合工程模擬計(jì)算的要求，目前被廣泛使用。

對(duì)于蒸氣云的燃燒爆炸過(guò)程，軟件FLUNET提供了預(yù)混燃燒模型、通用有限速率模型以及組分概率密度運(yùn)輸燃燒模型等多種模擬反應(yīng)的方法。其中，組分概率密度運(yùn)輸燃燒模型通過(guò)求解概率密度函數(shù)，模擬計(jì)算時(shí)對(duì)反應(yīng)速率的處理未進(jìn)行任何假設(shè)，嚴(yán)格考慮了化學(xué)反應(yīng)與湍流流動(dòng)的交互過(guò)程[14]。盡管采用此模型的計(jì)算量較大、計(jì)算時(shí)間較長(zhǎng)以及對(duì)網(wǎng)格質(zhì)量的要求交高，但從模擬精確性上來(lái)講是高于其他燃燒模型的，因此本文采用組分概率密度運(yùn)輸燃燒模型對(duì)鋼儲(chǔ)罐內(nèi)部蒸氣云爆炸過(guò)程進(jìn)行模擬。

軟件FLUENT采用有限體積法，通過(guò)求解描述每種組成物質(zhì)的流動(dòng)、擴(kuò)散和反應(yīng)源的守恒方程來(lái)模擬燃燒爆炸。對(duì)于壓力與速度耦合處理，采用半隱方法來(lái)求解壓力耦合方程的SIMPLE算法[15]。另外，爆炸流場(chǎng)的湍流流動(dòng)伴隨的化學(xué)反應(yīng)導(dǎo)致模型中質(zhì)量/動(dòng)量和物質(zhì)運(yùn)輸方程強(qiáng)烈耦合，在反應(yīng)流中獲得收斂解比較困難，本文采用亞松弛方法來(lái)處理這類(lèi)耦合問(wèn)題。

1.3 流固耦合界面的網(wǎng)格映射與數(shù)據(jù)傳遞

流固耦合界面是流場(chǎng)與結(jié)構(gòu)的交界面，由于弱耦合算法中的流固耦合作用只發(fā)生在此交界面上，所以耦合界面上的數(shù)據(jù)傳遞是雙向耦合計(jì)算的關(guān)鍵。在模擬過(guò)程中，結(jié)構(gòu)在流體荷載作用下發(fā)生變形或運(yùn)動(dòng)，結(jié)構(gòu)需將變形或運(yùn)動(dòng)通過(guò)耦合界面由CSD (Computational Structural Dynamics)網(wǎng)格傳遞給CFD網(wǎng)格，然后流場(chǎng)分析軟件在新的條件下求解流場(chǎng)，再把新的流體荷載通過(guò)耦合界面?zhèn)鬟f給CSD網(wǎng)格[16]。由耦合控制方程式(8)可以看出，流場(chǎng)和結(jié)構(gòu)在邊界上對(duì)應(yīng)的質(zhì)點(diǎn)之間應(yīng)該滿足作用力平衡條件與位移協(xié)調(diào)條件。然而，CSD和CFD對(duì)各自計(jì)算區(qū)域的網(wǎng)格劃分要求不同(對(duì)于鋼儲(chǔ)罐內(nèi)部蒸氣云爆炸來(lái)說(shuō)，流體網(wǎng)格比結(jié)構(gòu)網(wǎng)格要密一些)。這就會(huì)導(dǎo)致耦合界面上網(wǎng)格的不匹配，在多數(shù)情況下節(jié)點(diǎn)并不重合，無(wú)法直接傳遞交換數(shù)據(jù)，因此需要采用合理有效的映射方法進(jìn)行數(shù)據(jù)傳遞[17]。在雙向流固耦合分析中，耦合界面流場(chǎng)部分網(wǎng)格上的節(jié)點(diǎn)必須映射到耦合界面固體網(wǎng)格的節(jié)點(diǎn)上，以傳遞力參數(shù)；耦合界面固體部分網(wǎng)格上的節(jié)點(diǎn)必須映射到耦合界面流體網(wǎng)格的節(jié)點(diǎn)上，以傳遞位移參數(shù)。因此，流固耦合分析中一次完整的數(shù)據(jù)交換包含兩次在耦合界面處的映射[18]。

耦合界面上不同類(lèi)型網(wǎng)格之間的數(shù)據(jù)映射可以看成是數(shù)據(jù)在不同網(wǎng)格之間的雙向插值。ANSYS Workbench中的System Coupling耦合求解模塊提供了兩種插值方式，即Profile Preserving插值法和Profile Conservative插值法[19]。對(duì)于位移等參數(shù)的傳遞，一般采用Profile Preserving插值法；而對(duì)于力、動(dòng)量等參數(shù)的傳遞，一般采用Profile Conservative插值法[20]。

2 鋼儲(chǔ)罐動(dòng)力響應(yīng)的雙向流固耦合分析

2.1 流體域模型與固體域模型

本文研究對(duì)象為工程上較為常用的5 000 m3拱頂鋼儲(chǔ)罐，直徑20 m、罐壁高17.82 m，具體幾何尺寸由中國(guó)石油化工集團(tuán)提供[21]，圖1為儲(chǔ)罐結(jié)構(gòu)示意圖。

圖1 儲(chǔ)罐結(jié)構(gòu)示意圖

圖1中：D為鋼儲(chǔ)罐內(nèi)徑；R0為拱頂半徑(R0/D=1.2)；Hw為罐壁高度；Hr為頂蓋矢高；H為鋼儲(chǔ)罐結(jié)構(gòu)總高度，表1為鋼儲(chǔ)罐罐壁與頂蓋厚度。

表1 儲(chǔ)罐罐壁與頂蓋厚度

鋼儲(chǔ)罐內(nèi)部爆炸的雙向流固耦合分析模型由兩部分組成，即流體域(內(nèi)部爆炸流場(chǎng))和固體域(鋼儲(chǔ)罐結(jié)構(gòu))。雙向耦合分析需要同時(shí)在ANSYS Workbench中的Transient Structural模塊(ANSYS Mechanical)中創(chuàng)建固體域模型和在Fluid Flow模塊(FLUENT)中創(chuàng)建流體域模型。圖2為流體域耦合模型，即鋼儲(chǔ)罐內(nèi)部密閉空間。為了描述密閉空間內(nèi)可燃?xì)怏w爆炸的主要特征并便于求解，對(duì)模型做出以下簡(jiǎn)化假設(shè)：① 儲(chǔ)罐內(nèi)部無(wú)緊密排列的障礙物(如設(shè)備、隔板等)，也無(wú)儲(chǔ)液(LNG (Liquefied Natural Gas)、LPG (Liquefied Petroleum Gas)等)；② 點(diǎn)火前儲(chǔ)罐內(nèi)部是均勻的可燃?xì)怏w與空氣的混合物，處于常溫靜止?fàn)顟B(tài)；③ 可燃?xì)怏w、燃燒產(chǎn)物以及空氣均滿足理想氣體狀態(tài)方程；④ 燃燒氣體混合物的比熱容隨溫度變化，滿足混合規(guī)則；⑤ 可燃?xì)怏w燃燒反應(yīng)過(guò)程是單步不可逆的；⑥ 爆炸過(guò)程中不考慮儲(chǔ)罐內(nèi)部空間與外界的熱量傳遞輻射，容器壁面絕熱。起爆前，鋼儲(chǔ)罐內(nèi)部充滿處于常溫靜止?fàn)顟B(tài)的蒸氣云(即均勻混合的乙炔與空氣)，乙炔體積分?jǐn)?shù)為10%，罐內(nèi)初始溫度為300 K，初始?jí)毫?.1 MPa。起爆時(shí)刻，點(diǎn)火區(qū)域處于儲(chǔ)罐內(nèi)部中心位置，點(diǎn)火溫度為2 500 K，如圖3所示。流體域耦合模型采用k-ε湍流模型、組分概率密度運(yùn)輸燃燒模型、SIMPLE壓力/速度耦合算法、亞松弛等數(shù)值模型與算法進(jìn)行模擬。為了減少計(jì)算量，模型采用二分之一實(shí)體建模并在對(duì)稱面(B)施加對(duì)稱邊界條件；流體域頂面(C)和壁面(D)分別與儲(chǔ)罐頂蓋和罐壁接觸，是能夠隨儲(chǔ)罐運(yùn)動(dòng)變形的流體邊界而不再被考慮為剛性壁，故將其設(shè)置為流體耦合面；流體域底面(A)與固定的儲(chǔ)罐底板接觸不能產(chǎn)生變形，設(shè)置為固定流體邊界。

圖4為固體域耦合模型，模型同樣采用二分之一殼體建模并在儲(chǔ)罐對(duì)稱面(B)施加對(duì)稱邊界條件；儲(chǔ)罐頂蓋(C)和罐壁(D)分別與流體域頂面和壁面接觸，可在內(nèi)部爆炸流場(chǎng)的作用下產(chǎn)生變形，故設(shè)置為固體耦合面；而儲(chǔ)罐底板(A)置為固定邊界。

圖2 流體域耦合模型

圖3 儲(chǔ)罐內(nèi)部初始溫度場(chǎng)

圖4 固體域耦合模型

2.2 動(dòng)網(wǎng)格模型

流固耦合分析中流場(chǎng)計(jì)算包含十分復(fù)雜的湍流流動(dòng)與化學(xué)反應(yīng)，其計(jì)算收斂性與模擬結(jié)果的精確性均非常依賴于流場(chǎng)的網(wǎng)格質(zhì)量。由于組分概率密度運(yùn)輸燃燒模型的計(jì)算成本較高，流場(chǎng)網(wǎng)格數(shù)量不宜過(guò)多，因此需要一種計(jì)算效率與精度均可以得到保證的網(wǎng)格劃分方式。本文利用ANSYS ICEM CFD網(wǎng)格劃分工具，通過(guò)創(chuàng)建初始?jí)K、關(guān)聯(lián)曲線以及采用O-grid方式將鋼儲(chǔ)罐內(nèi)部呈半圓柱形且?guī)Ч绊數(shù)牧鲌?chǎng)進(jìn)行O型剖分后再進(jìn)行網(wǎng)格劃分，可以得到網(wǎng)格質(zhì)量的評(píng)價(jià)系數(shù)均在0.8以上的高質(zhì)量六面體結(jié)構(gòu)化網(wǎng)格(質(zhì)量評(píng)價(jià)系數(shù)越接近1，網(wǎng)格質(zhì)量越高；系數(shù)越接近0，網(wǎng)格質(zhì)量越低)。另外，網(wǎng)格尺寸對(duì)流場(chǎng)模擬計(jì)算也有一定影響。本文作者通過(guò)比較不同網(wǎng)格間距(0.1～0.5 m)模型的計(jì)算結(jié)果選擇了計(jì)算效率與精度均能夠得到保證的網(wǎng)格間距，即網(wǎng)格單元最大尺寸為0.3 m，網(wǎng)格單元總數(shù)為155 356。

另外，當(dāng)固體耦合面的位移參數(shù)傳遞到流體耦合面后，流體網(wǎng)格需要進(jìn)行相應(yīng)的調(diào)整。軟件FLUENT提供的動(dòng)網(wǎng)格模型能夠有效處理流場(chǎng)邊界運(yùn)動(dòng)變形引起的計(jì)算域改變的流動(dòng)問(wèn)題。動(dòng)網(wǎng)格算法主要包括鋪層法、彈性光順?lè)ㄒ约熬植恐貥?gòu)法。本文采用適用于四邊形、六面體或三棱柱網(wǎng)格的鋪層法，該算法可以根據(jù)流體域的擴(kuò)張或收縮來(lái)相應(yīng)地生成(拉伸)網(wǎng)格或合并(壓縮)網(wǎng)格。儲(chǔ)罐結(jié)構(gòu)部分采用四邊形網(wǎng)格劃分，單元網(wǎng)格尺寸為0.4 m，單元網(wǎng)格數(shù)為6 003。

2.3 材料模型與參數(shù)

相較于靜力作用，爆炸沖擊荷載作用下的儲(chǔ)罐結(jié)構(gòu)動(dòng)力響應(yīng)同時(shí)存在著塑性應(yīng)力硬化與高應(yīng)變率強(qiáng)化效應(yīng)。隨著應(yīng)變率的提高，鋼材的屈服極限和強(qiáng)度極限均有一定程度的提高而延伸率降低[22]，因而分析結(jié)構(gòu)在爆炸荷載下的彈塑性動(dòng)力響應(yīng)時(shí)，高應(yīng)變率效應(yīng)不可忽略。Cowper-Symonds模型通過(guò)在雙線性彈塑性模式基礎(chǔ)上加入冪指數(shù)形式的應(yīng)變率因子來(lái)縮放屈服應(yīng)力[23]，即同時(shí)考慮了塑性應(yīng)力硬化和應(yīng)變率強(qiáng)化效應(yīng)，材料本構(gòu)方程為

(9)

表2 C-S材料模型參數(shù)

值得注意的是，材料本構(gòu)與參數(shù)可以在ANSYS Workbench中的Engineering Data模塊進(jìn)行設(shè)置定義，如圖5所示。其中，A為靜屈服壓力σ0，B為塑性硬化模量Ep，n為塑性硬化指數(shù)，D為應(yīng)變率參數(shù)常量C，q為應(yīng)變率參數(shù)常量p。

2.4 雙向流固耦合的求解流程與解算設(shè)置

雙向流固耦合需將每一時(shí)間步的流場(chǎng)計(jì)算結(jié)果加載到結(jié)構(gòu)引起變形后，又反過(guò)來(lái)重新計(jì)算流場(chǎng)，如此反復(fù)迭代下去，如圖6所示。在ANSYS Workbench提供的雙向耦合分析模塊System Coupling中可以設(shè)置流固耦合面的數(shù)據(jù)交換(力、位移)、計(jì)算步長(zhǎng)、流場(chǎng)與固體的解算順序、求解總時(shí)長(zhǎng)等，如圖7所示。由于每一時(shí)間步的迭代均需要流體域、固體域以及數(shù)據(jù)傳遞的計(jì)算收斂，因此雙向耦合分析過(guò)程相當(dāng)復(fù)雜，計(jì)算量十分龐大，對(duì)計(jì)算設(shè)備性能要求較高，計(jì)算耗時(shí)長(zhǎng)，模擬時(shí)需要綜合考慮收斂性與計(jì)算成本的問(wèn)題。

圖5 ANSYS Workbench的Engineering Data模塊中的鋼儲(chǔ)罐材料參數(shù)

圖6 雙向流固耦合求解流程

(a) 流體域力傳遞給固體域

(b) 固體域位移傳遞給流體域

2.5 雙向流固耦合分析結(jié)果

圖8為起爆后鋼儲(chǔ)罐內(nèi)部壓力場(chǎng)的變化情況?？梢钥闯?，在爆炸沖擊波未到達(dá)儲(chǔ)罐結(jié)構(gòu)時(shí)，壓力前波陣面以爆源為中心呈褶皺的球形向外擴(kuò)張，流體耦合面沒(méi)有運(yùn)動(dòng)變形；在t=35 ms左右，沖擊波到達(dá)儲(chǔ)罐并發(fā)生碰撞反射，此時(shí)罐內(nèi)乙炔/空氣混合物基本燃燒完全，在沖擊波作用下頂蓋與罐壁開(kāi)始向外膨脹變形，由于耦合效應(yīng)，爆炸流場(chǎng)邊界也隨之向外膨脹；在t=45 ms左右，由于儲(chǔ)罐結(jié)構(gòu)塑性應(yīng)變到達(dá)失效應(yīng)變而耦合計(jì)算終止。在整個(gè)過(guò)程中，動(dòng)網(wǎng)格根據(jù)流場(chǎng)邊界的擴(kuò)張來(lái)相應(yīng)地拉伸(或生成)網(wǎng)格，以適應(yīng)流體計(jì)算區(qū)域的變化，確保流固交界面始終耦合且能夠有效傳遞交換數(shù)據(jù)。

(a) t=10 ms

(b) t=30 ms

(c) t=35 ms

(d) t=40 ms

(e) t=43 ms

(f) t=45 ms

圖9與圖10分別給出了儲(chǔ)罐結(jié)構(gòu)在不同時(shí)刻下的變形云圖與等效應(yīng)變?cè)茍D?？梢钥闯霎?dāng)爆炸沖擊波還未到達(dá)儲(chǔ)罐時(shí)(35 ms之前)，結(jié)構(gòu)會(huì)受到鄰近流體單元壓力上升而出現(xiàn)數(shù)量級(jí)較小的變形與應(yīng)變；當(dāng)爆炸沖擊波到達(dá)儲(chǔ)罐時(shí)，罐壁中間區(qū)域首先出現(xiàn)相對(duì)較大的變形與應(yīng)變；隨著爆炸沖擊波的繼續(xù)擴(kuò)張和反射，變形與塑性應(yīng)變區(qū)域由罐壁中央逐漸向儲(chǔ)罐頂蓋與底板發(fā)展；在43 ms時(shí)，儲(chǔ)罐結(jié)構(gòu)最大變形出現(xiàn)在頂蓋中心區(qū)域 (約0.7 m)，而最大等效應(yīng)變出現(xiàn)在頂蓋與罐壁的環(huán)向連接處(約0.206，十分接近失效應(yīng)變0.21)；爆炸流場(chǎng)繼續(xù)通過(guò)流固耦合作用于儲(chǔ)罐，在45 ms時(shí)，結(jié)構(gòu)因在頂蓋與罐壁環(huán)向連接處的最大等效應(yīng)變已超過(guò)失效應(yīng)變而發(fā)生破壞，此時(shí)最大變形仍位于頂蓋中心區(qū)域 (約1.03 m)；底板因與地面固接且建模時(shí)采用剛性體模型，因此在爆炸流場(chǎng)作用過(guò)程中無(wú)變形和應(yīng)變。

(a) t=10 ms

(b) t=30 ms

(c) t=35 ms

(d) t=40 ms

(e) t=43 ms

(f) t=45 ms

(a) t=10 ms

(b) t=30 ms

(c) t=35 ms

(d) t=40 ms

(e) t=43 ms

(f) t=45 ms

圖11、圖12分別給出了不同時(shí)刻儲(chǔ)罐頂蓋與罐壁兩個(gè)區(qū)域中最大位移時(shí)程曲線與最大變形速度時(shí)程曲線，可以看出頂蓋結(jié)構(gòu)相比于罐壁結(jié)構(gòu)在爆炸沖擊波的作用下產(chǎn)生了較大的動(dòng)力響應(yīng)；頂蓋區(qū)域最大變形(1.03 m左右)約為罐壁區(qū)域最大變形(0.56 m左右)1.8倍；頂蓋區(qū)域最大變形速度(230 m/s左右)約為罐壁區(qū)域最大變形速度(140 m/s左右)1.6倍。這是因?yàn)椋环矫鎸?shí)際工程中儲(chǔ)罐罐壁厚度要大于頂蓋；另一方面較平坦的頂蓋形式對(duì)于爆炸沖擊荷載更為敏感，相比之下呈圓柱形的罐壁結(jié)構(gòu)有著更好的抗爆性能。

圖11 儲(chǔ)罐頂蓋與罐壁最大位移時(shí)程曲線

圖12 儲(chǔ)罐頂蓋與罐壁最大變形速度時(shí)程曲線

3 鋼儲(chǔ)罐動(dòng)力響應(yīng)的解耦分析

3.1 爆炸荷載的分區(qū)加載

對(duì)于內(nèi)部爆炸作用下儲(chǔ)罐動(dòng)力響應(yīng)的研究，按是否考慮爆炸流場(chǎng)和儲(chǔ)罐結(jié)構(gòu)間的流固耦合效應(yīng)，可分為解耦方法和耦合方法。雙向流固耦合分析相較于解耦分析更符合儲(chǔ)罐內(nèi)爆的實(shí)際情況且所得計(jì)算結(jié)果更為精確，但其建模與計(jì)算過(guò)程也有如下一些不足：① 雙向耦合分析模型同時(shí)包含流體域模型、固體域模型以及耦合面數(shù)據(jù)傳遞模塊，建模過(guò)程十分繁瑣，對(duì)模型的調(diào)試也非常耗時(shí)；② 模型同時(shí)包含內(nèi)部爆炸流場(chǎng)網(wǎng)格與儲(chǔ)罐結(jié)構(gòu)網(wǎng)格，其網(wǎng)格數(shù)量龐大，對(duì)計(jì)算機(jī)硬件有較高的要求且計(jì)算耗時(shí)長(zhǎng)；③ 對(duì)于每一時(shí)間步的計(jì)算，固體域、流體域以及耦合面的數(shù)據(jù)傳遞都需要分別迭代計(jì)算收斂才能推進(jìn)到下一時(shí)間步，而內(nèi)部爆炸流場(chǎng)涉及到復(fù)雜的化學(xué)反應(yīng)并伴隨湍流流動(dòng)，儲(chǔ)罐結(jié)構(gòu)在爆炸沖擊波作用下的動(dòng)力響應(yīng)也涉及到幾何與材料的非線性，因此同時(shí)收斂相對(duì)困難，計(jì)算成本非常高；④ 由于耦合界面的存在，儲(chǔ)罐結(jié)構(gòu)的大變形或者失效破壞會(huì)引起流體耦合面的畸變，形成流體單元的負(fù)體積而導(dǎo)致數(shù)據(jù)傳遞失敗、模擬計(jì)算終止。此外，也有一些學(xué)者的研究結(jié)果表明殼壁變形對(duì)內(nèi)部爆炸流場(chǎng)影響小于1%，多數(shù)情況下解耦方法可得到近似程度較高的爆炸流場(chǎng)。作者在文獻(xiàn)[21]中將拱頂鋼儲(chǔ)罐結(jié)構(gòu)考慮為剛性壁，利用軟件FLUENT建立了能夠模擬內(nèi)部蒸氣云爆炸過(guò)程的CFD模型(所采用的湍流燃燒模型與數(shù)值方法與本文一致)，獲得了罐壁與頂蓋上任意位置處完整的爆炸荷載時(shí)程曲線。因此，本節(jié)將基于不考慮流固耦合效應(yīng)的內(nèi)部爆炸荷載對(duì)鋼儲(chǔ)罐結(jié)構(gòu)的動(dòng)力響應(yīng)進(jìn)行解耦分析。

由于中心起爆的爆炸荷載在儲(chǔ)罐內(nèi)表面對(duì)稱分布，可以將儲(chǔ)罐頂蓋與罐壁區(qū)域分別呈圓環(huán)狀劃分，圖13給出了四分之一儲(chǔ)罐結(jié)構(gòu)內(nèi)表面荷載分區(qū)情況，共18個(gè)荷載區(qū)域。在軟件FLUENT的計(jì)算過(guò)程中選擇每個(gè)區(qū)域內(nèi)中點(diǎn)處布置壓力觀測(cè)點(diǎn)(共18個(gè))，用于提取對(duì)應(yīng)位置的爆炸荷載時(shí)程曲線，然后均勻施加于該區(qū)域表面。圖14給出了其中四條典型的荷載分區(qū)曲線。在確定了各個(gè)分區(qū)的爆炸荷載時(shí)程曲線后，需要將其加載于鋼儲(chǔ)罐有限元模型。利用有限元前后處理軟件LS-PREPOST分別創(chuàng)建18個(gè)殼單元集合(關(guān)鍵字*SET_SHELL_LIST)與18條荷載分區(qū)曲線(關(guān)鍵字*DEFINE_CURVE)，并且一一對(duì)應(yīng)加載(關(guān)鍵字*LOAD_SHELL_SET)。圖15給出了在儲(chǔ)罐頂蓋與罐壁上分區(qū)加載18條內(nèi)部爆炸荷載時(shí)程曲線的過(guò)程，圖中選擇加載的區(qū)域?yàn)楹谏?/p>

圖13 儲(chǔ)罐結(jié)構(gòu)內(nèi)表面荷載分區(qū)

(a) 荷載分區(qū)曲線R1

(b) 荷載分區(qū)曲線R8

(c) 荷載分區(qū)曲線S5

(d) 荷載分區(qū)曲線S10

(a) 加載R1～R4

(b) 加載R1～R6

(c) 加載R1～R8

(a) 加載S1～S2

(b) 加載S1～S6

(c) 加載S1～S10

3.2 鋼儲(chǔ)罐動(dòng)力響應(yīng)的解耦分析結(jié)果

在有限元前后處理軟件LS-PREPOST中進(jìn)行了內(nèi)部爆炸荷載的分區(qū)加載后，保存關(guān)鍵字文件并交付LS-DYNA求解器進(jìn)行解耦計(jì)算。圖16和圖17分別給出了解耦分析時(shí)儲(chǔ)罐的變形云圖與等效塑性應(yīng)變?cè)茍D。與圖9和圖10對(duì)比可知，采用解耦分析所獲得的儲(chǔ)罐變形情況和塑性區(qū)域發(fā)展過(guò)程與耦合分析結(jié)果基本一致；隨著沖擊波的作用，變形與塑性應(yīng)變區(qū)域首先出現(xiàn)在罐壁中央并逐漸向頂蓋與底板發(fā)展；在43 ms時(shí)，結(jié)構(gòu)最大變形在頂蓋中心區(qū)域(約0.9 m)，而最大等效應(yīng)變出現(xiàn)在頂蓋與罐壁的環(huán)向連接處(約0.208，接近于失效應(yīng)變)；在45 ms時(shí)，頂蓋與罐壁的環(huán)向連接處的最大等效應(yīng)變已超過(guò)失效應(yīng)變，儲(chǔ)罐發(fā)生破壞，此時(shí)最大變形仍位于頂蓋中心區(qū)域(約1.34 m)。

由于解耦分析是將分區(qū)荷載曲線施加在儲(chǔ)罐單元上，結(jié)構(gòu)響應(yīng)的計(jì)算并不會(huì)像雙向耦合分析一樣在儲(chǔ)罐發(fā)生破壞時(shí)終止，而是刪除相應(yīng)的失效單元后繼續(xù)求解計(jì)算直到所設(shè)置的分析終結(jié)時(shí)間(即荷載時(shí)程曲線的終止時(shí)間90 ms)。但是所施加的分區(qū)荷載時(shí)程曲線后半段(45 ms以后)與實(shí)際情況中的爆炸荷載時(shí)程曲線會(huì)有較大差別，這是因?yàn)橐环矫鎯?chǔ)罐破壞后罐內(nèi)的高壓力會(huì)因破壞泄漏而逐漸下降，另一方面儲(chǔ)罐從初始狀態(tài)發(fā)展到結(jié)構(gòu)破壞的過(guò)程中已出現(xiàn)了較大的變形(耦合效應(yīng)作用下會(huì)引起爆炸荷載分布的改變)。因此，繼續(xù)求解計(jì)算所得的結(jié)構(gòu)響應(yīng)不能很好地描述儲(chǔ)罐破壞后的實(shí)際情況，故本節(jié)只給出了解耦分析中儲(chǔ)罐從初始狀態(tài)到出現(xiàn)失效破壞過(guò)程的計(jì)算結(jié)果。

(a) t=40 ms

(b) t=43 ms

(c) t=45 ms

(a) t=40 ms

(b) t=43 ms

(c) t=45 ms

3.3 解耦與耦合結(jié)果的對(duì)比分析

為了更加清楚地了解兩種分析方法所得結(jié)果的差異，選取了計(jì)算過(guò)程中儲(chǔ)罐結(jié)構(gòu)的最大位移和最大等效塑性應(yīng)變時(shí)程曲線進(jìn)行對(duì)比分析，如圖18所示。從圖中可以看出，兩種分析方法所得到的最大變形和最大等效塑性應(yīng)變時(shí)程曲線變化趨勢(shì)基本相同；當(dāng)爆炸沖擊波到達(dá)并作用于儲(chǔ)罐結(jié)構(gòu)后，解耦分析所得的儲(chǔ)罐最大位移與最大等效塑性應(yīng)變均超出雙向耦合分析結(jié)果。這是因?yàn)?，解耦分析中沒(méi)有考慮爆炸沖擊波與儲(chǔ)罐結(jié)構(gòu)的耦合效應(yīng)，而儲(chǔ)罐結(jié)構(gòu)的變形會(huì)降低沖擊波的反射超壓峰值，故解耦計(jì)算中加載于儲(chǔ)罐內(nèi)表面的爆炸荷載比實(shí)際情況的大，從而也導(dǎo)致該分析結(jié)果中儲(chǔ)罐結(jié)構(gòu)的變形與塑性應(yīng)變偏大。

4 結(jié) 論

本文利用ANSYS Workbench仿真平臺(tái)分別建立了基于軟件FLUENT計(jì)算的CFD模型和基于ANSYS計(jì)算的固體結(jié)構(gòu)力學(xué)模型，實(shí)現(xiàn)了內(nèi)部蒸氣體爆炸作用下鋼儲(chǔ)罐結(jié)構(gòu)動(dòng)力響應(yīng)過(guò)程的雙向流固耦合分析，同時(shí)進(jìn)行了解耦計(jì)算，并對(duì)兩種方法進(jìn)行了對(duì)比分析，得到以下主要結(jié)論：

(1) 利用ANSYS Workbench仿真平臺(tái)中的System Coupling模塊能夠有效地耦合流固交界面并在每一時(shí)間步準(zhǔn)確傳遞交換數(shù)據(jù)(位移、力)。

(2) 在內(nèi)部爆炸流場(chǎng)作用下，儲(chǔ)罐頂蓋與罐壁均會(huì)產(chǎn)生向外膨脹變形并在應(yīng)力較大位置出現(xiàn)塑性區(qū)；頂蓋結(jié)構(gòu)在爆炸流場(chǎng)作用下會(huì)產(chǎn)生更大的動(dòng)力響應(yīng)，相較之下呈圓柱形的罐壁結(jié)構(gòu)有著更好的抗爆性能。

(a) 最大變形

( b) 最大等效塑性應(yīng)變

(3) 儲(chǔ)罐結(jié)構(gòu)最大塑性應(yīng)變一般出現(xiàn)在頂蓋與罐壁的環(huán)向連接處，該區(qū)域會(huì)因塑性應(yīng)變最先超過(guò)失效應(yīng)變而發(fā)生破壞。

(4) 采用解耦分析所獲得的儲(chǔ)罐結(jié)構(gòu)變形情況和塑性區(qū)域發(fā)展過(guò)程與雙向流固耦合分析所得結(jié)果類(lèi)似，最大變形與失效破壞區(qū)域也與耦合分析結(jié)果一致；由于未考慮流固耦合效應(yīng)而被高估的爆炸荷載以及分區(qū)加載帶來(lái)的誤差，解耦分析求得的儲(chǔ)罐結(jié)構(gòu)變形位與塑性應(yīng)變等動(dòng)力響應(yīng)參數(shù)偏大。

(5) 解耦分析過(guò)程相對(duì)簡(jiǎn)單，收斂難度較小，計(jì)算成本較低，能夠較為方便地進(jìn)行內(nèi)部爆炸作用下的鋼儲(chǔ)罐結(jié)構(gòu)動(dòng)力響應(yīng)與薄弱危險(xiǎn)區(qū)域初步預(yù)測(cè)，在抗爆設(shè)計(jì)計(jì)算時(shí)也偏于安全保守；雙向耦合分析模型包含流體域、固體域模型以及耦合面數(shù)據(jù)傳遞模塊，建模過(guò)程繁瑣、模型調(diào)試費(fèi)時(shí)、網(wǎng)格數(shù)量龐大、迭代收斂緩慢，計(jì)算成本較高，但相比于解耦分析更符合鋼儲(chǔ)罐內(nèi)部蒸氣云爆炸的實(shí)際情況且所得計(jì)算結(jié)果更為精確。