基于流固耦合的立式拱頂儲罐油氣爆炸數(shù)值模擬

歐益宏， 李 潤， 袁廣強， 李國慶， 王世茂

(1. 陸軍勤務學院 油料系，重慶 401311; 2. 陸軍72489部隊，山東 煙臺 265301)

隨著石油化工事業(yè)的發(fā)展，儲油罐安全問題日益增多[1-3]。油罐區(qū)爆炸火災事故時有發(fā)生，容易造成重大的傷亡。因此，為確保庫區(qū)安全，需要對儲罐爆炸發(fā)生發(fā)展過程以及罐體結構響應進行分析研究。

國內(nèi)外許多學者針對儲罐內(nèi)部可燃氣體爆炸發(fā)生發(fā)展規(guī)律進行了研究，在不同容積模擬儲罐內(nèi)部進行了實驗[4-5]，對密閉空間內(nèi)部著火模式[6-7]、超壓荷載變化規(guī)律[8]進行了研究，同時采用不同燃燒模型[9-10]對不同容積[11]、罐頂結構[12-13]的大型儲罐結構內(nèi)部可燃氣體爆炸進行了CFD數(shù)值模擬,得到了不同因素對儲罐內(nèi)部火焰形態(tài)，爆炸沖擊載荷的影響。由于儲罐爆炸過程及其復雜，主要包括內(nèi)部油氣爆炸，爆炸沖擊壓力波的形成與傳播，沖擊波與儲罐結構的相互作用，目前的研究大多將壁面視為剛體，不考慮結構變形對流場發(fā)展情況的影響，事實上，儲罐爆炸過程是一個典型的流固耦合過程，針對耦合條件沖擊響應問題，胡可等[14-15]利用采用TNT當量模型對立式柱形容器內(nèi)部蒸氣云爆炸進行了數(shù)值模擬進行了研究，發(fā)現(xiàn)耦合效應對爆炸載荷的影響很小，在計算時可將結構視為剛體，Aune等[16-17]對薄金屬板進行了沖擊試驗以及數(shù)值模擬研究，發(fā)現(xiàn)耦合條件下的結構變形更小，且更符合工程實際，在計算時不能忽略耦合效應。在進行流固耦合仿真時，研究者大多采用TNT當量模型模擬蒸氣云爆炸，然而事實上，TNT爆源與可燃蒸氣云具有本質上的差異，無法有效模擬整個可燃蒸氣爆炸過程。因此，本文采用CFD模型對大型儲罐內(nèi)部爆炸進行了數(shù)值模擬，分別采用非耦合以及耦合方法對爆炸發(fā)展過程進行了計算，對兩種條件下的流場以及超壓荷載變化進行了對比分析，同時對耦合條件下的罐體結構響應進行了研究。以期為儲罐結構抗爆抑爆優(yōu)化提供一定參考。

1 數(shù)值模擬模型

1.1 控制方程

由于儲罐中油氣爆炸是一個具有強湍流流動性的復雜化學反應過程，因此在計算時采用RNGk-ε湍流模型來模擬爆炸過程，其中主要的控制方程為

連續(xù)性方程:

(1)

動量方程:

(2)

組分輸運方程:

(3)

能量方程:

k方程:

(5)

ε方程:

(7)

式中：ρ為密度；u為速度；P為壓力；e=CvT+fsHc為比內(nèi)能，E為總能量，k為湍動能，ε為湍動能耗散率；fs為第s種組分的質量分數(shù)；μt=Cμρk2/ε為湍流粘性系數(shù)；Rc為燃燒速率；Cμ取0.09；Cv為定容比熱；T為溫度；Hc為燃燒熱，Γ*=μt/(σ)*為湍流擴散系數(shù)，(σ)*為湍流普朗特數(shù)，δij為克羅內(nèi)克算子，i,j為坐標方向，C1、C2為常數(shù)，Rε為方程修正項。

1.2 燃燒模型

油氣爆炸過程為強湍流燃燒過程，本文采用有限速率/渦耗散模型來對爆炸過程進行模擬，在計算過程中，為減少計算成本，并提高模型的準確程度，將油氣燃燒化學反應模型簡化為以下兩個步驟：

反應1:

(8)

反應2:

(9)

式中：CnHm為反應物；O2為氧化劑；Fmid為中間產(chǎn)物，F(xiàn)end為最終產(chǎn)物；Q1、Q2為反應釋放的熱量。

