林全生，高 紅，宋 波

（1. 天津大學(xué) 化工學(xué)院，天津 300350；2. 公安部天津消防研究所，天津 300381）

［關(guān)鍵字］泡沫；數(shù)值模擬；黏度；沿程阻力系數(shù)；Fluent軟件

壓縮空氣泡沫滅火劑是一類主要用于撲救固體可燃物火災(zāi)的新型高效泡沫滅火劑，具備產(chǎn)生高穩(wěn)定泡沫的性能，可用于石油化工可燃性液體危險(xiǎn)品儲(chǔ)罐火災(zāi)時(shí)的比鄰儲(chǔ)罐防護(hù)［1-4］。該類滅火劑在固定式壓縮空氣泡沫系統(tǒng)中具有突出的滅火效果。但泡沫在管道運(yùn)輸過程中，以水為基礎(chǔ)的現(xiàn)有理論難以滿足要求，需要找到適合壓縮空氣泡沫流體流動(dòng)的計(jì)算方法，為工程設(shè)計(jì)提供準(zhǔn)確依據(jù)。針對泡沫，利用單相流對其進(jìn)行研究的較多。葛曉霞等［5］通過搭建實(shí)驗(yàn)平臺(tái)來觀察泡沫在管內(nèi)的流態(tài)，給出了壓縮空氣泡沫流型的類型，并通過理論分析得出了模型參數(shù)的算法，確定了摩阻系數(shù)及壓力損失的計(jì)算方程。Sun等［6］通過實(shí)驗(yàn)發(fā)現(xiàn)泡沫流體是剪切稀化的冪率流體，并得到了摩阻系數(shù)與廣義雷諾數(shù)的關(guān)系。陳旸等［7］開展了壓縮空氣泡沫管網(wǎng)輸送實(shí)驗(yàn)，研究了泡沫流量、氣液比、泡沫種類對輸送管路壓力變化的影響及摩擦阻力損失。Hu等［8］對氣液兩相流在90°彎管中的流動(dòng)情況進(jìn)行了模擬，得到了一定氣液比下雙相流彎管壓力分布規(guī)律。Supa-Amornkul等［9］分別采用單相流和雙相流對彎管壓降進(jìn)行數(shù)值模擬，氣體體積分?jǐn)?shù)在50%以上的情況下，沿長度方向單相流壓力分布與實(shí)驗(yàn)結(jié)果差別不大，相比之下，多相流并不能很好地預(yù)測彎管壓降。林全生等［10］利用上海消防總隊(duì)的實(shí)驗(yàn)數(shù)據(jù)，對泡沫在不同管徑、不同長度的管路中的流動(dòng)阻力損失進(jìn)行了數(shù)值模擬，并得出了范寧摩擦因子與廣義雷諾數(shù)的擬合關(guān)系式。Fluent數(shù)值模擬已在工程上有很多應(yīng)用［11-14］，它在消防泡沫輸運(yùn)方面的研究還有很大的探索空間。

本工作利用Fluent數(shù)值模擬計(jì)算軟件對壓縮空氣泡沫在管道內(nèi)的流動(dòng)進(jìn)行了數(shù)值模擬研究，考慮了發(fā)泡倍數(shù)和流量的影響，并與相關(guān)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行了對比，通過擬合得出沿程阻力系數(shù)的顯性表達(dá)式，為壓縮空氣泡沫的工程計(jì)算提供指導(dǎo)。

1 數(shù)值模擬計(jì)算方法

1.1 幾何模型和網(wǎng)格劃分

基于有限元方法，根據(jù)文獻(xiàn)［7］，管路設(shè)備尺寸如圖1所示，利用Solidworks軟件建立三維幾何模型。利用ANSYS Workbench16.0中的mesh模塊對幾何模型進(jìn)行四面體網(wǎng)格劃分，網(wǎng)格數(shù)量為113萬，并對幾何模型的壁面設(shè)置膨脹層進(jìn)行邊界層加密處理，如圖2所示。

圖1 實(shí)體管路示意圖［7］Fig.1 Diagram of solid piping［7］.

圖2 幾何模型網(wǎng)格劃分Fig.2 Mesh division of geometric model.

1.2 計(jì)算模型假設(shè)

葛曉霞等［5］對壓縮空氣泡沫進(jìn)行了流型的劃分，將泡沫在管內(nèi)的流動(dòng)形態(tài)按含氣率從低到高的順序依次劃分為泡狀流、塞狀流、層狀流、波狀流、彈狀流及環(huán)狀流。Rouhani等［15］對兩相流的流型進(jìn)行了劃分，認(rèn)為水平流動(dòng)主要包括平滑流、波浪流、柱塞流和彌散流等。消防領(lǐng)域應(yīng)用的泡沫氣體體積分?jǐn)?shù)一般大于80%，泡沫流動(dòng)過程中是液膜包圍著氣體運(yùn)動(dòng)，泡沫形態(tài)與Rouhani等提出來的彌散流相似，泡沫在管道內(nèi)的流動(dòng)是充分發(fā)展的，且管道內(nèi)的泡沫湍流程度較高，導(dǎo)致析出的液體重新混合，從而使泡沫一直保持一個(gè)較為均勻的狀態(tài)，不會(huì)發(fā)生兩相分離現(xiàn)象［10］。故而，假設(shè)泡沫在同一管道截面為均勻流動(dòng)，進(jìn)一步將泡沫簡化為均質(zhì)單相流體［10，16］。實(shí)際工況下，對于泡沫在管路中的流動(dòng)不研究熱量的傳遞，因此模擬過程中不計(jì)流動(dòng)過程中的能量損失。假設(shè)不存在壁面滑移現(xiàn)象。

1.3 模型的選擇及邊界條件的設(shè)定

湍流在實(shí)際工業(yè)中的存在較為普遍，實(shí)際工程應(yīng)用中，泡沫流速較大，大部分的泡沫流體都處于湍流狀態(tài)。Spalart-Allmaras（S-A）模型適用于雷諾數(shù)較低的湍流情況，因此采用S-A單方程模型模擬壓縮空氣泡沫在管道內(nèi)的流動(dòng)。此外網(wǎng)格需滿足黏性影響區(qū)流體流動(dòng)求解要求，故壁面處網(wǎng)格需要局部加密［10］。選用SIMPLE算法，針對文獻(xiàn)［7］中1%的A類泡沫氣液比7∶1和10∶1的兩種發(fā)泡倍數(shù)，在不同流量下的圓形管路中的流動(dòng)進(jìn)行模擬研究。本工作采用體積流量換算出的速度作為入口邊界條件，氣液比為7∶1時(shí)入口流量分別為846，1 660，2 420，3 300，3 930，5 060 L/min，氣液比為10∶1時(shí)入口流量分別為1 240，2 230，3 300，4 420 L/min。出口邊界條件為壓力出口，由于文獻(xiàn)［7］中管路出口有球形閥，球形閥出口邊界條件的壓力值為標(biāo)準(zhǔn)大氣壓。壁面邊界條件為壁面無滑移條件。

1.4 控制方程

本工作不考慮熱量的損失，不考慮能量方程。采用S-A單方程模型求解泡沫在管路中的流動(dòng)問題。基本控制方程包括連續(xù)性方程、動(dòng)量方程以及S-A方程。

連續(xù)性方程見式（1）:

式中，ρ為流體密度，kg/m3；t為時(shí)間，s；Vx為x方向的流體流動(dòng)速度，m/s；Vy為y方向的流體流動(dòng)速度，m/s。

動(dòng)量方程見式（2）:

式中，p為流體流動(dòng)中靜壓，Pa；ui和uj為速度張量，m/s；τij為應(yīng)力張量，Pa；ρgi與Fi分別是i方向上重力體積力和外部體積力，N；xi和xj為坐標(biāo)張量。

