压缩空气泡沫管路输送的模拟计算

林全生，高 红，宋 波

（1. 天津大学 化工学院，天津 300350；2. 公安部天津消防研究所，天津 300381）

［关键字］泡沫；数值模拟；黏度；沿程阻力系数；Fluent软件

压缩空气泡沫灭火剂是一类主要用于扑救固体可燃物火灾的新型高效泡沫灭火剂，具备产生高稳定泡沫的性能，可用于石油化工可燃性液体危险品储罐火灾时的比邻储罐防护［1-4］。该类灭火剂在固定式压缩空气泡沫系统中具有突出的灭火效果。但泡沫在管道运输过程中，以水为基础的现有理论难以满足要求，需要找到适合压缩空气泡沫流体流动的计算方法，为工程设计提供准确依据。针对泡沫，利用单相流对其进行研究的较多。葛晓霞等［5］通过搭建实验平台来观察泡沫在管内的流态，给出了压缩空气泡沫流型的类型，并通过理论分析得出了模型参数的算法，确定了摩阻系数及压力损失的计算方程。Sun等［6］通过实验发现泡沫流体是剪切稀化的幂率流体，并得到了摩阻系数与广义雷诺数的关系。陈旸等［7］开展了压缩空气泡沫管网输送实验，研究了泡沫流量、气液比、泡沫种类对输送管路压力变化的影响及摩擦阻力损失。Hu等［8］对气液两相流在90°弯管中的流动情况进行了模拟，得到了一定气液比下双相流弯管压力分布规律。Supa-Amornkul等［9］分别采用单相流和双相流对弯管压降进行数值模拟，气体体积分数在50%以上的情况下，沿长度方向单相流压力分布与实验结果差别不大，相比之下，多相流并不能很好地预测弯管压降。林全生等［10］利用上海消防总队的实验数据，对泡沫在不同管径、不同长度的管路中的流动阻力损失进行了数值模拟，并得出了范宁摩擦因子与广义雷诺数的拟合关系式。Fluent数值模拟已在工程上有很多应用［11-14］，它在消防泡沫输运方面的研究还有很大的探索空间。

本工作利用Fluent数值模拟计算软件对压缩空气泡沫在管道内的流动进行了数值模拟研究，考虑了发泡倍数和流量的影响，并与相关实验数据进行了对比，通过拟合得出沿程阻力系数的显性表达式，为压缩空气泡沫的工程计算提供指导。

1 数值模拟计算方法

1.1 几何模型和网格划分

基于有限元方法，根据文献［7］，管路设备尺寸如图1所示，利用Solidworks软件建立三维几何模型。利用ANSYS Workbench16.0中的mesh模块对几何模型进行四面体网格划分，网格数量为113万，并对几何模型的壁面设置膨胀层进行边界层加密处理，如图2所示。

图1 实体管路示意图［7］Fig.1 Diagram of solid piping［7］.

图2 几何模型网格划分Fig.2 Mesh division of geometric model.

1.2 计算模型假设

葛晓霞等［5］对压缩空气泡沫进行了流型的划分，将泡沫在管内的流动形态按含气率从低到高的顺序依次划分为泡状流、塞状流、层状流、波状流、弹状流及环状流。Rouhani等［15］对两相流的流型进行了划分，认为水平流动主要包括平滑流、波浪流、柱塞流和弥散流等。消防领域应用的泡沫气体体积分数一般大于80%，泡沫流动过程中是液膜包围着气体运动，泡沫形态与Rouhani等提出来的弥散流相似，泡沫在管道内的流动是充分发展的，且管道内的泡沫湍流程度较高，导致析出的液体重新混合，从而使泡沫一直保持一个较为均匀的状态，不会发生两相分离现象［10］。故而，假设泡沫在同一管道截面为均匀流动，进一步将泡沫简化为均质单相流体［10，16］。实际工况下，对于泡沫在管路中的流动不研究热量的传递，因此模拟过程中不计流动过程中的能量损失。假设不存在壁面滑移现象。

1.3 模型的选择及边界条件的设定

湍流在实际工业中的存在较为普遍，实际工程应用中，泡沫流速较大，大部分的泡沫流体都处于湍流状态。Spalart-Allmaras（S-A）模型适用于雷诺数较低的湍流情况，因此采用S-A单方程模型模拟压缩空气泡沫在管道内的流动。此外网格需满足黏性影响区流体流动求解要求，故壁面处网格需要局部加密［10］。选用SIMPLE算法，针对文献［7］中1%的A类泡沫气液比7∶1和10∶1的两种发泡倍数，在不同流量下的圆形管路中的流动进行模拟研究。本工作采用体积流量换算出的速度作为入口边界条件，气液比为7∶1时入口流量分别为846，1 660，2 420，3 300，3 930，5 060 L/min，气液比为10∶1时入口流量分别为1 240，2 230，3 300，4 420 L/min。出口边界条件为压力出口，由于文献［7］中管路出口有球形阀，球形阀出口边界条件的压力值为标准大气压。壁面边界条件为壁面无滑移条件。

1.4 控制方程

本工作不考虑热量的损失，不考虑能量方程。采用S-A单方程模型求解泡沫在管路中的流动问题。基本控制方程包括连续性方程、动量方程以及S-A方程。

连续性方程见式（1）:

式中，ρ为流体密度，kg/m3；t为时间，s；Vx为x方向的流体流动速度，m/s；Vy为y方向的流体流动速度，m/s。

动量方程见式（2）:

式中，p为流体流动中静压，Pa；ui和uj为速度张量，m/s；τij为应力张量，Pa；ρgi与Fi分别是i方向上重力体积力和外部体积力，N；xi和xj为坐标张量。

上式中应力张量形式可表达为式（3）:

式中，τij是由于流体中分子黏性而在微元体表面产生的黏性应力τ的分量；μa为流体黏性系数；δij为克罗内克δ符号。

图3 气液比为7∶1时文献［7］实验值与模拟值结果的对比Fig.3 The comparison of experimental value in literature［7］ and simulated result when the gas-liquid ratio is 7∶1.

1.5 Fluent数值模拟

根据文献［7］中的实验数据，需要得出模拟过程中第1个压力传感器与第6个压力传感器处的压力值。利用Fluent软件在传感器1和传感器6处建立截面，提取出该位置的压力值。利用数据点1与数据点6的压力值的差值作为阻力损失，并分别将数据点1和数据点6的数值模拟值与实验值进行比较。

2 结果与讨论

2.1 压缩空气泡沫压力损失模拟分析

将文献［7］中阻力损失数据与模拟结果进行对比。气液比分别为7∶1和10∶1时的压力对比及阻力损失误差见图3和图4，泡沫流量分别为3 000 L/min和5 000 L/min时的压力对比及阻力损失误差见图5和图6。由图3～6可看出，模拟结果的相对误差基本在10%以内，数值模拟的结果在工程误差允许的范围内。

图4 气液比为10∶1时文献［7］实验值与模拟值结果的对比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时不同发泡倍数下文献［7］实验值与模拟值结果的对比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时不同发泡倍数下文献［7］实验值与模拟值结果对比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 数值模拟过程的黏度分析

由文献［7］中表观黏度随剪切速率的变化关系可以得知，由于实验条件的限制只能测出剪切速率230 s-1以下的泡沫黏度，然而实际工程中泡沫流体速度较大，泡沫流体的剪切速率也较大。通过计算可以得出泡沫流体在较高流速下的剪切速率，具体方法如下。

非牛顿流体在管内流动的基本假设:流体已充分发展，管壁处不发生滑移；不考虑流变性与时间的关系；流体在管内等温流动，且流动边界直径不变。根据以上假设，有非牛顿流体管内流动公式，即Rabinowitsch-Mooney方程，见式（5）:

式中，γw为管壁处的剪切速率，s-1；τw为管壁处的剪切应力，Pa；u为流体的流动速度，m/s；D为管路内径，m；n为流变指数。

由Rabinowitsch-Mooney方程可知，只要确定了流变指数n的值，就可以把流体速度和剪切速率联系起来。将文献［7］中的流量转换成速度u，通过对文献［7］中黏度和剪切速率的关系进行拟合得出流变指数n，由式（5）即可以计算出实验条件下的剪切速率。根据陈旸等［7］的研究，压缩空气A类泡沫的流变规律符合Herschel-Bulkley模型，见式（6）:

公式两边同时除以剪切速率γ，得式（7）:

式中，K为流变系数；τ为剪切应力，Pa；τ0为剪切速率为零时的剪切应力，Pa；γ为剪切速率，s-1。从而计算得出高剪切速率下的泡沫流体表观黏度。

泡沫流体是非牛顿流体，具有剪切变稀的性质。通过对不同气液比及不同流量下的泡沫在直径不变的管路中流动进行数值模拟，最终误差在10%之后再确定该条件下的黏度，从而得出模拟条件下不同气液比及不同流量下压缩空气泡沫的黏度。将计算所得泡沫流体的表观黏度与模拟条件得出的表观黏度进行对比，结果见图7。

图7 计算值与模拟值的对比Fig.7 Comparison of calculated and simulated values.

由图7可看出，模拟与计算所得表观黏度值吻合相对较好。因此Rabinowitsch-Mooney方程能够较准确地预测高流速下的泡沫剪切速率。

2.3 沿程阻力系数分析

根据1.3节中的边界条件，对泡沫在管路中的流动进行数值模拟。水平管路总的压力损失Δp可以表示为式（8）:

式中，Δpf为沿程阻力损失，Pa；Δpa为流体流动时加速度压力降，Pa。

当流体流动是充分发展的稳态流动时，流体的加速度压力降是可以忽略不计的。将式（8）与达西公式结合可以得出沿程阻力系数与压力损失的关系，见式（9）:

式中，λ为沿程阻力系数；L为管路长度，m。

从式（9）可知，沿程阻力系数、管道长径比以及泡沫密度、流速都是影响流动压降的重要因素。Sun等［6］总结了沿程阻力系数的关系式，见式（10）:

式中，T为温度，℃。

在不考虑温度的情况下，沿程阻力系数可表示为雷诺数的函数，即λ=f（Re），如图8所示。

对计算结果进行拟合，得到式（11）。

对于直径为80 mm的管路，气液比为7∶1和10∶1的不同发泡倍数的泡沫沿程阻力系数与雷诺数的关系，它们的理论计算值与数值模拟值具有较好的一致性。林全生等［10］对泡沫流量不变、发泡倍数不变，改变管路直径的情况进行了研究，得出了沿程阻力系数与广义雷诺数的关系，但是没有考虑发泡倍数与流量的影响，因此仅适用于发泡倍数为6.7倍的情况，对于泡沫发泡倍数变化的情况则不适用。根据文献［7］中的数据，对不同发泡倍数的泡沫在直径不变的管路中的流动进行了研究，所以得出的沿程阻力系数与广义雷诺数的关系仅适用于管径不变、泡沫发泡倍数与泡沫流量变化的情况。

图8 沿程阻力系数与雷诺数的关系Fig.8 The relation between drag coefficient along the path(λ) and Re.

2.4 不同流量及不同发泡倍数下出口压力的变化

气液比为7∶1和10∶1时泡沫出口压力随流量的变化如图9所示。

图9 不同气液比下出口压力随流量的变化Fig.9 The variation relationship between outlet pressure and flow rate under different gas-liquid ratios.

从图9可发现，在实验所选用的流量下，出口压力随流量的增加而增大，这是因为泡沫流量增加，泡沫流体对出口的冲击力增大，泡沫出口压力增大。泡沫出口压力与测试点处的压力具有相同的变化趋势。这对实际工程应用具有一定的指导意义。

3 结论

1）数值模拟所得管路中的压力损失与实验值具有较好的一致性，模拟结果误差在10%以内。

2）根据Rabinowitsch-Mooney方程，可以计算出高流速下的剪切速率；再根据黏度随剪切速率的变化关系，计算得出高流速下的泡沫黏度。根据Rabinowitsch-Mooney方程计算出的泡沫黏度与数值模拟所得泡沫黏度吻合较好。

3）利用达西公式分析并建立了沿程阻力系数与雷诺数的关系式λ= 2.998 3/Re0.5909，气液比7∶1和10∶1的不同发泡倍数的泡沫沿程阻力系数与雷诺数的关系的文献实验值与数值模拟值具有较好的一致性。

4）模拟结果表明，管路出口处的压力随流量的增加而增加，并且在不同测压点具有相同的变化趋势，这对实际工程应用具有一定的指导意义。

符 号 说 明

CS-A模型中的常数，C=0.622

D管路内径，m

Fi外部体积力，N

G湍流黏性产生项，kg/（m3·s）

K流变系数

L管路长度，m

n流变指数

p流体流动中静压，Pa

Δp水平管路总的压力损失，Pa

Δpa加速度压力降，Pa

Δpf沿程压力降，Pa

Re雷诺数

T温度，℃

t时间，s

u流体的流动速度，m/s

ui速度张量，m/s

uj速度张量，m/s

vxx方向的流体流动速度，m/s

vyy方向的流体流动速度，m/s

xi坐标张量

xj坐标张量

Y壁面阻挡与黏性阻尼引起的湍流黏性的减少项，kg/（m3·s）

γ剪切速率，s-1

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

δij克罗内克δ符号

λ沿程阻力系数

μ流体黏度，Pa·s

μa流体黏性系数

ρ流体密度，kg/m3

ρgi重力体积力，N

σS-A湍流模型中的常数，σ=2/3

τ剪切应力，Pa

τ0剪切速率为零时的剪切应力，Pa

τw管壁处的剪切应力，Pa

τij应力张量，Pa

参 考 文 献

［1］ 陆强，包志明，陈涛，等. A类泡沫灭火剂灭火性能试验研究［J］. 消防科学与技术，2013，32（2）:177-179.

［2］ 初迎霞，王海明，乔启宇，等. 压缩空气泡沫系统在林火扑救中的应用［J］. 林业机械与木工设备，2004，32（8）:33-36.

［3］ 傅学成，叶宏烈，包志明，等. A类泡沫灭火剂的发展与瞻望［J］. 消防科学与技术，2008，27（8）:590-592.

［4］ 傅学成，包志明，陈涛，等. 压缩空气泡沫的隔热防护性能研究［J］. 消防科学与技术，2009，28（3）:204-207.

［5］ 葛晓霞，高健，赵昊. 压缩空气泡沫管内流动特性分析［J］.消防科学与技术，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］. 消防科学与技术，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类泡沫在水平管道内流动研究［J］. 消防科学与技术，2017，36（9）:1265-1268.

［11］ 张仂，谷芳. 基于计算流体力学热流固耦合仿真的换热器折流板结构优化［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］. 矿业研究与开发，2016，36（1）:94-98.

［14］ 王新民，丁德强，吴亚斌，等. 膏体充填管道输送数值模拟与分析［J］. 中国矿业，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 类泡沫水平管路压降试验及数值模拟［J］. 化工学报，2018，69（10）:4184-4193.