1.3 固體控制方程

油氣爆炸超壓引發(fā)結構振動、位移的控制方程為

(10)

式中：Ms為固體質量矩陣；Cs為阻尼矩陣；KS為剛度矩陣；r為固體結構位移；τs為固體結構受到的應力。

1.4 耦合面控制方程

在流體與固體的耦合面有位移、熱流量、溫度、應力相等，

(11)

式中:q為熱流量；T為溫度；下標f為流體；s為固體。

1.5 計算模型

本文分析對象采用1 000 m3拱頂儲罐，具體數(shù)據(jù)由中石化設計院提供，由于儲罐結構具有軸對稱性，因此采用二分之一建模，計算模型分為流場與固體結構兩部分，為簡化模型，儲罐內(nèi)部爆炸強度會隨著儲液高度增加而降低[11]。因此，為進行保守計算，在模型簡化時，假設流場內(nèi)部無石油等液體以計算最大破壞效果,使其結果能用于指導實際工程應用。同時假設儲罐結構除薄殼外無隔板、加強圈等其他裝置，儲罐結構參數(shù)，如表1所示。

表1 儲罐結構參數(shù)Tab.1 Tank structure parameters

1.5.1 流體模型

如圖1，劃分網(wǎng)格時，分別采用單元最大尺寸為0.2 m、0.15 mm、0.1 mm的網(wǎng)格進行了網(wǎng)格無關性檢驗,發(fā)現(xiàn)0.15 mm網(wǎng)格與0.1 mm網(wǎng)格計算結果相差不大，同時，考慮到動網(wǎng)格變形，流體內(nèi)部模型最終采用單元最大尺寸為0.15 m的四面體網(wǎng)格進行劃分，總網(wǎng)格數(shù)量約為46萬。整個計算區(qū)域采用有限體積法對區(qū)域進行離散，壓力速度耦合采用piso算法。初始區(qū)域溫度設置為300 K，初始油氣濃度設為1.7%(化學當量比1.05[18])壓力設置為大氣壓，即p=0,在罐頂中心處設置壓力測點PT1，罐頂-壁環(huán)向連接處設置壓力測點PT2，罐底-環(huán)向連接處設置壓力測點PT3，在儲罐左側壁面近底部0.8 m處設置一個半徑為0.5 m的半球形區(qū)域，利用patch功能將該區(qū)域的的溫度設為2 000 K，模擬人孔處著火，計算區(qū)域的邊界條件設為壁面封閉、絕熱、無滲透。

圖1 流體區(qū)域網(wǎng)格Fig.1 Fluid area grid

1.5.2固體模型

(12)

表2 模型材料參數(shù)Tab.2 Model material parameters

圖2 固體域網(wǎng)格Fig.2 Solid domain grid

1.6 耦合參數(shù)及動網(wǎng)格設置

建立儲罐爆炸雙向耦合模型需要分別在流體域與固體域設置耦合面，在本文模型中，由于考慮的是儲罐內(nèi)部起爆，因此考慮將儲油罐罐頂以及環(huán)向壁面以及相對應的流場面設為流固耦合面，由于計算流場邊界存在變形問題，因此采用動網(wǎng)格模型中的彈簧光順以及網(wǎng)格重構法對流場邊界進行更新，彈簧光順的思路為將網(wǎng)格節(jié)點之間的連線近似為彈簧，通過計算節(jié)點力的平衡方程得到各節(jié)點光順后的位置。方程為

(13)

然而，只使用彈簧光順法時，當邊界位移遠大于局部網(wǎng)格尺寸時，網(wǎng)格質量會下降甚至出現(xiàn)負體積網(wǎng)格，為了解決此問題，本文采用網(wǎng)格重構法將這些超出網(wǎng)格尺寸標準的網(wǎng)格收集起來，并在該網(wǎng)格基礎上進行局部網(wǎng)格重構并重新對其進行質量評估。

流場計算采用fluent軟件，結構計算采用 transientstructure模塊進行計算，流固數(shù)據(jù)交互在workbench中的systemcoupling模塊進行。在數(shù)據(jù)交互過程中，對位移采用界面保形插值法進行傳輸，對壓力采用守恒插值進行傳輸[20]。

1.7 模型驗證

