Application of MacCormack algorithm in transient flow simulation of long-distance pipeline
-
摘要: 为了提高长输管道仿真精度和速率, 基于MacCormack格式基本原理和流动方程建立了长输液体管道水力瞬变流动仿真模型。以某管道为例, 分别采用该仿真算法和特征线法模拟了阀门关闭和流量增加引起的瞬变流动, 并探讨了MacCormack格式中不同边界条件处理方法及Courant常数值对模拟结果的影响。结果表明: MacCormack格式可用于精确仿真长输管道瞬变流动过程, 采用特征线法求解预估层边界参数后, 其振荡幅度、收敛速度均优于特征线仿真方法; 同时, 采用特征线法处理MacCormack格式边界条件较线性外插法更易收敛。研究成果可为长输管道瞬变流动仿真提供参考。
-
关键词:
- 长输管道 /
- 瞬变流 /
- 模拟 /
- MacCormack /
- 边界条件
Abstract: To improve the simulation precision and speed of long-distance pipelines, the simulation model of hydraulic transient flow in long-distance liquid pipelines was established based on the basic principle of MacCormack format and flow equation.Take one pipeline as an example, in which the transient flow caused by valve closing and flow increasing was simulated by the simulation algorithm and the characteristic method respectively.The processing method for different boundary conditions in MacCormack format and the influence of Courant constant on the simulation results were discussed.The results show that, the MacCormack format can be used to accurately simulate the transient flow process of long-distance pipelines and, after the characteristic method was used to solve the estimated layer boundary parameters, its oscillation amplitude and convergence rate are better than that obtained by the characteristic simulation method.Moreover, the boundary condition of MacCormack format processed by characteristic method is conducive to converge than by linear extrapolation method.The research results can provide reference for transient flow simulation of long-distance pipelines.-
Keywords:
- long-distance pipeline /
- transient flow /
- simulation /
- MacCormack /
- boundary condition
-
阀门关闭、流量或压力突变会引起长输管道中流场瞬变, 求解该瞬变流动常见的数值计算方法有特征线法[1]、隐式差分法[2]等。尽管上述算法在油气长输管道流动分析中已经发挥了重要作用, 但随着中国大型油气管网系统建设进程不断推进, 积极探索更多算法将有助于进一步提高长输管道仿真精度和速率, 为流场瞬变行为分析提供多样化数值模拟方法。
MacCormack格式, 又称预估-校正算法, 1969年由MacCormack首次提出[3], 具有时间和空间的二阶精度, 其离散格式具有格式简单、便于理解与程序实现、占用电脑内存少、计算时间短等优点[4-5]。该格式可用于求解可压缩和不可压缩流体非定常流动, 先后在航空航天[4, 6-8]、武器装备[9-10]、河道水流[11-12]、数值传热[13-14]、地震地区管道应力[15-16]、气体钻井[17-18]等方面有所应用, 但在长输管道模拟方面应用较少。为此, 采用MacCormack格式推导了长输液体管道瞬变流动数值求解格式, 模拟了阀门关闭和流量突增引发的瞬变过程。与特征线法进行对比分析, 同时探讨了边界条件选用方法及CFL(Courant-Friedrichs-Lewy)系数对模拟结果的影响, 可为长输管道瞬变流动仿真提供参考。
1. 基于MacCormack算法的液体管道瞬变流动数值求解格式
1.1 MacCormack基本格式
流体流动微分方程可表示为矩阵形式[19], 其表达式为:
![](/fileYQCY/journal/article/yqcy/2020/5/href="yqcy-39-5-542-E1.png")
(1) 式中: U、F、H分别为微分方程的解向量、通量项、源项; t为时间, s; x为空间长度, m。
MacCormack格式算法(图 1), 由预估步和校正步两部分组成, 其表达式分别为:
![](/fileYQCY/journal/article/yqcy/2020/5/href="yqcy-39-5-542-E2.png")
(2) ![](/fileYQCY/journal/article/yqcy/2020/5/href="yqcy-39-5-542-E3.png")
(3) 式中: Uin、F in、Hin分别为矩阵U、F、H在第i节点、第n时刻的最终值;
由此得到MacCormack格式的最终值表达式为：
![](/fileYQCY/journal/article/yqcy/2020/5/href="yqcy-39-5-542-E4.png")
(4) 当使用上述格式求解时, 在空间处理上预估步采用向前差分、校正步采用向后差分。可见, MacCormack格式在时间和空间方向均具有二阶精度, 且该格式也可采用预估步向后差分、校正步向前差分的格式求解, 仿真精度相同[5]。
1.2 MacCormack数值求解格式
以等温输油管道为例, 建立瞬变流动模型的MacCormack格式, 其流动微分方程表达式为:
![](/fileYQCY/journal/article/yqcy/2020/5/href="yqcy-39-5-542-E5.png")
(5) 式中: ρ为密度, kg/m3; v为流速, m/s; p为压力, Pa; k为流体体积弹性系数, Pa; θ为管道倾角, rad; λ为水力摩阻因数; D为管道内直径, m; g为当地重力加速度, m/s2。
根据MacCormack格式的基本原理可建立上述偏微分方程组的MacCormack预估-校正仿真模型, 其预估步、校正步、平均值的表达式分别为:
![](/fileYQCY/journal/article/yqcy/2020/5/href="yqcy-39-5-542-E6.png")
(6) ![](/fileYQCY/journal/article/yqcy/2020/5/href="yqcy-39-5-542-E7.png")
(7) ![](/fileYQCY/journal/article/yqcy/2020/5/href="yqcy-39-5-542-E8.png")
(8) MacCormack格式空间和时间步长的选取遵循CFL稳定性准则[20], 其表达式为:
![](/fileYQCY/journal/article/yqcy/2020/5/href="yqcy-39-5-542-E9.png")
(9) 式中: C为Courant常数, C ≤ 1;a为流场中任一点的压力波速, m/s。
1.3 边界条件
边界上无法同时满足预估-校正格式中的向前和向后差分条件。MacCormack格式边界条件可由特征线法[21]和线性外插法得到[5]。在计算过程中, 采用体积流量Q代替流速v, 压头H代替压力p。以上游边界流量Q1、下游边界压头HN+1固定的一段管道(图 2, P、S、R为管道上的特征节点)为例, 采用2种方法可得到不同的边界条件数值求解格式。
(1) 采用特征线法求解得到管道入口和出口处的特征线格式, 其表达式为:
![](/fileYQCY/journal/article/yqcy/2020/5/href="yqcy-39-5-542-E10.png")
(10) 式中: H1、HS、HR、HN+1分别为入口、S点、R点及出口压头, m; Q1、QS、QR、QN+1分别为入口、S点、R点及出口体积流量, m3/h; aS、aR分别为S点、R点压力波速, m/s; λS、λR分别为S点、R点水力摩阻因数; θS、θR分别为S点、R点管道倾角, rad; A为管道截面积, m2。
(2) 采用线性外插法求解得到管道入口和出口处的线性外插格式, 其表达式为:
![](/fileYQCY/journal/article/yqcy/2020/5/href="yqcy-39-5-542-E11.png")
(11) 2. 实例分析
2.1 水击计算
以文献[22]中管道为例进行出口阀门关闭时的水击计算分析。该管道内径为0.2 m, 当量管长3.05 km, 等分3段(节点编号分别为0、1、2、3), 入口、出口压头分别为60 m、20 m。
在经典MacCormack格式中, 差分格式〔式(2)〕、特征线格式〔式(10)〕分别用于求取边界条件的预估值和最终解。考虑到预估-校正格式特征, 提出将特征线法同时运用于求解边界条件的预估值、最终解。因此, 分别采用MacCormack格式和特征线法求解管道出口阀门关闭后的流体水力瞬变过程, 得到节点1和节点3的压头随时间的变化曲线(图 3)。其中蓝线表示采用MacCormack格式、预估层边界采用向前差分格式〔式(2)〕; 红线表示采用MacCormack格式、预估层边界采用特征线法处理; 黑线是采用文献[22]中特征线法仿真结果。同时, 采用SPS软件模拟该工况, 得到不同仿真方法的计算结果(表 1)。可见, 用MacCormack格式、特征线法的计算结果与SPS仿真结果基本一致，求解结果可信度高。当边界条件的预估值采用差分方法时，MacCormack格式与特征线法仿真结果收敛速度、振幅基本一致，未表现出优越性，但采用特征线法处理预估层边界值，可大幅度提高MacCormack格式收敛速率。优化后的MacCormack格式显示出较原MacCormack格式以及特征线法明显的优越性，仿真结果振幅更小、收敛速度更快。其原因是：相较于特征线法求解，MacCormack格式在时间和空间上均能保证二阶精度；对比MacCormack格式的两种边界处理方式，虽然边界值的最终解都是由特征线法得到，但边界预估值将影响中间节点校正层求解，进而影响仿真结果；预估层边界值采用特征线法处理，有助于在构造中间节点的校正格式时采用更精确的解，加速收敛进程。因此，以下研究过程中均采用特征线法处理MacCormack预估层边界。
表 1 某管道阀门关闭后节点1~节点3压头随时间变化的不同仿真方法计算结果对比
MacCormack格式的时间和空间步长需满足CFL条件。Courant常数值(C值)变化会改变时空网格密集程度, 进而影响收敛进程, 适宜的C值需通过反复仿真试验确定[5]。采用MacCormack格式不同C值求解得到管道关阀后入口流量、出口压头随时间的变化曲线(图 4)。可见, 不同C值条件下的仿真收敛速度、振幅及稳定值差异均较小。因此, 对于该仿真实例, 可取C=1.00。
边界上的未知量采用线性外插法求解时，由于需要同时层相邻两个节点参数，因此必须提前求解该时层中间节点参数。由于瞬变流动本构方程的非线性特征，如果直接运用内部节点最终计算结果线性插值确定边界流动参数，将会直接导致数值求解结果发散。对此，将线性外插法进行改进：对给定边界，其预估值即为给定值；对待求边界，其预估值或校正值采用线性外插；以式（4）确定最终值。采用MacCormack格式且边界处理条件分别为特征线边界、线性外插边界、改进的线性外插边界，模拟节点1的流量、压头随时间的变化趋势（图 5）。可见，采用改进的线性外插边界能够保持MacCormack格式的精度优势，强化给定边界条件对仿真过程的影响，可降低线性外插的误差，大幅提高了线性外插处理边界条件的效果，但相较于采用特征线法求解边界条件的MacCormack格式，收敛速度过低。
2.2 流量变化
保证管道入口压头不变, 出口流量突增20%, 分别采用优化后的MacCormack格式(预估层边界采用特征线法求解)和特征线法求解该工况下的水力瞬变过程, 得到管道上各节点待求流量、压头随时间的变化曲线(图 6)。同时, 采用SPS软件模拟该工况, 得到不同仿真方法的计算结果(表 2)。可见, MacCormack格式和特征线法的计算结果与SPS仿真结果基本一致，从而验证了MacCormack算法的正确性。虽然运用特征线法和优化后的MacCormack格式计算所得到的各节点流量、压头稳定值一致，但后者的振幅、收敛速度明显优于前者。
表 2 采用不同仿真方法得到的某管道出口流量突增20% 后各节点处压头、流量稳定值对比
采用优化后的MacCormack格式不同C值得到管道出口流量突增20%后节点1处的流量、压头随时间的变化曲线(图 7)。可见, 降低C值对流量计算无影响, C值减小虽不能加速压头数据收敛, 但可改善仿真结果振荡情况, 促使压头曲线适度平滑。
采用MacCormack格式且边界处理条件分别为特征线边界、线性外插边界、改进的线性外插边界, 求解节点1的流量、压头随时间的变化趋势(图 8)。可见, 采用特征线法处理MacCormack格式的边界条件, 仿真效果最佳, 边界条件采用改进后的线性外插方法, 克服了未改进时的发散特性, 与特征线法处理效果逼近。
3. 结论
通过对MacCormack格式在长输管道瞬变流动仿真中的应用研究, 得到以下结论:
(1) 提出优化后的MacCormack格式(预估层边界值采用特征线法处理)和改进线性外插边界的MacCormack格式, 将其与特征线法、经典MacCormack格式(预估层边界采用差分法处理)、边界线性外插的MacCormack格式进行对比, 通过实例计算表明优化后的MacCormack格式效果最好; 改进线性外插边界的MacCormack格式能够克服改进前仿真结果的发散问题, 但收敛效果不及优化后的MacCormack格式。同时, 改进后的线性外插方法应用效果不稳定。
(2) 通过实例证明了Courant常数值对优化后的MacCormack格式不敏感。水击过程可取C=1.00;流量增加过程可取C=0.95。
(3) 研究成果可为长输管道瞬变流动仿真提供数值方法借鉴。今后可在已建立的水力瞬变格式基础上, 增加热力参数的MacCormack求解格式, 构造基于预估-校正原理的长输管道水热力耦合数值求解新方法。
-
[1] 刘恩斌, 彭善碧, 李长俊. 液体管道瞬变流理论及应用[M]. 北京: 石油工业出版社, 2017: 7-32. LIU E B, PENG S B, LI C J. Theory and application of transient flow in liquid pipeline[M]. Beijing: Petroleum Industry Press, 2017: 7-32.
[2] 孙建国, 王寿喜. 气体管网的动态仿真[J]. 油气储运, 2001, 20(8): 18-21. https://www.cnki.com.cn/Article/CJFDTOTAL-YQCY200108005.htm SUN J G, WANG S X. Transient simulations for gas pipeline networks[J]. Oil & Gas Storage and Transportation, 2001, 20(8): 8-21. https://www.cnki.com.cn/Article/CJFDTOTAL-YQCY200108005.htm
[3] MACCORMACK R W. The effect of viscosity in hypervelocity impact cratering[J]. Journal of Spacecraft and Rockets, 2003, 40(5): 757-763. doi: 10.2514/2.6901
[4] 李松波. 用Mac Cormack二步格式计算无粘流场的数值实践[J]. 空气动力学学报, 1983(4): 1-9. https://www.cnki.com.cn/Article/CJFDTOTAL-KQDX198304000.htm LI S B. Numerical practice for calculating inviscid flow field using Mac Cormack two-step scheme[J]. Acta Aerodynamica Sinica, 1983(4): 1-9. https://www.cnki.com.cn/Article/CJFDTOTAL-KQDX198304000.htm
[5] ANDERSON J R. 计算流体力学入门[M]. 姚朝辉, 周强, 译. 北京: 清华大学出版社, 2010: 175-178. ANDERSON J R. Computational fluid dynamics: the basics with applications[M]. Translated by YAN C H and ZHOU Q. Beijing: Tsinghua University Press, 2010: 175-178.
[6] HIXON R, TURKEL E. Compact implicit Mac Cormack-type schemes with high accuracy[J]. Journal of Computational Physics, 2000, 158(1): 51-70. doi: 10.1006/jcph.1999.6406
[7] SELLE A, FEDKIW R, KIM B M, LIU Y J, ROSSIGNAC J. An unconditionally stable Mac Cormack method[J]. Journal of Scientific Computing, 2008, 35(2/3): 350-371.
[8] 田辉, 蔡国飙, 孙再庸. N2O/HTPB固液火箭发动机喷管两相流计算[J]. 航空动力学报, 2008, 23(6): 1146-1150. https://www.cnki.com.cn/Article/CJFDTOTAL-HKDI200806032.htm TIAN H, CAI G B, SUN Z Y. Computation of nozzle two-phase flow in a N2O/HTPB hybrid rocket motor[J]. Journal of Aerospace Power, 2008, 23(6): 1146-1150. https://www.cnki.com.cn/Article/CJFDTOTAL-HKDI200806032.htm
[9] 陈龙淼, 钱林方. 复合材料身管膛内气流与固壁传热分析[J]. 南京理工大学学报(自然科学版), 2008, 32(3): 327-332. doi: 10.3969/j.issn.1005-9830.2008.03.015 CHEN L M, QIAN L F. Heat transfer between gas and inner surface of composite material barrels[J]. Journal of Naming University of Science and Technology(Natural Science), 2008, 32(3): 327-332. doi: 10.3969/j.issn.1005-9830.2008.03.015
[10] 李海庆, 张小兵, 李筱炜, 袁亚雄. 激光多点点火二维两相流数值模拟[J]. 兵工学报, 2012, 33(3): 257-260. https://www.cnki.com.cn/Article/CJFDTOTAL-BIGO201203002.htm LI H Q, ZHANG X B, LI X W, YUAN Y X. Two-dimensional two-phase flow numerical simulation of multi-point laser ignition[J]. Acta Armamentarii, 2012, 33(3): 257-260. https://www.cnki.com.cn/Article/CJFDTOTAL-BIGO201203002.htm
[11] AMOKRANE H, VILLENEUVE J P. A numerical method for solving the water flow equation in unsaturated porous media[J]. Ground Water, 1996, 34(4): 666-674. doi: 10.1111/j.1745-6584.1996.tb02054.x
[12] INAN A, BALAS L. Numerical modeling of extended mild slope equation with modified Mac Cormack method[J]. WSEASTransactions on Fluid Mechanics, 2009, 4(1): 14-23.
[13] GLASS D E, ÖZIŞIK M N, MCRAE D S, VICK B. Hyperbolic heat conduction with temperature-dependent thermal conductivity[J]. Journal of Applied Physics, 1986, 59(6): 1861-1865. doi: 10.1063/1.336413
[14] MALINOWSKI L, BIELSKI S. Transient temperature field in a parallel-flow three-fluid heat exchanger with the thermal capacitance of the walls and the longitudinal walls conduction[J]. Applied Thermal Engineering, 2009, 29(5): 877-883.
[15] MARTINO D G, ACUNTO B D, FONTANA N, GIUGNI M. Dynamic response of continuous buried pipes in seismic areas[C]. Chicago: Conference of Pipeline 2006: Service to the Owner, 2006: 1-9.
[16] CORRADO V, D ACUNTO B, FONTANA N, GIUGNI M. Estimation of dynamic strains in finite end-constrained pipes in seismic areas[J]. Mathematical and Computer Modelling, 2009, 49(3/4): 789-797.
[17] 唐贵, 程宏英, 孟英峰. 气体钻井过程中的瞬态流动分析[J]. 钻采工艺, 2006, 29(2): 5-6, 19. https://www.cnki.com.cn/Article/CJFDTOTAL-ZCGY200602002.htm TANG G, CHENG H Y, MENG Y F. Transient flow analysis in gas drilling[J]. Drilling & Production Technology, 2006, 29(2): 5-6, 19. https://www.cnki.com.cn/Article/CJFDTOTAL-ZCGY200602002.htm
[18] 唐贵, 孟英峰, 王存新, 舒梅. 气体钻井过程中井内不稳定流动与现场试验研究[J]. 天然气工业, 2005, 25(2): 75-77. https://www.cnki.com.cn/Article/CJFDTOTAL-TRQG20050200O.htm TANG G, MENG Y F, WANG C X, SHU M. Borehole gas unsteady flow in the process of gas drilling and its field test study[J]. Natural Gas Industry, 2005, 25(2): 75-77. https://www.cnki.com.cn/Article/CJFDTOTAL-TRQG20050200O.htm
[19] 李长俊, 贾文龙. 油气管道多相流[M]. 北京: 化学工业出版社, 2015: 188-195. LI C J, JIA W L. The multiphase flow in oil and gas pipeline[M]. Beijing: Chemical Industry Press, 2015: 188-195.
[20] 邓志彬. 天然气管道堵塞检测理论方法研究[D]. 成都: 西南石油大学, 2016: 1-139. DENG Z B. The partial blockage detection theory and method for the natural gas pipeline[D]. Chengdu: Southwest Petroleum University, 2016: 1-139.
[21] AMARA L, BERREKSI A, ACHOUR B. Adapted Mac Cormack finite-differences schemes for water hammer simulation[J]. Journal of Civil Engineering and Science, 2013, 2(4): 226-233.
[22] 蒋仕章, 蒲家宁. 水力瞬变特征线法和隐式差分法的对比分析[J]. 油气储运, 2001, 20(1): 12-14. https://www.cnki.com.cn/Article/CJFDTOTAL-YQCY200101003.htm JIANG S Z, PU J N. Correlation analysis of characteristic method and implicit difference method for fluid transient in long distance pipeline[J]. Oil & Gas Storage and Transportation, 2001, 20(1): 12-14. https://www.cnki.com.cn/Article/CJFDTOTAL-YQCY200101003.htm
下载:
计量
- 文章访问数:
- HTML全文浏览量: 0
- PDF下载量: