基于查表插值算法的高压CO...2管道泄漏瞬态特性数值模拟_鲁寨军.pdf

文章编号:1009-6094(2023)02-0442-09基于查表插值算法的高压 CO2管道泄漏瞬态特性数值模拟*鲁寨军1，2，3，王甲强1，2，3，张梓轩1，2，3，孙永龙1，2，3，王根达1，2，3，刘东润1，2，3(1 中南大学轨道交通安全教育部重点实验室，长沙 410075;2 轨道交通安全关键技术国际合作联合实验室，长沙 410075;3 轨道交通列车安全保障技术国家地方联合工程研究中心，长沙 410075)摘要:为准确高效地模拟高压 CO2管道泄漏的瞬态特性，基于 Fluent 仿真平台，利用用户自定义真实气体模型(UserDefined Real Gas Model，UDRGM)和用户自定义函数(UserDefined Function，UDF)，结合查表法和双线性插值法建立CO2的真实气体模型，并将压力驱动的 Lee 模型通过用户自定义函数嵌入 Fluent 求解器来模拟 CO2的非平衡相变过程，建立了高压 CO2管道泄漏的非平衡相变数值模型。通过与Botros 等的试验数据进行对比分析，验证了该模型的准确性。在此基础上，对比了上述模拟方法与编译 SW(Span-Wagner)状态方程模拟方法的精度和效率，最后使用本模型研究了不同初始压力对高压 CO2管道泄漏瞬态特性的影响。结果表明:两种模拟方法精度接近，最大相差为 7.37%，但提出的模拟方法效率明显优于编译 SW 状态方程的模拟方法，计算时间相较缩短约 86.9%;初始压力为 11.27 MPa 的最大总出口质量流量比 4.36 MPa 的大 7.24 kg/s，而对应的气相 CO2出口质量流量却低 0.14 kg/s;初始压力为 11.27MPa 和 4.36 MPa 对应的近场最大射流流速分别为 155.12m/s 和 175.50 m/s;管内初始压力越低，泄漏后到达亚稳态时过热度越大，管内及泄漏口附近相变程度越剧烈，气相出口质量流量越大，近场射流峰值速度越大。这项研究可为工业规模管道泄漏的模拟和管道泄漏的三维模拟提供更高效的方式。关键词:安全工程;环境工程;二氧化碳;非平衡相变;计算流体力学;管道泄漏;气液两相流;数值模拟中图分类号:X937文献标志码:ADOI:10.13637/j issn 1009-6094.2021.1814*收稿日期:20211013作者简介:鲁寨军，教授，博士，博士生导师，从事高速列车碰撞安全保护及二氧化碳相变特性研究，qlzjzd csu edu cn;刘东润(通信作者)，助理研究员，博士，从事高速列车空气动力学研究，csuldr foxmail com。0引言碳捕捉与封存(Carbon Capture and Storage，CCS)是实现碳达峰和碳中和最直接的技术方式，而管道运输是 CCS 技术中输运二氧化碳(CO2)的最主要方式1。为了提高输运效率，CO2通常以高压形式输送，过程中难免发生泄漏事故。而泄漏时 CO2流动的瞬态特性对事故的预防、发展及造成后果预测至关重要。因此，开展高压 CO2管道泄漏特性研究，对 CO2高压输运管道设计及管道泄漏风险预测评估具有重要意义。对于 CO2管道泄漏特性的研究主要分为试验研究和数值仿真研究。但由于试验研究设计和实施条件要求高、难度大、成本高，数值仿真研究成为近年来的研究热点。目前国内外学者已经进行了大量的数值仿真研究并取得了丰硕的成果，其中有大量的研究工作围绕在气体物性模型和相变模型等决定数值仿真精度的两个方面。2003 年，Oke 等2 首次基于均匀平衡模型(Homogeneous Equilibrium Model，HEM)建立了一个高压管道穿孔数值模型，通过对比发现，在预测泄漏质量流量方面数值计算结果与试验结果吻合良好。2017 年，Zheng 等3 建立了CO2减压过程的时间松弛相变模型(HomogenousRelaxation Model，HRM)，成功预测了 CO2管内减压过程中 的 液气 相 变 过 程 和 气固 相 变 过 程。Giacomelli 等4 基于 Mixture 多相流模型对 CO2管道泄漏进行大量仿真并对比了 HEM 和 Mixture 多相流模型的相关精度。以上研究表明，HEM 认为气液两相之间存在热力学平衡和力学平衡，由于其忽略了相滑移和非平衡相变，预测精度并不高;而 HRM 中引入了松弛系数来量化非平衡现象对相变质量源相的影响，其模拟结果具有更高的精度;Mixture 多相流模型可以考虑气液相滑移及松弛系数来量化非平衡现象对相变质量源相的影响，同样具有较高的仿真精度。Flechas 等5 对 比 了 分 别 采 用 PR(Peng-Robinson)和 SW(Span-Wagner)方程作为真实气体状态方程的数值模型在 CO2管道泄漏质量流量预测精度上的差异。2019 年，刘斌等6 使用精确的状态 方 程 PR、GERG(GroupeEuropeendeRecherches Gazieres)建立了一维计算流体动力学管道破裂减压预测模型，来实现更好的源强度估计。Xiao 等7 8、鲁寨军等9 利用 Fluent 中户自定义真实气体模型和用户自定义函数将 SW 状态方程嵌入 Fluent 求解器中建立了管道全尺寸泄漏模型。以上研究发现，真实气体状态方程明显优于理想气体状态方程，SW 状态方程和 GERG2008 状态方程244第 23 卷第 2 期2023 年 2 月安全 与 环 境 学 报Journal of Safety and EnvironmentVol 23No 2Feb，2023较 Peng-Robinson 状态方程具有更高的精度，但 SW 状态方程在纯 CO2的计算中更加优异。由于 SW 状态方程的复杂性，且其为密度和温度的显式函数，使用其进行热力学参数计算必须依赖于给定的密度和温度。而 Fluent 平台中的用户自定义真实气体模型中的热力学性质方程要求为压力和温度的显式方程，因此在编译 SW 状态方程时必须将其显式化，编译过程工作量大、难度高、收敛性差。并且为确保显式化过程热力学性质的计算精度，其在计算过程中需要设置很小时间步长，CPU 计算时间很长。因此对于工业规模管道泄漏的模拟和管道泄漏的三维模拟有较大的挑战。因此，本文基于 Fluent 中用户自定义真实气体模型和用户自定义函数，结合 NIST REFPROP(美国国家标准技术研究所开发的工质物性计算软件https:/webbook nist gov/chemistry/fluid)中 CO2数据，通过查表法结合双线性插值法将 CO2液相、气相、超临界相及亚稳态相的热力学性质载入 Fluent中。并基于 Mixture 多相流模型和 Lee 相变模型，建立 CO2管道泄漏的二维非平衡相变数值模型。首先与 Botros 等10 的试验数据进行对比分析，验证该模型的准确性，其次从计算精度和计算时间两个方面与编译 SW 状态方程的模拟方法进行对比，最后使用该模型分别研究管道初始压力为 11.27 MPa 和4.36 MPa 的液态 CO2的泄漏特性，对比不同压力对泄漏过程中管内、泄漏口和管外近场泄漏特性的影响。1仿真模型和方法1.1多相流模型本文中使用的数值模型主要基于混合模型，该模型适用于非均匀多相流的模拟。连续性方程可以表示如下。mt+(mum)=0(1)式中m和 um分别为混合密度和混合质量平均速度，t 表示时间。混合密度 m由下式计算得出。m=3i=1ii(2)式中i、i分别是第 i 相的体积分数和密度，本文中混合物分为 3 相 CO2气相、CO2液相和空气相，分别用下标 1、2、3 来表示。因此，1+2+3=1(3)混合相的质量平均速度 um定义如下。um=3i=1iiuim(4)式中ui是第 i 相速度。混合相的动量方程采用以下形式。(mum)t+(mumum)=p+fe(5)式中p 是压力，fe是有效动量通量，可以通过以下公式计算。fe=m(um+uTm)+mg+F+3i=1iiud，iud，()i(6)式中m是混合物黏度，g 是重力加速度，F 为外部体积力，T 为温度，ud，i表示第 i 相的漂移速度(driftvelocity)，定义为ud，i=ui um(7)混合相的能量守恒方程采用以下形式。t3i=1(iiEi)+3i=1(iui(iEi+p)=(keT+SE)(8)式中Ei是第 i 相的能量，ke是有效热导率，SE是体积源项。在混合模型中，考虑了流体的可压缩性。式(8)中的 Ei可以表示如下。Ei=hi+u2i2pi(9)其中hi是第 i 相的显焓。由连续性方程导出的气体体积分数与传质源之间的关系如下。(11)t+(11u1)=(11ud，1)+Sm(10)式中1、1、u1分别是 CO2气相的体积分数、密度和速度;ud，1是 CO2气相的漂移速度;Sm是气体的传质源项，Sm与 CO2的蒸发和冷凝过程有关。1.2非平衡相变模型当高压 CO2从管道泄漏时，过程中可能会发生气 液相变(蒸发)或液 气相变(冷凝)，Lee 模型是较为常见的模型。原始的 Lee 模型11 是以饱和温度为判断相变条件的，但由于 CO2管道泄漏过程是压力驱动的液气非平衡相变，以饱和压力作为相变发生的判断依据更准确7。即当 p psat时，发生气 液相变。m12=1211psat p1psat(11)3442023 年 2 月鲁寨军，等:基于查表插值算法的高压 CO2管道泄漏瞬态特性数值模拟Feb，2023当 p psat时，发生液 气相变。m21=2122psat p1psat(12)式中下标 sat 表示饱和状态参数;m21是从液相到气体的传质量，m12是从气相到液相的传质量;21和12分别是蒸发因子和冷凝因子，该系数用于控制相变传质过程的快慢。1.3CO2热力学性质模型由于高压 CO2管道泄漏过程中涉及多态相变和传质传热，为了能够更加精确地模拟其泄漏的瞬态特性，必须使用精确的真实气体方程来描述 CO2的热力学性质。而 Fluent 平台内置的真实气体方程如PR 状态方程描述密相和部分液相的准确性较差，内置的 NIST real gas models 定义的真实气体模型无法适用于多相流模型12。因此大部分学者在 Fluent平台中采用用户自定义真实气体模型来定义可压缩流体的热物性质，但其在多相流的仿真中只能使用一个用户自定义真实气体模型12，如气相或液相。因此本文采用用户自定义真实气体模型和用户自定义函数来分别建立 CO2液相和气相的热力学性质模型，见图 1，其中 为密度，h 为焓，s 为熵，c 为声速，cp为比定压热容，M 为相对分子质量，为黏度，kc为导热系数，d/dT 为密度 对温度 T 的导数，d/dp为密度 对压力 p 的导数，dh/dT 为焓 h 对温度 T的导数，dh/dp 为焓 h 对压力 p 的导数。由于在管道运输过程中，管内 CO2主要以高压液态的形式输送，管道破裂后发生的主要为液相向气相的相变，因此将液态设置为主相，使用用户自定义真实气体模型来定义，气相为次相，使用用户自定义函数来定义。虽然用户自定义函数中定义比定压热容和焓值的宏仅为温度的函数，而实际计算过程中流体比定压热容和焓值由压力和温度共同决定，二者存在一定误差。且用户自定义函数中忽略了密度及焓对压力和温度的梯度项，但根据 Giacomelli 等4 的研究结果，这些误差对最终的计算结果影响很小。本文预先使用 NIST REFPROP 数据库将 CO2的热物性质建立成关于压力、温度的三维数据表。Giacomelli 等4 和 Romei 等13 对查找表间隔步长进行的灵敏度分析显示，可以在不需要过度细化的情况下实现足够的精度水平。因此，本文采用压力间隔为 0.1 MPa，温度间隔为 0.5 K，共 17 个数据表。该数据表涵盖了 CO2液相、超临界相、气相及亚稳态气相和亚稳态液相。并且为了确保计算的顺利进行，采用外插法将亚稳态液和亚稳态气性质估算至各自亚稳界线(spinodal