上式中應(yīng)力張量形式可表達(dá)為式（3）:

式中，τij是由于流體中分子黏性而在微元體表面產(chǎn)生的黏性應(yīng)力τ的分量；μa為流體黏性系數(shù)；δij為克羅內(nèi)克δ符號。

圖3 氣液比為7∶1時(shí)文獻(xiàn)［7］實(shí)驗(yàn)值與模擬值結(jié)果的對比Fig.3 The comparison of experimental value in literature［7］ and simulated result when the gas-liquid ratio is 7∶1.

1.5 Fluent數(shù)值模擬

根據(jù)文獻(xiàn)［7］中的實(shí)驗(yàn)數(shù)據(jù)，需要得出模擬過程中第1個(gè)壓力傳感器與第6個(gè)壓力傳感器處的壓力值。利用Fluent軟件在傳感器1和傳感器6處建立截面，提取出該位置的壓力值。利用數(shù)據(jù)點(diǎn)1與數(shù)據(jù)點(diǎn)6的壓力值的差值作為阻力損失，并分別將數(shù)據(jù)點(diǎn)1和數(shù)據(jù)點(diǎn)6的數(shù)值模擬值與實(shí)驗(yàn)值進(jìn)行比較。

2 結(jié)果與討論

2.1 壓縮空氣泡沫壓力損失模擬分析

將文獻(xiàn)［7］中阻力損失數(shù)據(jù)與模擬結(jié)果進(jìn)行對比。氣液比分別為7∶1和10∶1時(shí)的壓力對比及阻力損失誤差見圖3和圖4，泡沫流量分別為3 000 L/min和5 000 L/min時(shí)的壓力對比及阻力損失誤差見圖5和圖6。由圖3～6可看出，模擬結(jié)果的相對誤差基本在10%以內(nèi)，數(shù)值模擬的結(jié)果在工程誤差允許的范圍內(nèi)。

圖4 氣液比為10∶1時(shí)文獻(xiàn)［7］實(shí)驗(yàn)值與模擬值結(jié)果的對比Fig.4 The comparison of experimental value in literature［7］ and simulated result when the gas-liquid ratio is 10∶1.

圖5 流量為3 000 L/min時(shí)不同發(fā)泡倍數(shù)下文獻(xiàn)［7］實(shí)驗(yàn)值與模擬值結(jié)果的對比Fig.5 The comparison of experimental value in literature［7］ and simulated result at different foaming multiples with flow rate of 3 000 L/min.

圖6 流量為5 000 L/min時(shí)不同發(fā)泡倍數(shù)下文獻(xiàn)［7］實(shí)驗(yàn)值與模擬值結(jié)果對比Fig.6 The comparison of experimental value in literature［7］ and simulated result at different foaming multiples with flow rate of 5 000 L/min.

2.2 數(shù)值模擬過程的黏度分析

由文獻(xiàn)［7］中表觀黏度隨剪切速率的變化關(guān)系可以得知，由于實(shí)驗(yàn)條件的限制只能測出剪切速率230 s-1以下的泡沫黏度，然而實(shí)際工程中泡沫流體速度較大，泡沫流體的剪切速率也較大。通過計(jì)算可以得出泡沫流體在較高流速下的剪切速率，具體方法如下。

非牛頓流體在管內(nèi)流動(dòng)的基本假設(shè):流體已充分發(fā)展，管壁處不發(fā)生滑移；不考慮流變性與時(shí)間的關(guān)系；流體在管內(nèi)等溫流動(dòng)，且流動(dòng)邊界直徑不變。根據(jù)以上假設(shè)，有非牛頓流體管內(nèi)流動(dòng)公式，即Rabinowitsch-Mooney方程，見式（5）:

式中，γw為管壁處的剪切速率，s-1；τw為管壁處的剪切應(yīng)力，Pa；u為流體的流動(dòng)速度，m/s；D為管路內(nèi)徑，m；n為流變指數(shù)。

由Rabinowitsch-Mooney方程可知，只要確定了流變指數(shù)n的值，就可以把流體速度和剪切速率聯(lián)系起來。將文獻(xiàn)［7］中的流量轉(zhuǎn)換成速度u，通過對文獻(xiàn)［7］中黏度和剪切速率的關(guān)系進(jìn)行擬合得出流變指數(shù)n，由式（5）即可以計(jì)算出實(shí)驗(yàn)條件下的剪切速率。根據(jù)陳旸等［7］的研究，壓縮空氣A類泡沫的流變規(guī)律符合Herschel-Bulkley模型，見式（6）:

公式兩邊同時(shí)除以剪切速率γ，得式（7）:

式中，K為流變系數(shù)；τ為剪切應(yīng)力，Pa；τ0為剪切速率為零時(shí)的剪切應(yīng)力，Pa；γ為剪切速率，s-1。從而計(jì)算得出高剪切速率下的泡沫流體表觀黏度。

泡沫流體是非牛頓流體，具有剪切變稀的性質(zhì)。通過對不同氣液比及不同流量下的泡沫在直徑不變的管路中流動(dòng)進(jìn)行數(shù)值模擬，最終誤差在10%之后再確定該條件下的黏度，從而得出模擬條件下不同氣液比及不同流量下壓縮空氣泡沫的黏度。將計(jì)算所得泡沫流體的表觀黏度與模擬條件得出的表觀黏度進(jìn)行對比，結(jié)果見圖7。

圖7 計(jì)算值與模擬值的對比Fig.7 Comparison of calculated and simulated values.

由圖7可看出，模擬與計(jì)算所得表觀黏度值吻合相對較好。因此Rabinowitsch-Mooney方程能夠較準(zhǔn)確地預(yù)測高流速下的泡沫剪切速率。

2.3 沿程阻力系數(shù)分析

根據(jù)1.3節(jié)中的邊界條件，對泡沫在管路中的流動(dòng)進(jìn)行數(shù)值模擬。水平管路總的壓力損失Δp可以表示為式（8）:

式中，Δpf為沿程阻力損失，Pa；Δpa為流體流動(dòng)時(shí)加速度壓力降，Pa。

當(dāng)流體流動(dòng)是充分發(fā)展的穩(wěn)態(tài)流動(dòng)時(shí)，流體的加速度壓力降是可以忽略不計(jì)的。將式（8）與達(dá)西公式結(jié)合可以得出沿程阻力系數(shù)與壓力損失的關(guān)系，見式（9）:

式中，λ為沿程阻力系數(shù)；L為管路長度，m。

從式（9）可知，沿程阻力系數(shù)、管道長徑比以及泡沫密度、流速都是影響流動(dòng)壓降的重要因素。Sun等［6］總結(jié)了沿程阻力系數(shù)的關(guān)系式，見式（10）:

式中，T為溫度，℃。

在不考慮溫度的情況下，沿程阻力系數(shù)可表示為雷諾數(shù)的函數(shù)，即λ=f（Re），如圖8所示。

對計(jì)算結(jié)果進(jìn)行擬合，得到式（11）。

對于直徑為80 mm的管路，氣液比為7∶1和10∶1的不同發(fā)泡倍數(shù)的泡沫沿程阻力系數(shù)與雷諾數(shù)的關(guān)系，它們的理論計(jì)算值與數(shù)值模擬值具有較好的一致性。林全生等［10］對泡沫流量不變、發(fā)泡倍數(shù)不變，改變管路直徑的情況進(jìn)行了研究，得出了沿程阻力系數(shù)與廣義雷諾數(shù)的關(guān)系，但是沒有考慮發(fā)泡倍數(shù)與流量的影響，因此僅適用于發(fā)泡倍數(shù)為6.7倍的情況，對于泡沫發(fā)泡倍數(shù)變化的情況則不適用。根據(jù)文獻(xiàn)［7］中的數(shù)據(jù)，對不同發(fā)泡倍數(shù)的泡沫在直徑不變的管路中的流動(dòng)進(jìn)行了研究，所以得出的沿程阻力系數(shù)與廣義雷諾數(shù)的關(guān)系僅適用于管徑不變、泡沫發(fā)泡倍數(shù)與泡沫流量變化的情況。

圖8 沿程阻力系數(shù)與雷諾數(shù)的關(guān)系Fig.8 The relation between drag coefficient along the path(λ) and Re.

2.4 不同流量及不同發(fā)泡倍數(shù)下出口壓力的變化

氣液比為7∶1和10∶1時(shí)泡沫出口壓力隨流量的變化如圖9所示。

圖9 不同氣液比下出口壓力隨流量的變化Fig.9 The variation relationship between outlet pressure and flow rate under different gas-liquid ratios.

從圖9可發(fā)現(xiàn)，在實(shí)驗(yàn)所選用的流量下，出口壓力隨流量的增加而增大，這是因?yàn)榕菽髁吭黾樱菽黧w對出口的沖擊力增大，泡沫出口壓力增大。泡沫出口壓力與測試點(diǎn)處的壓力具有相同的變化趨勢。這對實(shí)際工程應(yīng)用具有一定的指導(dǎo)意義。

3 結(jié)論

1）數(shù)值模擬所得管路中的壓力損失與實(shí)驗(yàn)值具有較好的一致性，模擬結(jié)果誤差在10%以內(nèi)。

2）根據(jù)Rabinowitsch-Mooney方程，可以計(jì)算出高流速下的剪切速率；再根據(jù)黏度隨剪切速率的變化關(guān)系，計(jì)算得出高流速下的泡沫黏度。根據(jù)Rabinowitsch-Mooney方程計(jì)算出的泡沫黏度與數(shù)值模擬所得泡沫黏度吻合較好。

3）利用達(dá)西公式分析并建立了沿程阻力系數(shù)與雷諾數(shù)的關(guān)系式λ= 2.998 3/Re0.5909，氣液比7∶1和10∶1的不同發(fā)泡倍數(shù)的泡沫沿程阻力系數(shù)與雷諾數(shù)的關(guān)系的文獻(xiàn)實(shí)驗(yàn)值與數(shù)值模擬值具有較好的一致性。

4）模擬結(jié)果表明，管路出口處的壓力隨流量的增加而增加，并且在不同測壓點(diǎn)具有相同的變化趨勢，這對實(shí)際工程應(yīng)用具有一定的指導(dǎo)意義。

符 號 說 明

CS-A模型中的常數(shù)，C=0.622

D管路內(nèi)徑，m

Fi外部體積力，N

G湍流黏性產(chǎn)生項(xiàng)，kg/（m3·s）

K流變系數(shù)

L管路長度，m

n流變指數(shù)

p流體流動(dòng)中靜壓，Pa

Δp水平管路總的壓力損失，Pa

Δpa加速度壓力降，Pa

Δpf沿程壓力降，Pa

Re雷諾數(shù)

T溫度，℃

t時(shí)間，s

u流體的流動(dòng)速度，m/s

ui速度張量，m/s

uj速度張量，m/s

vxx方向的流體流動(dòng)速度，m/s

vyy方向的流體流動(dòng)速度，m/s

xi坐標(biāo)張量

xj坐標(biāo)張量

Y壁面阻擋與黏性阻尼引起的湍流黏性的減少項(xiàng)，kg/（m3·s）

γ剪切速率，s-1

γw管壁處的剪切速率，s-1

δij克羅內(nèi)克δ符號

λ沿程阻力系數(shù)

μ流體黏度，Pa·s

μa流體黏性系數(shù)

ρ流體密度，kg/m3

ρgi重力體積力，N

σS-A湍流模型中的常數(shù)，σ=2/3

τ剪切應(yīng)力，Pa

τ0剪切速率為零時(shí)的剪切應(yīng)力，Pa

τw管壁處的剪切應(yīng)力，Pa

τij應(yīng)力張量，Pa

參 考 文 獻(xiàn)

［1］ 陸強(qiáng)，包志明，陳濤，等. A類泡沫滅火劑滅火性能試驗(yàn)研究［J］. 消防科學(xué)與技術(shù)，2013，32（2）:177-179.

［2］ 初迎霞，王海明，喬啟宇，等. 壓縮空氣泡沫系統(tǒng)在林火撲救中的應(yīng)用［J］. 林業(yè)機(jī)械與木工設(shè)備，2004，32（8）:33-36.

［3］ 傅學(xué)成，葉宏烈，包志明，等. A類泡沫滅火劑的發(fā)展與瞻望［J］. 消防科學(xué)與技術(shù)，2008，27（8）:590-592.

［4］ 傅學(xué)成，包志明，陳濤，等. 壓縮空氣泡沫的隔熱防護(hù)性能研究［J］. 消防科學(xué)與技術(shù)，2009，28（3）:204-207.

［5］ 葛曉霞，高健，趙昊. 壓縮空氣泡沫管內(nèi)流動(dòng)特性分析［J］.消防科學(xué)與技術(shù)，2016，35（10）:1408-1410.

［6］ Sun Xiao，Liang Xiaobing，Wang Shuzhong，et al. Experimental study on the rheology of CO2viscoelastic surfactant foam fracturing fluid［J］. J Pet Sci Eng，2014，119（3）:104-111.

［7］ 陳旸，陳濤，胡成，等. 壓縮空氣泡沫傳輸阻力損失研究［J］. 消防科學(xué)與技術(shù)，2017，36（10）:1418-1420.

［8］ Hu Dongfang，Huang Zhengliang，Sun Jingyuan，et al. Numerical simulation of gas-liquid flow through a 90° duct bend with a gradual contraction pipe［J］. J Zhejiang Univ Sci A，2017，18（3）:212-224.

［9］ Supa-Amornkul S，Steward F R，Lister D H. Modeling twophase flow in pipe bends［J］. J Pressure Vessel Technol，2005，127（2）:204-209.

［10］ 林全生，張猛，宋波，等. 壓縮空氣A類泡沫在水平管道內(nèi)流動(dòng)研究［J］. 消防科學(xué)與技術(shù)，2017，36（9）:1265-1268.

［11］ 張仂，谷芳. 基于計(jì)算流體力學(xué)熱流固耦合仿真的換熱器折流板結(jié)構(gòu)優(yōu)化［J］. 石油化工，2012，41（6）:682-687.

［12］ Pravakar M，Khakhar D V. Hydraulic resistance of rigid polyurethane foams. I. Effect of different surfactants on foam structure and properties［J］. J Appl Polym Sci，2004，93（6）:2821-2829.

［13］ 甘德清，高鋒，陳超，等. 管道輸送高濃度全尾砂充填料漿的阻力損失研究［J］. 礦業(yè)研究與開發(fā)，2016，36（1）:94-98.

［14］ 王新民，丁德強(qiáng)，吳亞斌，等. 膏體充填管道輸送數(shù)值模擬與分析［J］. 中國礦業(yè)，2006，15（7）:57-59.

［15］ Rouhani S Z，Sohal M S. Two-phase flow patterns:A review of research results［J］. Prog Nucl Energy，1983，11（3）:219-259.

［16］ 王勇凱，高紅，宋波，等. 壓縮空氣 A 類泡沫水平管路壓降試驗(yàn)及數(shù)值模擬［J］. 化工學(xué)報(bào)，2018，69（10）:4184-4193.