為驗證CFD模型的有效性，本文對已開展的模擬儲罐油氣爆炸試驗進行了數(shù)值模擬并進行了對比，圖3為模擬儲罐實驗臺架，儲罐底部半徑為0.5 m,全高為0.8 m,儲罐容積為625 L,其中P1、P2、P3為壓力傳感器設置位置，P1置于罐頂中部，P2置于罐體壁面中部距罐底0.4 m,P3置于罐底中部，點火點置于人孔附近，距罐底0.1 m。具體實驗細節(jié)可以參考文獻[8],圖4為不同測點實驗與仿真的超壓變化對比圖，可以看到，模擬儲罐油氣爆炸數(shù)值模擬與實驗的超壓上升以及振蕩趨勢基本一致，但是由于數(shù)值模擬中將壁面設為絕熱條件，在罐內(nèi)油氣完全燃燒后沒有熱量散失，因此超壓一直保持穩(wěn)定不變，在實際實驗中，由于儲罐壁面的熱耗散，導致超壓在達到超壓峰值以后會逐漸減弱。

圖3 模擬儲罐臺架Fig.3 Simulation tank bench

總的來說，該CFD模型能夠較為準確的模擬罐內(nèi)油氣爆炸過程，因此可將其應用于實際的立式拱頂油罐油氣爆炸中。

圖4 實驗與仿真超壓時序圖Fig.4 Experimental and simulation overvoltage timing diagram

2 結果分析與討論

2.1 耦合結果分析

2.1.1火焰發(fā)展

圖5為非耦合(即將罐壁示為剛體)以及耦合條件下儲罐內(nèi)部爆炸的火焰發(fā)展情況，可以看到，從整體上來看，兩種條件下火焰發(fā)展模式基本相同，著火起爆后，火焰以層流速度緩慢發(fā)展，由于著火位置位于容器壁面靠近底部的位置，因此火焰以半球型向四周傳播，此時，由于火焰面積較小，單位時間內(nèi)參與反應的未燃油氣較少，火焰結構發(fā)展平穩(wěn)，隨著儲罐內(nèi)部反應的進行，火焰逐步湍流化，火焰鋒面面積逐漸增大，隨著火焰與罐壁的逐漸接觸，由于罐壁面的局部擾動導致火焰湍流化加劇，火焰反應進一步加劇，隨著儲罐內(nèi)部反應的進行，火焰逐步蔓延至整個儲罐內(nèi)部。

圖5 火焰發(fā)展情況Fig.5 Flame development

圖6為火焰?zhèn)鞑ニ俣葧r序圖，火焰?zhèn)鞑ニ俣炔捎梦墨I[21]的計算方法，結合不同時間點火焰結構進行對比可以看到，經(jīng)過前期短暫的層流火焰燃燒后，在無約束狀態(tài)下，火焰?zhèn)鞑ニ俣仍趬毫Σㄅc火焰的耦合激勵下迅速上升，并于300 ms時達到峰值65.05 m/s,隨后，火焰右側壁面開始干擾火焰及壓力波的傳播，在反射波的影響下，火焰?zhèn)鞑ニ俣乳_始下降，隨著燃燒的逐步進行，儲罐內(nèi)部超壓逐漸升高，從耦合條件下儲罐壁面發(fā)生變形，內(nèi)部流場邊界發(fā)生變化，導致其火焰發(fā)展速度比非耦合條件下的火焰發(fā)展速度稍慢，可以看到，非耦合條件下，720 ms時儲罐內(nèi)部已經(jīng)燃燒完全，而耦合條件下780 ms時才完全燃燒。

圖6 耦合與非耦合條件下火焰?zhèn)鞑ニ俣葧r序圖Fig.6 Coupling and non-coupled flame propagation velocity under the timing diagram

2.1.2超壓趨勢

圖7為非耦合以及耦合條件下超壓時序圖，為進行分析，分別取罐頂中心處PT1,罐頂-壁環(huán)向連接處PT2，罐底-環(huán)向連接處PT3三個測點的超壓時程圖進行對比，可以看到，儲罐內(nèi)部不同測點超壓趨勢基本相同，均呈現(xiàn)出緩慢增長-急劇上升-壓力振蕩-穩(wěn)定不變的變化趨勢，結合圖5可以看出，由于儲罐內(nèi)部初期處于層流燃燒狀態(tài)，燃燒區(qū)域較小，因此超壓速度上升較為緩慢，隨著燃燒的進行，火焰面積逐漸變大，已燃氣體受熱膨脹，火焰鋒面前未燃氣體被壓縮，形成前驅壓力波，由于火焰前驅壓力波的傳播，導致波前大量未燃氣體被預熱，形成高溫高壓的未燃氣體區(qū)域，當火焰?zhèn)鞑ブ猎搮^(qū)域時，未燃氣體被迅速引燃，燃燒進一步加劇，導致超壓急劇上升并到達第一個超壓峰值，在到達第一個超壓峰值后，由于壓力波在火焰前峰與壁面間往復碰撞，反射，因而形成了超壓振蕩現(xiàn)象，隨著儲罐內(nèi)部可燃氣體燃燒完全，罐內(nèi)不再產(chǎn)生熱量，由于在計算時將壁面設為絕熱條件，因此罐內(nèi)超壓逐漸保持不變。

圖7 非耦合以及耦合條件下超壓時序圖Fig.7 Overcurrent timing diagram for uncoupled and coupled conditions

將不同條件下超壓關鍵參數(shù)對比分析如表3所示，將不同測點進行對比，可以看到，相比罐頂中心點PT1，環(huán)向連接處測點PT2、PT3由于壓力波的反射以及反射波的折向傳播更為劇烈，因此超壓峰值更大,反之升壓速率方面，由于PT1達到最大超壓峰值的時間更短，因此升壓速率更快，為2.92 MPa/s。在石油化工行業(yè)中，儲罐若先從底部破壞會造成巨大財產(chǎn)損失以及次生災害的發(fā)生，因此，PT2與PT3超壓參數(shù)對比是實際工程中關心的重點，可以看到頂壁連接處超壓峰值最大，為Pmax=1.51 MPa，但底壁連接處振蕩脈動更為劇烈。

表3 超壓關鍵參數(shù)對比Tab.3 Overpressure key parameters comparison

對比耦合與非耦合條件下各項超壓參數(shù)，可以看到，耦合條件下超壓峰值平均比非耦合條件下小10%左右，同時耦合條件下振蕩脈動較小，說明耦合效應能減小脈沖速度以及爆炸載荷。該結論在Aune[16-17]的研究中得到了證實。

2.1.3 耦合效應分析

為深入分析耦合效應對火焰?zhèn)鞑ヒ约氨ㄝd荷的削弱原因，本文繪制了流場矢量跡線圖并進行了對比，可以看到，一方面，耦合效應導致儲罐變形，罐內(nèi)流場空間增大，未燃氣體被稀釋，未燃氣體分子與氧分子碰撞反應機率降低，參與反應的活化分子數(shù)量減少，燃燒反應速率變慢，進而導致傳播速度變慢，燃燒速率的減慢也會導致升壓速率的減緩。同時，相對非耦合條件，耦合條件下儲罐內(nèi)部體積的增大也使得已燃氣體對未燃氣體的壓縮效應減弱，導致壓力波的傳播速度減慢。另一方面，耦合效應下，流場壁面變得更加圓滑，這對流場以及壓力波都有較大影響，由于非耦合下的壁面邊界更為尖銳，因此其對于內(nèi)部流場的干擾影響更為激烈，結合流場圖，可以看到在500 ms時，非耦合條件下，流場內(nèi)部右上角，左下角已經(jīng)形成渦旋結構，火焰鋒面湍流效應較為明顯，而此時，耦合條件下反應區(qū)內(nèi)尚無明顯的渦旋結構形成，已燃區(qū)域渦旋結構的形成會對內(nèi)部流場形成強烈的擾動，進而引發(fā)火焰鋒面失穩(wěn)導致湍流火焰加速，導致超壓峰值以及振蕩脈動進一步提升。

圖8 非耦合與耦合條件流場矢量跡線圖Fig.8 Uncoupled and coupled flow field vector trace

2.2 儲罐結構響應分析

相比非耦合條件，耦合條件分析還能得到儲罐結構響應情況，圖9為儲罐結構變形圖(0.5倍變形顯示)，結合儲罐頂部/壁面最大位移時程曲線可以看到，儲罐變形可以分為無明顯變形-頂部拱起-壁面膨脹-變形減緩四個階段，結合圖5可以看出，儲罐內(nèi)部著火爆炸初期，由于燃燒較為緩慢，爆炸產(chǎn)生的超壓較小，因此儲罐各處基本上無變形，隨著燃燒的逐漸進行，超壓逐漸上升，儲罐拱頂和壁面先后出現(xiàn)膨脹變形，400 ms時頂-壁連接處由于拱頂與壁面的膨脹變形出現(xiàn)而拉伸并出現(xiàn)褶皺變形，由圖10可以看到，儲罐不同區(qū)域變形時程曲線均顯現(xiàn)出平穩(wěn)不變-急劇增長-緩慢變化的區(qū)域，由于儲罐環(huán)向壁面較罐頂具有較好抗爆性能，因此罐頂較環(huán)向壁面變形發(fā)生更快，400 ms時頂部最大變形(0.834 m)為罐壁最大變形(0.472 m)的1.77倍，750 ms時材料變形開始減緩，但頂部最大變形(2.57 m)仍為環(huán)向壁面(2.36 m)的1.09倍，因此儲罐頂部較壁面更易被破壞。

圖9 結構變形圖Fig.9 Structural deformation diagram

圖10 頂部及環(huán)向壁面最大位移時程圖Fig.10 Top and ring wall maximum displacement time history

圖11為儲罐von Mises等效應力云圖分布情況，可以看到，爆炸初期，由于只有著火點附近有壓力變化，因此最大等效應力首先出現(xiàn)于著火點處，隨后，最大等效應力出現(xiàn)于罐壁-罐頂連接處并逐漸增大，出現(xiàn)該現(xiàn)象的原因一方面是因為該區(qū)域超壓峰值由于壓力波的反射、匯聚比其他區(qū)域更大，另一方面是由于頂部與環(huán)向區(qū)域的變形會對該位置產(chǎn)生較大的拉應力，400 ms時連接處最大應力為388 MPa,已超過材料的靜屈服應力345 MPa，材料已進入塑型區(qū)域，極易被破壞，隨著超壓的急劇上升，罐頂、儲罐環(huán)向壁面變形加劇，罐頂中心、罐壁環(huán)向區(qū)域等效應力均超過了材料的屈服極限，說明此時該區(qū)域極易被破壞。800 ms時，最大應力區(qū)域出現(xiàn)于頂-壁連接處以及底-壁連接處附近，在實際工程中，由于儲罐拱頂一般會加肋，罐壁環(huán)向會增設抗風圈以及加強圈。結合2.1.2所得結論，需要對頂-壁連接處進行弱化連接處理，使得爆炸時該處出現(xiàn)斷裂破壞進而產(chǎn)生泄壓，同時對底部進行加強，以防止油罐爆炸災害所導致的罐底開裂的重大事故發(fā)生。

圖11 儲罐爆炸應力云圖Fig.11 Tank explosion stress diagram

3 結 論

本文對1 000 m3拱頂儲罐結構內(nèi)部油氣爆炸過程進行了數(shù)值模擬并進行了模型實驗驗證，對耦合與非耦合條件下的爆炸的流場、超壓變化數(shù)據(jù)進行了對比分析，同時，對耦合條件下儲罐結構動力響應進行了分析。得到以下結論：

(1)非耦合條件與耦合條件下火焰結構變化趨勢基本相同，但由于耦合條件下儲罐壁面發(fā)生變形，因此耦合條件下的火焰?zhèn)鞑ニ俣壬月?。耦合條件下儲罐內(nèi)部完全燃燒時間比非耦合條件下慢60 ms。

(2)非耦合與耦合條件下各位置的超壓變化趨勢基本相同，均呈現(xiàn)出緩慢增長-急劇上升-壓力振蕩-穩(wěn)定不變的趨勢，其中非耦合條件下罐頂-壁連接處超壓峰值最大,為1.51 MPa，但耦合條件較非耦合條件超壓峰值低10%左右，且脈動較小，說明耦合效應能減小爆炸載荷10%左右以及部分振蕩脈動。

(3)爆炸時，儲罐頂部變形最大，最大應力區(qū)域出現(xiàn)在頂-壁，底-壁連接處，考慮到實際工程情況，需要對頂-壁連接處進行弱化連接處理，使得爆炸時該處出現(xiàn)斷裂破壞進而產(chǎn)生泄壓，同時對底部進行加強，以防止油罐爆炸災害所導致的罐底開裂的重大事故發(fā)生。