胡 可， 趙 陽， 王 鋼， 郭艷軍

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

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

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

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

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

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

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

(1)

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

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

(2)

(3)

(4)

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

(5)

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

(6)

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

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

(7)

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

耦合控制方程為

(8)

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

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

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

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

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

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

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

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

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

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

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

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

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

表1 儲罐罐壁與頂蓋厚度

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

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

圖2 流體域耦合模型

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

圖4 固體域耦合模型

2.2 動網(wǎng)格模型

流固耦合分析中流場計(jì)算包含十分復(fù)雜的湍流流動與化學(xué)反應(yīng)，其計(jì)算收斂性與模擬結(jié)果的精確性均非常依賴于流場的網(wǎng)格質(zhì)量。由于組分概率密度運(yùn)輸燃燒模型的計(jì)算成本較高，流場網(wǎng)格數(shù)量不宜過多，因此需要一種計(jì)算效率與精度均可以得到保證的網(wǎng)格劃分方式。本文利用ANSYS ICEM CFD網(wǎng)格劃分工具，通過創(chuàng)建初始塊、關(guān)聯(lián)曲線以及采用O-grid方式將鋼儲罐內(nèi)部呈半圓柱形且?guī)Ч绊數(shù)牧鲌鲞M(jìn)行O型剖分后再進(jìn)行網(wǎng)格劃分，可以得到網(wǎng)格質(zhì)量的評價系數(shù)均在0.8以上的高質(zhì)量六面體結(jié)構(gòu)化網(wǎng)格(質(zhì)量評價系數(shù)越接近1，網(wǎng)格質(zhì)量越高；系數(shù)越接近0，網(wǎng)格質(zhì)量越低)。另外，網(wǎng)格尺寸對流場模擬計(jì)算也有一定影響。本文作者通過比較不同網(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提供的動網(wǎng)格模型能夠有效處理流場邊界運(yùn)動變形引起的計(jì)算域改變的流動問題。動網(wǎng)格算法主要包括鋪層法、彈性光順法以及局部重構(gòu)法。本文采用適用于四邊形、六面體或三棱柱網(wǎng)格的鋪層法，該算法可以根據(jù)流體域的擴(kuò)張或收縮來相應(yīng)地生成(拉伸)網(wǎng)格或合并(壓縮)網(wǎng)格。儲罐結(jié)構(gòu)部分采用四邊形網(wǎng)格劃分，單元網(wǎng)格尺寸為0.4 m，單元網(wǎng)格數(shù)為6 003。

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

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

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

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

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

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

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

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

圖8為起爆后鋼儲罐內(nèi)部壓力場的變化情況。可以看出，在爆炸沖擊波未到達(dá)儲罐結(jié)構(gòu)時，壓力前波陣面以爆源為中心呈褶皺的球形向外擴(kuò)張，流體耦合面沒有運(yùn)動變形；在t=35 ms左右，沖擊波到達(dá)儲罐并發(fā)生碰撞反射，此時罐內(nèi)乙炔/空氣混合物基本燃燒完全，在沖擊波作用下頂蓋與罐壁開始向外膨脹變形，由于耦合效應(yīng)，爆炸流場邊界也隨之向外膨脹；在t=45 ms左右，由于儲罐結(jié)構(gòu)塑性應(yīng)變到達(dá)失效應(yīng)變而耦合計(jì)算終止。在整個過程中，動網(wǎng)格根據(jù)流場邊界的擴(kuò)張來相應(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分別給出了儲罐結(jié)構(gòu)在不同時刻下的變形云圖與等效應(yīng)變云圖?？梢钥闯霎?dāng)爆炸沖擊波還未到達(dá)儲罐時(35 ms之前)，結(jié)構(gòu)會受到鄰近流體單元壓力上升而出現(xiàn)數(shù)量級較小的變形與應(yīng)變；當(dāng)爆炸沖擊波到達(dá)儲罐時，罐壁中間區(qū)域首先出現(xiàn)相對較大的變形與應(yīng)變；隨著爆炸沖擊波的繼續(xù)擴(kuò)張和反射，變形與塑性應(yīng)變區(qū)域由罐壁中央逐漸向儲罐頂蓋與底板發(fā)展；在43 ms時，儲罐結(jié)構(gòu)最大變形出現(xiàn)在頂蓋中心區(qū)域 (約0.7 m)，而最大等效應(yīng)變出現(xiàn)在頂蓋與罐壁的環(huán)向連接處(約0.206，十分接近失效應(yīng)變0.21)；爆炸流場繼續(xù)通過流固耦合作用于儲罐，在45 ms時，結(jié)構(gòu)因在頂蓋與罐壁環(huán)向連接處的最大等效應(yīng)變已超過失效應(yīng)變而發(fā)生破壞，此時最大變形仍位于頂蓋中心區(qū)域 (約1.03 m)；底板因與地面固接且建模時采用剛性體模型，因此在爆炸流場作用過程中無變形和應(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分別給出了不同時刻儲罐頂蓋與罐壁兩個區(qū)域中最大位移時程曲線與最大變形速度時程曲線，可以看出頂蓋結(jié)構(gòu)相比于罐壁結(jié)構(gòu)在爆炸沖擊波的作用下產(chǎn)生了較大的動力響應(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í)際工程中儲罐罐壁厚度要大于頂蓋；另一方面較平坦的頂蓋形式對于爆炸沖擊荷載更為敏感，相比之下呈圓柱形的罐壁結(jié)構(gòu)有著更好的抗爆性能。

圖11 儲罐頂蓋與罐壁最大位移時程曲線

圖12 儲罐頂蓋與罐壁最大變形速度時程曲線

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

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

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

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

圖13 儲罐結(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 鋼儲罐動力響應(yīng)的解耦分析結(jié)果

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

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

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

4 結(jié) 論

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

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

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

(a) 最大變形

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

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

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

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