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

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

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

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

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

1 數值模擬模型

1.1 控制方程

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

連續性方程:

(1)

動量方程:

(2)

組分輸運方程:

(3)

能量方程:

k方程:

(5)

ε方程:

(7)

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

1.2 燃燒模型

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

反應1:

(8)

反應2:

(9)

式中：CnHm為反應物；O2為氧化劑；Fmid為中間產物，Fend為最終產物；Q1、Q2為反應釋放的熱量。

1.3 固體控制方程

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

(10)

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

1.4 耦合面控制方程

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

(11)

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

1.5 計算模型

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

表1 儲罐結構參數Tab.1 Tank structure parameters

1.5.1 流體模型

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

圖1 流體區域網格Fig.1 Fluid area grid

1.5.2固體模型

(12)

表2 模型材料參數Tab.2 Model material parameters

圖2 固體域網格Fig.2 Solid domain grid

1.6 耦合參數及動網格設置

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

(13)

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

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

1.7 模型驗證

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

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

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

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

2 結果分析與討論

2.1 耦合結果分析

2.1.1火焰發展

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

圖5 火焰發展情況Fig.5 Flame development

圖6為火焰傳播速度時序圖，火焰傳播速度采用文獻[21]的計算方法，結合不同時間點火焰結構進行對比可以看到，經過前期短暫的層流火焰燃燒后，在無約束狀態下，火焰傳播速度在壓力波與火焰的耦合激勵下迅速上升，并于300 ms時達到峰值65.05 m/s,隨后，火焰右側壁面開始干擾火焰及壓力波的傳播，在反射波的影響下，火焰傳播速度開始下降，隨著燃燒的逐步進行，儲罐內部超壓逐漸升高，從耦合條件下儲罐壁面發生變形，內部流場邊界發生變化，導致其火焰發展速度比非耦合條件下的火焰發展速度稍慢，可以看到，非耦合條件下，720 ms時儲罐內部已經燃燒完全，而耦合條件下780 ms時才完全燃燒。

圖6 耦合與非耦合條件下火焰傳播速度時序圖Fig.6 Coupling and non-coupled flame propagation velocity under the timing diagram

2.1.2超壓趨勢

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

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

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

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

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

2.1.3 耦合效應分析

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

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

2.2 儲罐結構響應分析

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

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

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

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

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

3 結 論

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

(1)非耦合條件與耦合條件下火焰結構變化趨勢基本相同，但由于耦合條件下儲罐壁面發生變形，因此耦合條件下的火焰傳播速度稍慢。耦合條件下儲罐內部完全燃燒時間比非耦合條件下慢60 ms。

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

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