CN113962131A - 一种高效模拟大型天然气管网流动传热的方法 - Google Patents
一种高效模拟大型天然气管网流动传热的方法 Download PDFInfo
- Publication number
- CN113962131A CN113962131A CN202111304959.0A CN202111304959A CN113962131A CN 113962131 A CN113962131 A CN 113962131A CN 202111304959 A CN202111304959 A CN 202111304959A CN 113962131 A CN113962131 A CN 113962131A
- Authority
- CN
- China
- Prior art keywords
- pipeline
- fluid
- pipe network
- grid
- equation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 title claims abstract description 56
- 238000000034 method Methods 0.000 title claims abstract description 38
- 239000003345 natural gas Substances 0.000 title claims abstract description 28
- 239000011159 matrix material Substances 0.000 claims abstract description 44
- 238000004088 simulation Methods 0.000 claims abstract description 10
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 9
- 239000012530 fluid Substances 0.000 claims description 44
- 230000008859 change Effects 0.000 claims description 16
- 230000004323 axial length Effects 0.000 claims description 8
- 230000001788 irregular Effects 0.000 claims description 7
- 230000001133 acceleration Effects 0.000 claims description 6
- 238000002939 conjugate gradient method Methods 0.000 claims description 3
- 125000004122 cyclic group Chemical group 0.000 claims description 3
- 230000005484 gravity Effects 0.000 claims description 2
- 238000004364 calculation method Methods 0.000 abstract description 7
- 238000010586 diagram Methods 0.000 description 7
- 230000008569 process Effects 0.000 description 3
- 238000010276 construction Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 239000007789 gas Substances 0.000 description 2
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 238000004134 energy conservation Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 239000003209 petroleum derivative Substances 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 230000001052 transient effect Effects 0.000 description 1
- 238000011144 upstream manufacturing Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/14—Pipes
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/08—Thermal analysis or thermal optimisation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Computing Systems (AREA)
- Fluid Mechanics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开一种高效模拟大型天然气管网流动传热的方法,涉及天然气管网数据模拟领域,本方法在使用SIMPLE算法来求解控制方程的基础上,结合离散矩阵的稀疏性与管网结构特点,利用循环三对角矩阵求解算法与开源向量库,构建了管网中管道节点分开模拟的方法,从而加快了计算速度。本方法结合离散矩阵的稀疏性与管网结构特点,利用循环三对角矩阵求解算法与开源向量库,构建了管网中管道节点分开模拟的方法,从而加快了计算速度。
Description
技术领域
本发明涉及天然气管网数值模拟领域,具体为一种高效模拟大型天然气管网流动传热的方法。
背景技术
油气领域关系着我国的军民生活,随着经济的快速发展,我国对石油天然气类能源的需求越来越大,因此安全可靠的天然气管网建设更是关系着我国能源转型的关键因素之一。而建设安全可靠的天然气管网建设离不开准确高效的模拟气体在管网内的流动状态,如压力、流量、温度等。目前的管网模拟算法存在计算量大(牛顿迭代法或线性化方法等,需要计算雅可比矩阵,导致计算量过大)、离散系数矩阵求解速度慢(离散系数矩阵为大型不规则稀疏矩阵,通过迭代法求解)等因素,导致整体模拟速度慢,无法应用于大型天然气管网的模拟需求。
发明内容
本发明的目的是提出一种高效模拟大型天然气管网流动传热的的方法,本方法结合离散矩阵的稀疏性与管网结构特点,利用循环三对角矩阵求解算法与开源向量库,构建了管网中管道节点分开模拟的方法,从而加快了计算速度。
为达到上述目的,本发明采用以下技术方案:
一种高效模拟大型天然气管网流动传热的方法,包括以下步骤:
获取天然气管网的管网参数,该管网参数包括时间、管道轴向长度、管道直径、管道倾角、流入或流出管道的流体质量源、以及流体与环境之间的总换热系数;
建立天然气管网的一维控制方程,包括连续方程、动量方程和能量方程,并利用高斯定理将控制方程转换为积分形式;
根据管道的轴向长度,通过构建主网格和交错网格对天然气管网的所有管道进行离散;
对离散后的管道采用一阶迎风格式,将积分形式的控制方程在每个网格上进行离散,包括动量方程在交错网格上离散,连续方程和能量方程在主网格上离散;
将所有网格上的离散的控制方程中的系数根据网格顺序排列,得到一个类三对角的系数矩阵;
将天然气管网中的所有管道所对应的系数矩阵进行分块,每条管道所对应的系数矩阵分为一个块;
按照块提取每条管道所对应的系数矩阵,得到类循环三对角矩阵;并基于该类循环三对角矩阵进行快速求解,得到管道的两个端点和内部点的函数变化解;
利用函数变化解删除分块的系数矩阵中所有普通三对角矩阵元素,得到不规则稀疏矩阵;并采用稳定双共轭梯度法来求解该不规则稀疏矩阵,得到每条管道的端点值;
将每条管道的端点值代入到所述函数变化解中,得到每条管道的内部点值,由每条管道的端点值和内部点值构成天然气管网的解,实现对天然气管网流动传热模拟。
进一步地,构建主网格和交错网格的方法为:从每条管道的起点开始,以一预设长度为单位对整条管道进行划分,则划分出来的网格即为主网格;以两个相邻的主网格的中点为两边界,将该两边界内的区域作为一个交错网格。
进一步地,当在管道交汇处构建主网格和交错网格时,主网格设置在管道交汇处,交错网格设置在管道内部并与主网格错开半个网格。
进一步地,基于向量化函数库,迭代计算流体与管道之间的摩阻系数。
进一步地,利用SIMPLE算法,通过时间推进方法求解管网中流体速度、流体压力、流体速度、流体内能、流体焓以及流体温度参数随时间的变化数据。
本发明方法通过采用分离算法求解天然气管网控制方程,避免了耦合算法需要计算雅可比矩阵导致的计算量过大问题。再结合管网离散方程的特点,利用循环三对角矩阵快速求解算法将管网中的节点与管道解耦,在可以充分利用追赶法快速求解的特性的基础上,保证了管网中管道之间的正确联系。最后再利用开源向量化工具,加快非四则运算的速度,从而进一步加快求解速度。通过以上方法,形成了高效快速的大型天然气管网模拟方法。实验结果表明,本方法可以大幅提高模拟速度,其模拟速度是商业软件SPS的5-40倍,能够满足现场的需求。可应用于大型天然气管网的瞬态计算,服务于工业中的工艺设计、生产效率改进及节能降耗等方面。
附图说明
图1是使用交错网格对管网中管道进行离散图。
图2是管网节点采用离散网格图。
图3是类三对角的稀疏矩阵图。
图4是将稀疏矩阵进行分块图。
图5是类循环三对角矩阵图。
图6是管网拓扑结构图。
图7是管网入口流量变化趋势对比图。
图8是不同网格下的商业软件SPS耗时与本发明模拟的耗时比图。
具体实施方式
为使本发明的技术方案能更明显易懂,特举实施例并结合附图详细说明如下。
由于管道在轴向长度远大于管道直径,因此建立一维控制方程,包括连续方程、动量方程与能量方程:
连续方程为:动量方程为:能量方程为:其中t表示时间,x表示管道轴向长度,ρ表示流体密度,u表示流体速度,表示流入或流出管道的流体质量源,p表示流体压力,g表示重力加速度,θ表示管道倾角,λ表示流体与管道之间的摩阻系数,D表示管道直径,E表示流体内能,H表示流体焓,k表示流体与环境之间的总换热系数,Te表示环境温度,T表示流体温度,A表示管道横截面积。
流体与管道之间的摩阻系数λ采用柯列勃洛克公式:
其中Ke为管道绝对粗糙度,Re为雷诺数。
根据管道轴向长度,构建图1所示的主网格和交错网格对管网中所有管道进行离散,图1中,P表示当前节点,W表示西边节点,E表示东边节点,w表示西边界面,e表示东边界面。构建网格的方法为:从管道起点开始,采用一预设长度例如100~1000m为单位对整个管道进行划分,得到覆盖整条管道的虚拟网格,即为主网格(图1中相邻两条实线之间为一个主网格);将主网格中的中点(例如P点)为边界,则每两个相邻中点表示一个交错网格。
管网节点离散采用如图2所示离散网格。将主网格设置在管道交汇处,交错网格设置在管道内部(即不设置在管道交汇处),来避免速度在节点处出现多于两个方向的问题,即只有向前或向后的两个流动方向。图2中,实线框表示主网格,虚线框表示交错网格,数字表示网格界面,字母表示主网格中心。
在管道离散的基础上,采用一阶迎风格式(网格界面处需要但不在界面处而在节点处的量,采用来流方向上游节点处的值来代替)将积分形式的控制方程在每个网格上进行离散。其中动量方程在交错网格(见图1实线网格)上离散;将连续方程与能量方程在主网格(见图1虚线网格)上离散。
将所有网格上的离散的控制方程中的系数根据网格顺序排列,可得到一个类三对角的系数矩阵,如图3所示。
根据每条管道将图3所示的系数矩阵进行分块,即将管网中的每条管道对应的方程系数分为一块,如图4所示。
提取某条管道对应的系数矩阵(如图4中的某个虚线框),可得如图5所示的类循环三对角矩阵(标准三对角矩阵中基础上,第一行与最后一行有非零元素在非三对角上),其特征是除两个端点外,剩余矩阵为普通三对角矩阵。因此基于类循环三对角矩阵的快速求解算法,将两个端点(起点和终点)作为未知量,求解出管道两个端点(起点和终点)与内部点(除起点和终点以外的点)的函数变化解。求解过程如下:设某条管道的某个解用[x1,x2,…,xn-1,xn]来表示,且x1、xn为两个端点的解,则通过基于类循环三对角矩阵求解(求解内部三对角方程)可得到内部点的解为:x2=s2+a2x1+b2xn,…,xn-1=sn-1+an-1x1+bn-1xn,其中s表示x1=0、xn=0的解、a表示x1=1、xn=0的解、b表示x1=0、xn=1的解。
随后,删除图4所示的所有普通三对角矩阵元素,得到矩阵阶数大幅降低从而使得矩阵条件数也大幅降低的不规则稀疏矩阵,并采用稳定双共轭梯度法(BiCGStab)来求解此不规则稀疏矩阵,最终得到每条管道的端点值。然后代入上一步得到的管道两个端点与内部点的函数变化解x2=s2+a2x1+b2xn,…,xn-1=sn-1+an-1x1+bn-1xn中,得到内部点的值,最终得到整个管网的解。
由于在计算动量方程中的摩阻系数λ(摩阻项)时用到了隐式方程柯列勃洛克公式需要用到迭代法来求解,且该方程中又含有对数、开方等耗时操作,再结合在计算过程中,存在大量相同的操作,因此在计算摩阻时,利用基于硬件加速的开源向量化函数库(如Python中的Numpy等)来进一步加速求解速度,例如引入相应的向量化函数库,并将需要对数、开方的所有数据组成一个新数组并传入向量化函数库对应的函数来计算所需数据的对数、开方。
利用SIMPLE算法,通过时间推进方法求解管网中流体速度、流体压力、流体速度、流体内能、流体焓、流体温度等参数随时间的变化数据。
下面通过某管网算例来说明本发明的提速效果:图6为管网拓扑结构图,图7为模拟结果对比图,可以看出模拟结果基本相同,误差在可接受范围内;图8为加速比,可以看出加速效果非常显著,在5~40倍左右。结合现场实际应用时网格步长可取500m(500m对应图8中的第一个点)可知,本发明比商业软件快40倍左右,完全可以满足工业需求。
以上实施例仅用以说明本发明的技术方案而非对其进行限制,本领域的普通技术人员可以对本发明的技术方案进行修改或者等同替换,本发明的保护范围以权利要求所述为准。
Claims (10)
1.一种高效模拟大型天然气管网流动传热的方法,其特征在于,包括以下步骤:
获取天然气管网的管网参数,该管网参数包括时间、管道轴向长度、管道直径、管道倾角、流入或流出管道的流体质量源、以及流体与环境之间的总换热系数;
建立天然气管网的一维控制方程,包括连续方程、动量方程和能量方程,并利用高斯定理将控制方程转换为积分形式;
根据管道的轴向长度,通过构建主网格和交错网格对天然气管网的所有管道进行离散;
对离散后的管道采用一阶迎风格式,将积分形式的控制方程在每个网格上进行离散,包括动量方程在交错网格上离散,连续方程和能量方程在主网格上离散;
将所有网格上的离散的控制方程中的系数根据网格顺序排列,得到一个类三对角的系数矩阵;
将天然气管网中的所有管道所对应的系数矩阵进行分块,每条管道所对应的系数矩阵分为一个块;
按照块提取每条管道所对应的系数矩阵,得到类循环三对角矩阵;并基于该类循环三对角矩阵进行快速求解,得到管道的两个端点和内部点的函数变化解;
利用函数变化解删除分块的系数矩阵中所有普通三对角矩阵元素,得到不规则稀疏矩阵;并采用稳定双共轭梯度法来求解该不规则稀疏矩阵,得到每条管道的端点值;
将每条管道的端点值代入到所述函数变化解中,得到每条管道的内部点值,由每条管道的端点值和内部点值构成天然气管网的解,实现对天然气管网流动传热模拟。
4.如权利要求3所述的方法,其特征在于,基于向量化函数库,迭代计算流体与管道之间的摩阻系数λ。
6.如权利要求1所述的方法,其特征在于,利用SIMPLE算法,通过时间推进方法求解管网中流体速度、流体压力、流体速度、流体内能、流体焓以及流体温度参数随时间的变化数据。
7.如权利要求1所述的方法,其特征在于,构建主网格和交错网格的方法为:从每条管道的起点开始,以一预设长度为单位对整条管道进行划分,则划分出来的网格即为主网格;以两个相邻的主网格的中点为两边界,将该两边界内的区域作为一个交错网格。
8.如权利要求7所述的方法,其特征在于,预设长度为100~1000m。
9.如权利要求1或8所述的方法,其特征在于,当在管道交汇处构建主网格和交错网格时,主网格设置在管道交汇处,交错网格设置在管道内部并与主网格错开半个网格。
10.如权利要求1所述的方法,其特征在于,求取管道的两个端点和内部点的函数变化解的方法为:对于某条管道的函数变化解[x1,x2,…,xn-1,xn],其中x1、xn为两个端点的函数变化解,通过对类循环三对角矩阵求解,得到内部点的函数变化解即:x2=s2+a2x1+b2xn,…,xn-1=sn-1+an-1x1+bn-1xn,其中s表示x1=0,xn=0的解,a表示x1=1,xn=0的解,b表示x1=0,xn=1的解。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202111304959.0A CN113962131B (zh) | 2021-11-05 | 2021-11-05 | 一种高效模拟大型天然气管网流动传热的方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202111304959.0A CN113962131B (zh) | 2021-11-05 | 2021-11-05 | 一种高效模拟大型天然气管网流动传热的方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN113962131A true CN113962131A (zh) | 2022-01-21 |
CN113962131B CN113962131B (zh) | 2024-04-30 |
Family
ID=79469375
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202111304959.0A Active CN113962131B (zh) | 2021-11-05 | 2021-11-05 | 一种高效模拟大型天然气管网流动传热的方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN113962131B (zh) |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
NO20121424A1 (no) * | 2012-11-27 | 2014-05-28 | Sinvent As | Fremgangsmåte for simulering av flerfasefase fluidstrømninger i rørledninger |
CN104133958A (zh) * | 2014-07-28 | 2014-11-05 | 浙江中控软件技术有限公司 | 一种复杂管网模拟仿真计算方法及装置 |
US20150261893A1 (en) * | 2014-04-22 | 2015-09-17 | China University Of Petroleum - Beijing | Method and apparatus for determining pipeline flow status parameter of natural gas pipeline network |
CN106777753A (zh) * | 2016-12-29 | 2017-05-31 | 中国科学院工程热物理研究所 | 一种管网内外传热耦合仿真方法 |
CN110826188A (zh) * | 2019-10-14 | 2020-02-21 | 北京石油化工学院 | 一种基于gpu加速的天然气管网水力参数仿真方法 |
-
2021
- 2021-11-05 CN CN202111304959.0A patent/CN113962131B/zh active Active
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
NO20121424A1 (no) * | 2012-11-27 | 2014-05-28 | Sinvent As | Fremgangsmåte for simulering av flerfasefase fluidstrømninger i rørledninger |
US20150261893A1 (en) * | 2014-04-22 | 2015-09-17 | China University Of Petroleum - Beijing | Method and apparatus for determining pipeline flow status parameter of natural gas pipeline network |
CN104133958A (zh) * | 2014-07-28 | 2014-11-05 | 浙江中控软件技术有限公司 | 一种复杂管网模拟仿真计算方法及装置 |
CN106777753A (zh) * | 2016-12-29 | 2017-05-31 | 中国科学院工程热物理研究所 | 一种管网内外传热耦合仿真方法 |
CN110826188A (zh) * | 2019-10-14 | 2020-02-21 | 北京石油化工学院 | 一种基于gpu加速的天然气管网水力参数仿真方法 |
Non-Patent Citations (2)
Title |
---|
郑建国;陈国群;宋飞;艾慕阳;赵佳丽;: "大型天然气管网仿真模型及其求解技术研究", ***仿真学报, no. 06, 8 June 2012 (2012-06-08) * |
郑鹏宇;马群;叶永进;: "天然气管道不稳定流动数值模拟计算", 广东化工, no. 18, 30 September 2015 (2015-09-30) * |
Also Published As
Publication number | Publication date |
---|---|
CN113962131B (zh) | 2024-04-30 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN104461677B (zh) | 一种基于cfd和fem技术的虚拟热试验方法 | |
CN106503396B (zh) | 基于有限差分法与有限体积法耦合的多维水力***瞬变模拟方法 | |
CN111259547B (zh) | 一种用于综合能源***运行控制的天然气气路建模方法 | |
CN105302997A (zh) | 一种基于三维cfd的液柱分离-弥合水锤的模拟方法 | |
CN107947245B (zh) | 考虑天然气***约束的等值最优潮流模型构建方法 | |
CN110826188B (zh) | 一种基于gpu加速的天然气管网水力参数仿真方法 | |
Wansaseub et al. | Optimal U-shaped baffle square-duct heat exchanger through surrogate-assisted self-adaptive differential evolution with neighbourhood search and weighted exploitation-exploration | |
Luo et al. | Adaptive edge-based finite element schemes for the Euler and Navier-Stokes equations on unstructured grids | |
CN111859645B (zh) | 冲击波求解的改进musl格式物质点法 | |
WO2024016621A1 (zh) | 反应堆试验模型的规模确定方法、装置和计算机设备 | |
CN110276131B (zh) | 基于多项式响应面模型的翼身融合水下滑翔机外形优化方法 | |
CN115079592A (zh) | 一种船舶核动力装置热力***管网仿真方法 | |
CN115017808A (zh) | 一种基于改进蝴蝶算法优化hkelm的管道冲蚀预测方法 | |
CN110826261A (zh) | 一种基于fluent的埋地燃气管道泄漏模拟方法 | |
CN113962131A (zh) | 一种高效模拟大型天然气管网流动传热的方法 | |
CN117371354A (zh) | 一种天然气管网瞬稳态仿真方法、装置、设备及介质 | |
CN111950135A (zh) | 一种基于网络解耦的电-气互联***概率能流计算方法 | |
CN111783309A (zh) | 基于内部守恒的蒸汽供热网络动态仿真方法 | |
CN115795715A (zh) | 一种用于高温气冷堆换热装置热工水力的仿真方法及*** | |
CN115935566A (zh) | 天然气管网的模拟仿真方法、***、存储介质和电子设备 | |
CN113468785A (zh) | 一种高温超导电缆波纹管内液氮流动的模型仿真方法 | |
CN114117819A (zh) | 一种热蒸汽网络稳态仿真方法、***、设备及介质 | |
CN114936522B (zh) | 一种天然气***动态建模方法 | |
Liu et al. | Hull Resistance Performance Optimization Based on Fitting Surface Deformation | |
CN117217121B (zh) | 基于分布式并行topso-ewoa-de算法的岩体力学参数反演方法及*** |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |