CN102779198A - 天然气管网输送系统动态仿真系统建立方法 - Google Patents

Abstract

本发明是一种天然气管网输送系统动态仿真系统建立方法。涉及管道系统技术领域。它是将管段及非管元件数学模型联立，离散化后进行线性化处理，形成工程上等价、可直接求解的动态仿真系统线性数值模型；其流程为：建立天然气管网输送系统动态仿真数学模型；天然气管网输送系统动态仿真数学模型离散化；离散化各节点方程后重排序；线性化处理；得到天然气管网输送系统动态仿真数值模型。本发明建立任意拓扑结构准确度和计算效率高。

Description

天然气管网输送系统动态仿真系统建立方法

技术领域

本发明是一种天然气管网输送系统动态仿真系统建立方法。涉及管道系统技术领域。

背景技术

天然气管网输送系统仿真包括稳态仿真和动态仿真两方面，稳态仿真是进行管网设计的有力工具，动态仿真能够模拟管网运行参数随时间的变化，可以针对不同的工艺设计方案进行瞬态计算，保障有效的调度管理和及时发现并处理突发事故。天然气管网输送系统动态仿真与稳态仿真相比求解难度更大，应用范围更广。而能否对管网系统进行较好的动态仿真，关键在于建立能够准确完整地描述管网动态运行时各种状态的数学模型。

2005年“天津大学论文”中“燃气管网系统安全性及仿真的理论分析与应用研究”公开了一种以三大守恒方程为基础建立等温和非等温条件下的燃气管网稳动态仿真理论模型，通过对比所得结果进行比较和分析，得出其工程的可应用性，进而为燃气管网的优化提供参考。

完整的天然气管网输送系统包括管段和非管元件等。经典的天然气管网输送系统动态仿真数值模型已经可以较为准确描述灌完动态运行时的各种状况，但整体求解也比较困难，计算效率较低。对于建立满足生产应用的动态仿真数值模型，必须权衡模型准确度和计算效率两方面问题。

发明内容

本发明的目的是发明一种建立任意拓扑结构准确度和计算效率高的天然气管网输送系统动态仿真系统建立方法。

本发明提出的适用于任意拓扑结构天然气管网输送系统动态仿真系统建立方法，是将管段及非管元件数学模型联立，离散化后进行线性化处理，形成工程上等价、可直接求解的动态仿真系统线性数值模型。

本发明的流程如下(见图1)：

建立天然气管网输送系统动态仿真数学模型；

天然气管网输送系统动态仿真数学模型离散化；

离散化各节点方程后重排序；

线性化处理；

得到天然气管网输送系统动态仿真数值模型。

所述建立天然气管网输送系统动态仿真数学模型是：

天然气管网输送系统主要包括管段和阀门、压缩机以及分支点等非管元件；

1.管段的控制方程

1)质量守恒控制方程

(Aρ)_t+(Aρv)_x＝0，0≤x≤L；t≥0    (1)

式中：L为管长，m；x为距离，m；t为时间，s；A管道截面积，m²；ρ为气体密度，kg/m³；v为气体流速，m/s。

2)动量守恒控制方程

$v_t+{\mathrm{vv}}_x+\frac{\mathrm{Px}}{\mathrm{\ρ}}+{\mathrm{gh}}_x+\frac{f}{{2D}_i}v\left|v\right|=\mathrm{0,0}\≤x\≤L;t\≥0---\left(2\right)$

式中：f为摩阻系数，无单位；Di为管道内径，m；g为重力加速度，m/s²；h为管段高程，m。

3)能量守恒控制方程

${\mathrm{\ρc}}_v(T_t+{\mathrm{vT}}_x)=-T{\left(\frac{\∂P}{\∂T}\right)}_{\mathrm{\ρ}}{v}_x+\mathrm{\ρ}\frac{f}{{2D}_i}{\left|v\right|}^3-\frac{{4U}_w}{D_i}(T-T_g),0\≤x\≤L,t\≥0---\left(3\right)$

式中：p为气体压力，Pa；T为气体温度，℃；U_w为总传热系数，W/(m².℃.)；Tg为土壤温度，℃；c_v为气体比热容，J/(kg.℃)；其它变量意义同质量守恒控制方程和动量守恒控制方程。

2.阀门的控制方程

1)质量守恒控制方程

A_upρ_upv_up＝A_dwρ_dwv_dw              (4)

式中下标up表示上游变量，dw表示下游变量。

2)压降特性方程

$A_{\mathrm{up}}{\mathrm{\ρ}}_{\mathrm{up}}{V}_{\mathrm{up}}=C_g\sqrt{(P_{\mathrm{up}}-P_{\mathrm{dw}})P_{\mathrm{dw}}}---\left(5\right)$

式中：C_g为阀门阻抗系数。

3)能量守恒方程

$T_{\mathrm{dw}}=T_{\mathrm{up}}+(P_{\mathrm{up}}-P_{\mathrm{dw}})\frac{1}{C_P}[\frac{T}{{\mathrm{\ρ}}^2}\frac{{\left(\frac{\∂P}{\∂T}\right)}_{\mathrm{\ρ}}}{{\left(\frac{\∂P}{\∂\mathrm{\ρ}}\right)}_T}-\frac{1}{\mathrm{\ρ}}]---\left(6\right)$

式中c_p为定压比热容，kJ/(kmol·K)。

3.离心式压缩机的控制方程

1)质量守恒控制

A_dwρ_dwv_dw-A_upρ_upv_up＝M_fuel    (7)

式中：M_fuel为燃料消耗。

2)增压特性方程

p_dw＝εp_up                      (8)

式中：ε为压比。

3)升温特性方程

$T_{\mathrm{dw}}=T_{\mathrm{up}}[1+({\mathrm{\ϵ}}^{\frac{m-1}{m}}-1)/\mathrm{\η}]---\left(9\right)$

式中：m为多方指数，η为绝热效率。

${\mathrm{\ϵ}}^2=a_1{\left(\frac{n}{n_0}\right)}^2+b_1\left(\frac{n}{n_0}\right)Q+c_1{Q}^2---\left(10\right)$

4.分支点的控制方程

设分支管为m进n出型。

1)连续性方程

$\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\mathrm{Inpour}{M}_i=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{Outpour}{M}_j$

2)动量方程

$\∀1\≤i\≤m,$ $\∀1\≤j\≤n,$

Inpourp_i＝Outpourp_j

3)能量方程

m＝1时， $\∀1\≤j\≤n,$

InpourT＝OutpourT_j

m＞1时， $\∀2\≤j\≤n,$

OutpourT₁＝OutpourT_j

$\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\mathrm{Inpour}{C}_{p_i}{M}_i{T}_i=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\mathrm{Outpour}{C}_{p_j}{M}_j{T}_j$

Cp定压比热容，P为压力，M质量流量，T温度；

所述天然气管网输送系统动态仿真数学模型离散化主要是根据管网拓扑顺序对管段进行离散，将非管元件视为边界离散；

天然气管网输送系统动态仿真计算较为复杂，要想求得方程组的解析解相当困难，只能退而求其次求得方程组的数值解，这需要将数学模型变成数值模型；离散化就是将连续问题的解用一组离散要素来表征而近似求解的方法，天然气管网输送系统动态仿真数学模型离散化是整个仿真系统数值模型建立的关键一步，这一步将工艺流程建模的计算机表示(有向图)和描述气体在管道中流动的三大控制方程、气体状态方程、各类设备的特性方程以及外部边界的状态方程联系起来转换为一个数值系统进行求解；

1.管段的离散化

根据天然气管网拓扑排序的顺序按照一定的空间步长依次对每个管段进行空间离散，得到一系列计算节点，如图2所示。针对每个计算节点，将其流动的质量守恒、动量守恒、能量守恒三大控制方程分别用数值形式表示即完成管段的离散化；

2.边界的离散化

天然气管网输送系统中的边界包括外部边界和内部边界两种；外部边界主要指气源和分输点等，而内部边界主要指阀门、压缩机和分支点等；对于所有边界均视为一个计算节点，外部边界和内部边界状态方程离散化不同的地方在于，外部边界节点状态方程只是上游或者下游中的一种情况，而内部边界节点状态方程则包括上下游节点间的约束关系；

上述离散化处理过程属于经典的方法，具体过程不在详述；将以上每个计算节点所得的控制方程及边界约束关系式联立起来形成一个封闭的方程组；原则上到此一般建模方法中已完成天然气管网的数值离散，管段和边界等节点状态方程排列顺序不影响离散方程组解的存在性；但不同的排列顺序对于求解时间和效率的影响很大，为了便于维护和计算并提高计算效率，本发明提出建立数学模型方法中增加了按照以下顺序对联立后方程组进行排序的工作：

(1)管段

按照有向图的深度优先搜索算法对元件进行拓扑排序，按照在拓扑序列中出现的顺序，依次对管段的空间离散计算节点进行排序；

(2)内部边界建立的计算节点间的约束关系式；

(3)除分支点外外部边界建立的计算节点间的约束关系式；

(4)分支点建立的计算节点间的约束关系式；

所述天然气管网输送系统动态仿真数值模型线性化是：

对于离散化后形成的动态仿真熟悉模型的线性化处理最早见于国外对管段控制方程的处理，国内未见类似处理，且本方法之前国内外均未见对非管元件方程线性化处理；线性化的处理方法是将数值模型中的非线性部分线性化，主要优点是提高了计算速度快，虽然线性化后损失了部分计算精度，但从提出本建模方法目的考虑——建立一个满足生产实际应用以及长时间范围内模拟管网运行等需求的天然气管网输送系统动态仿真数值模型，线性化处理还是适合的；

1.管段流动方程线性化处理

首先对离散化后偏微分方程组按自变量重写，都写成以p、v、T为自变量的方程，形式如下：

u_t+Au_x＝F    (11)

下面将非线性项Au_x进行线性化处理。根据下式

A(v)w＝A(v₀)w+[D_vA(v₀)；w₀](v-v₀)+O(||v-v₀||²+||w-w₀||²)

其中使用[D_uA(x，v)；z]来表示位于(i，j)的元素为

${[D_{u}A(x,v);z]}_{i,j}=\underset{l=1}{\overset{n}{\mathrm{\Σ}}}{D}_{j+1}{A}_{i,l}(x,v)z_l---\left(12\right)$

的n×n维矩阵，

$D_u{H}_k=(\frac{{\∂H}_k}{{\∂u}_1},\frac{{\∂H}_k}{{\∂u}_2},L,\frac{{\∂H}_k}{{\∂u}_n})=(D_2{H}_k,D_3{H}_k,L,D_{n+1}{H}_k)---\left(13\right)$

舍掉高阶项便得到管道流动方程线性化处理后的结果；

2.边界特性方程线性化处理

边界特性方程即主要非管元件状态方程；在常见天然气管网输送系统动态仿真模型数值求解方法中，如特征线法和隐式差分法等，直接将偏微分方程转化成差分方程，然后用数值迭代法求解，优点是计算精度高而且计算较灵活，缺点是求解速度慢；而本发明是针对管道生产实际使用的仿真软件，对计算时间和效率有较高要求，这也是考虑对各方程进行线性化处理的初衷；边界特性方程线性化处理的思想和过程与管段流动方程线性化一致，也是先将各方程按基本变量重写，然后利用泰勒展开舍掉高阶无穷小量，保留线性项；

经过以上“线性化”处理后，原有的非线性方程组转化成为一个在工程上等价的线性方程组，可以直接进行求解，且线性方程组的系数具有明显的物理意义，实现了整个水力、热力系统的“线性化”。

本发明提出的适用于任意拓扑结构天然气管网输送系统动态仿真数值模型建立方法，将管段及非管元件数学模型联立，离散化后对描述天然气在管道内的流动方程和通过非管元件(气源、分输站、阀门、压缩机、换热器)时的动力特性方程做线性化预处理，并按一定原则对节点方程进行重排序，将描述大型天然气管网输送工艺的非线性方程组转换为线性方程组，将复杂的非线性问题转化为线性问题求解。

本发明可对任意拓扑结构的天然气管网输送系统进行建模，模型精度和计算复杂度满足生产实际应用水平。

本建模方法适用于国内现存的任意拓扑结构的天然气管网输送系统动态仿真建模，并可推广到液体管道输送系统的动态仿真领域，满足生产应用前景广阔。

附图说明

图1天然气管网输送系统动态仿真系统建立流程图

图2管段离散化示意图

图3天然气管网拓扑结构图

具体实施方式

实施例.以本例来说明本发明的具体实施方式并对本发明作进一步的说明。本例是一实验方法，其流程如图1所示。

1.本实施例针对的天然气管网拓扑结构如下图3所示。外部边界情况如表2所示，分支点情况如表1所示。

表1  实施例外部边界情况表

序号	外部边界	关联元件
			1	ES001	BV001，BV002，BV003
2	ED002	BV018
			3	ED003	BV019

表2  实施例分支点情况表

序号	分支点	上游关联元件	下游关联元件
				1	ND010	BV004，BV005	BV007，BV008
2	ND001	ES001	BV001，BV002，BV003
				3	ND014	BV009	ICC002，IRV002
4	ND017	BV010	ICC003，IRV003
				5	ND015	ICC002，IRV002	CV001
6	ND018	ICC003，IRV003	CV002
				7	ND007	BV001，BV013	PL001
8	ND004	BV011，BV012	BV013
				9	ND006	PL001	BV014，BV015
10	ND019	BV016，BV017	BV018
				11	ND005	BV014	BV017，PL002

2.动态仿真数学模型的离散化

将控制方程和状态方程按照分支、外部边界、分支点的顺序依次离散，此步骤属于经典方法，下面只给出关键的离散符号说明。

(1)设N为正整数，记Δx＝1/N，

x_j＝jΔx，j＝0，1，L，N

$x_{j+1/2}=(j+\frac{1}{2})\mathrm{\Δx},j=\mathrm{0,1},L,N-1$

其中，x_j：空间节点，x_j+1/2：中间节点。

对任意函数u(x)，记u_j＝u(x_j)，u_j+1/2＝u(x_j+1/2)，

(2)设M为正整数，记Δt＝1/M，

tⁿ＝nΔt，n＝0，1，L，M

t^n+θ＝(n+θ)Δt，n＝0，1，L，M-1

对任意函数u(t)，θ∈[0，1]，记uⁿ＝u(tⁿ)，u^n+θ＝u(t^n+θ)，u^n，θ＝θuⁿ⁺¹+(1-θ)uⁿ

(3)对任意函数u(x)，定义对x的偏导函数：

${\∂}_x{u}_{j+1/2}=\frac{u_{j+1}-u_j}{\mathrm{\Δx}},$

(4)对任意函数u(t)，θ∈[0，1]，定义对t的偏导函数：

${\∂}_t{u}^{n+\mathrm{\θ}}=\frac{u^{n+1}-u^n}{\mathrm{\Δt}}$

3.对离散化后各节点方程按照如下拓扑顺序重新排列：

ES001，ND001，BV003，ND009，BV005，BV002，ND008，BV004，ND010，BV008，ND013，BV010，ND017，IRV003，ICC003，ND018，CV002，ND003，BV012，BV007，ND012，BV009，ND014，IRV002，ICC002，ND015，CV001，ND002，BV011，BV001，ND007，BV013，ND004，PL001，ND006，BV015，ND016，BV016，BV014，ND005，PL002，ND021，BV019，ND022，ED003，BV017，ND019，BV018，ND020，ED002

4.线性化处理：

实施例中式(11)中各项具体为：

$u=\left[\begin{array}{c}\mathrm{Ap}\\ v\end{array}\right],$ $A=\left[\begin{array}{cc}v& A\frac{\mathrm{\ρ}}{{\mathrm{\ρ}}_p}\\ \frac{1}{\mathrm{A\ρ}}& v\end{array}\right],$ $F=\left[\begin{array}{c}0\\ -\frac{f}{{2D}_i}v\left|v\right|\end{array}\right].$

根据式(12)和(13)，离散之后的动态仿真数学模型经过线性化处理变为如下形式的线性方程组：

$(\frac{1}{2\mathrm{Vt}}{I}_n+\frac{A\left(U_{i+1/2}^j\right)}{2\mathrm{Vx}}+\frac{1}{4}G(.)-\frac{1}{4}\frac{\∂F\left(U_{i+1/2}^j\right)}{{\∂u}_i})U_{i+1}^{j+1}$

$+(\frac{1}{2\mathrm{Vt}}{I}_n-\frac{A\left(U_{i+1/2}^j\right)}{2\mathrm{Vx}}+\frac{1}{4}G(.)-\frac{1}{4}\frac{\∂F\left(U_{i+1/2}^j\right)}{{\∂u}_i})U_i^{j+1}$

$=F\left(U_{i+1/2}^j\right)+(\frac{1}{2\mathrm{Vt}}{I}_n-\frac{A\left(U_{i+1/2}^j\right)}{2\mathrm{Vx}}-\frac{1}{4}G(.)+\frac{1}{4}\frac{\∂F\left(U_{i+1/2}^j\right)}{{\∂u}_i})U_{i+1}^j+\′$

$(\frac{1}{2\mathrm{Vt}}{I}_n+\frac{A\left(U_{i+1/2}^j\right)}{2\mathrm{Vx}}-\frac{1}{4}G(.)+\frac{1}{4}\frac{\∂F\left(U_{i+1/2}^j\right)}{{\∂u}_i})U_i^j+(G(.)-\frac{\∂F\left(U_{i+1/2}^j\right)}{{\∂u}_i})U_{i+1/2}^j$

其中：0≤i≤N-1，0≤j≤M-1，

$G(.)=\left[\begin{array}{cc}(1-\frac{\mathrm{\ρ}{\left({\mathrm{\ρ}}_P\right)}_p}{{\left({\mathrm{\ρ}}_p\right)}^2})\·\frac{v_{i+1}^j-v_i^j}{\mathrm{\Δx}}& A\frac{p_{i+1}^j-p_i^j}{\mathrm{\Δx}}\\ -\frac{{\mathrm{\ρ}}_P}{A^2{\mathrm{\ρ}}^2}\·A\frac{p_{i+1}^j-p_i^j}{\mathrm{\Δx}}& \frac{v_{i+1}^j-v_i^j}{\mathrm{\Δx}}\end{array}\right]$

$\frac{\∂F}{\∂u}=\left[\begin{array}{cc}0& 0\\ -\frac{f_p}{2A\left(D_i\right)}v\left|v\right|& -\frac{f}{\left(D_i\right)}\left|v\right|-\frac{f_v}{2\left(D_i\right)}v\left|v\right|\end{array}\right].$

用

表示质量守恒方程中

项的系数，

表示动量守恒方程中

项的系数，

表示边界约束方程中

的系数，整理后最终得到的线性方程组系数矩阵如下式所示，最后两行表示边界条件。

其中：

$C_{p_i^{j+1}}=\frac{1}{2\mathrm{\Δt}}-\frac{v_{i+1/2}^j}{2\mathrm{\Δx}}+\frac{1}{4}\·(1-\frac{{\mathrm{\ρ}\left({\mathrm{\ρ}}_p\right)}_p}{{\left({\mathrm{\ρ}}_p\right)}^2}){|}_{(p_{i+1/2}^j,v_{i+1/2}^j)}\frac{v_{i+1}^j{v}_i^j}{\mathrm{\Δx}}$

$C_{v_i^{j+1}}=-\frac{1}{2\mathrm{\Δx}}\·A\frac{\mathrm{\ρ}}{{\mathrm{\ρ}}_p}{|}_{(p_{i+1/2}^j,v_{i+1/2}^j)}+\frac{1}{4}\·A\frac{p_{i+1}^j-P_i^j}{\mathrm{\Δx}}$

$C_{p_{i+1}^{j+1}}=\frac{1}{2\mathrm{\Δt}}+\frac{v_{i+1/2}^j}{2\mathrm{\Δx}}+\frac{1}{4}\·(1-\frac{\mathrm{\ρ}{\left({\mathrm{\ρ}}_p\right)}_p}{{\left({\mathrm{\ρ}}_p\right)}^2}){|}_{(p_{i+1/2}^j,v_{i+1/2}^j)}\frac{v_{i+1}^j-v_i^j}{\mathrm{\Δx}}$

$C_{v_{i+1}^{j+1}}=\frac{1}{2\mathrm{\Δx}}\·A\frac{\mathrm{\ρ}}{{\mathrm{\ρ}}_p}{|}_{(p_{i+1/2}^j,v_{i+1/2}^j)}+\frac{1}{4}\·A\frac{p_{i+1}^j-p_i^j}{\mathrm{\Δx}}$

$M_{p_i^{j+1}}=-\frac{1}{2\mathrm{\Δx}}\·\frac{1}{A{\mathrm{\ρ}}_{i+1/2}^j}-\frac{1}{4}\·\left(\frac{{\mathrm{\ρ}}_p}{A^2{\mathrm{\ρ}}^2}\right){|}_{(p_{i+1/2}^j,v_{i+1/2}^j)}A\frac{p_{i+1}^j-p_i^j}{\mathrm{\Δx}}+\frac{1}{4}\·\left(\frac{\left|v\right|v}{2\mathrm{AD}}{f}_p\right){|}_{(p_{i+1/2}^j,v_{i+1/2}^j)}$

$M_{v_{\text{i}}^{j+1}}=\frac{1}{2\mathrm{\Δt}}-\frac{v_{i+1/2}^j}{2\mathrm{\Δx}}+\frac{1}{4}\·\frac{v_{i+1}^j{v}_i^j}{\mathrm{\Δx}}+\frac{1}{4}\·(\frac{\left|v\right|}{\left(D_i\right)}f+\frac{v\left|v\right|}{2\left(D_i\right)}{f}_v){|}_{(p_{i+1/2}^j,v_{i+1/2}^j)}$

$M_{p_{i+1}^{j+1}}=\frac{1}{2\mathrm{\Δx}}\·\frac{1}{\mathrm{A\ρ}}{|}_{(p_{i+1/2}^j,v_{i+1/2}^j)}+\frac{1}{4}\·(-\frac{{\mathrm{\ρ}}_p}{A^2{\left(\mathrm{\ρ}\right)}^2}){|}_{(p_{i+1/2}^j,v_{i+1/2}^j)}\·A\frac{p_{i+1}^j-p_i^j}{\mathrm{\Δx}}+\frac{1}{4}\·\left(\frac{\left|v\right|v}{2\mathrm{AD}}{f}_p\right){|}_{(p_{i+1/2}^j,v_{i+1/2}^j)}$

$M_{v_{i+1}^{j+1}}=\frac{1}{2\mathrm{\Δt}}+\frac{1}{2\mathrm{\Δx}}\·v_{i+1/2}^j+\frac{1}{4}\·\frac{v_{i+1}^j-v_i^j}{\mathrm{\Δx}}+\frac{1}{4}\·(\frac{\left|v\right|}{\left(D_i\right)}f+\frac{v\left|v\right|}{2\left(D_i\right)}{f}_v){|}_{(p_{i+1/2}^j,v_{i+1/2}^j)}$

$B_{p_0^{j+1}}=1,$ $B_{p_{N-1}^{j+1}}=1,$ 0≤i≤N-1，0≤j≤M-1。

记线性方程组右端常数项形式为：[C₀ C₁ C₂ L  C_2N-3  C_2N-2  C_2N-1]^T，其中偶数项和奇数项分别为(0≤k≤N-1)：

$C_{2k}=\frac{1}{2\mathrm{\Δt}}{p}_{i+1}^j-\frac{1}{2\mathrm{\Δx}}(v_{i+1/2}^j{p}_{i+1}^j+\frac{{2\mathrm{\ρ}}_p}{\mathrm{\ρ}}{v}_{i+}^j){|}_{(p_{i+1/2}^j,v_{i+1/2}^j)}$

$-\frac{1}{4}\left(2(1-\frac{\mathrm{\ρ}{\left({\mathrm{\ρ}}_p\right)}_p}{{\left({\mathrm{\ρ}}_p\right)}^2})\right){|}_{(p_{i+1/2}^j,v_{i+1/2}^j)}\frac{v_{i+1}^j-v_i^j}{\mathrm{\Δx}}\·p_{i+1}^j+\frac{p_{i+1}^j-p_i^j}{\mathrm{\Δx}}\·v_{i+1}^j$

$+\frac{1}{2\mathrm{\Δt}}{p}_i^j+\frac{1}{2\mathrm{\Δx}}(v_{i+1/2}^j{p}_i^j+\frac{{2\mathrm{\ρ}}_p}{\mathrm{\ρ}}{v}_i^j){|}_{(p_{i+1}^j,v_{i+1/2}^j)}$

$-\frac{1}{4}(2(1-\frac{\mathrm{\ρ}{\left({\mathrm{\ρ}}_p\right)}_p}{{\left({\mathrm{\ρ}}_p\right)}^2}){|}_{(p_{i+1/2}^j,v_{i+1/2}^j)}\frac{v_{i+1}^j-v_i^j}{\mathrm{\Δx}}\·p_i^j+\frac{p_{i+1}^j-p_i^j}{\mathrm{\Δx}}\·v_i^j)$

$+(2(1-\frac{\mathrm{\ρ}{\left({\mathrm{\ρ}}_p\right)}_p}{{\left({\mathrm{\ρ}}_p\right)}^2}){|}_{(p_{i+1/2}^j,v_{i+1/2}^j)}\frac{v_{i+1}^j-v_i^j}{\mathrm{\Δx}}\·p_{i+1/2}^j+\frac{p_{i+1}^j-p_i^j}{\mathrm{\Δx}}\·v_{i+1/2}^j)$

$C_{2k+1}=(-\frac{\left(t\right)}{2\left(D_i\right)}v|v\left|\right){|}_{(p_{i+1/2}^j,v_{i+1/2}^j)}+\frac{1}{2\mathrm{\Δt}}{v}_{i+1}^j-\frac{1}{2\mathrm{\Δx}}(\frac{1}{\mathrm{\ρ}}{p}_{i+1}^j+{\mathrm{vv}}_{i+1}^j){|}_{(p_{i+1/2}^j,v_{i+1/2}^j)}$

$-\frac{1}{4}((-\frac{{\mathrm{\ρ}}_p}{{\mathrm{\ρ}}^2}){|}_{(p_{i+1/2}^j,v_{i+1/2}^j)}\frac{p_{i+1}^j-p_i^j}{\mathrm{\Δx}}\·p_{i+1}^j+\frac{v_{i+1}^j+v_i^j}{\mathrm{\Δx}}\·v_{i+1}^j)$

$-\frac{1}{4}(\frac{v\left|v\right|}{2\left(D_i\right)}{f}_p\·p_{i+1}^j+(\frac{\left|v\right|}{\left(D_i\right)}f+\frac{v\left|v\right|}{2\left(D_i\right)}{f}_v)\·v_{i+1}^j){|}_{(p_{i=1/2}^j,v_{i+1/2}^j)}$

$+\frac{1}{2\mathrm{\Δt}}{v}_i^j+\frac{1}{2\mathrm{\Δx}}(\frac{1}{\mathrm{\ρ}}{p}_i^j+{\mathrm{vv}}_i^j){|}_{(p_{i+1/2}^j,v_{i+1/2}^j)}$

$-\frac{1}{4}((-\frac{{\mathrm{\ρ}}_p}{{\mathrm{\ρ}}^2}){|}_{(p_{i+1/2}^j,v_{i+1/2}^j)}\frac{p_{i+1}^j-p_i^j}{\mathrm{\Δx}}\·p_i^j+\frac{v_{i+1}^j-v_i^j}{\mathrm{\Δx}}\·v_i^j)$

$-\frac{1}{4}(\frac{v\left|v\right|}{2\left(D_i\right)}{f}_p\·p_i^j+(\frac{\left|v\right|}{\left(D_i\right)}f+\frac{v\left|v\right|}{2\left(D_i\right)}{f}_v)\·v_i^j){|}_{(p_{i+1/2}^j,v_{i+1/2}^j)}$

$+((-\frac{{\mathrm{\ρ}}_p}{{\mathrm{\ρ}}^2}){|}_{(p_{i+1/2}^j,v_{i+1/2}^j)}\frac{p_{i+1}^j-p_i^j}{\mathrm{\Δx}}\·p_{i+1/2}^j+\frac{v_{i+1}^j-v_i^j}{\mathrm{\Δx}}\·v_{i+1/2}^j)$

$+(\frac{v\left|v\right|}{2\left(D_i\right)}{f}_p\·p_{i+1/2}^j+(\frac{\left|v\right|}{\left(D_i\right)}f+\frac{v\left|v\right|}{2\left(D_i\right)}{f}_v)\·v_{i+1/2}^j){|}_{(p_{i+1/2}^j,v_{i+1/2}^j)}$

至此，实施例中对应的天然气管网输送系统动态仿真系统建立完毕。

本例对任意拓扑结构的天然气管网输送系统进行建模，模型精度和计算复杂度满足生产实际应用水平。

Claims (4)

1.一种天然气管网输送系统动态仿真系统建立方法，其特征是将管段及非管元件数学模型联立，离散化后进行线性化处理，形成工程上等价、可直接求解的动态仿真系统线性数值模型；

其流程为：

建立天然气管网输送系统动态仿真数学模型；

天然气管网输送系统动态仿真数学模型离散化；

离散化各节点方程后重排序；

线性化处理；

得到天然气管网输送系统动态仿真数值模型。

2.根据权利要求1所述的天然气管网输送系统动态仿真系统建立方法，其特征是所述建立天然气管网输送系统动态仿真数学模型是：

天然气管网输送系统主要包括管段和阀门、压缩机以及分支点非管元件；

1)管段的控制方程

(1)质量守恒控制方程

(Aρ)_t+(Aρv)_x＝0，0≤x≤L；t≥0    (1)

式中：L为管长，m；x为距离，m；t为时间，s；A管道截面积，m²；ρ为气体密度，kg/m³；v为气体流速，m/s。

(2)动量守恒控制方程

$v_t+{\mathrm{vv}}_x+\frac{p_x}{\mathrm{\ρ}}+{\mathrm{gh}}_x+\frac{f}{{2D}_i}v\left|v\right|=\mathrm{0,0}\≤x\≤L;t\≥0---\left(2\right)$

式中：f为摩阻系数，无单位；D_i为管道内径，m；g为重力加速度，m/s²；h为管段高程，m。

(3)能量守恒控制方程

${\mathrm{\ρc}}_v(T_t+{\mathrm{vT}}_x)=-T{\left(\frac{\∂p}{\∂T}\right)}_{\mathrm{\ρ}}{v}_x+\mathrm{\ρ}\frac{f}{{2D}_i}|v{|}^3-\frac{{4U}_w}{D_i}(T-T_g),0\≤x\≤L,t\≥0---\left(3\right)$

式中：p为气体压力，Pa；T为气体温度，℃；U_w为总传热系数，W/(m².℃.)；Tg为土壤温度，℃；c_v为气体比热容，J/(kg.℃)；其它变量意义同质量守恒控制方程和动量守恒控制方程。

2)阀门的控制方程

(1)质量守恒控制方程

A_upρ_upv_up＝A_dwρ_dwv_dw          (4)

式中下标up表示上游变量，dw表示下游变量。

(2)压降特性方程

$A_{\mathrm{up}}{\mathrm{\ρ}}_{\mathrm{up}}{v}_{\mathrm{up}}=C_g\sqrt{(p_{\mathrm{up}}-p_{\mathrm{dw}})p_{\mathrm{dw}}}---\left(5\right)$

式中：C_g为阀门阻抗系数。

(3)能量守恒方程

$T_{\mathrm{dw}}=T_{\mathrm{up}}+(P_{\mathrm{up}}-P_{\mathrm{dw}})\frac{1}{C_p}[\frac{T}{{\mathrm{\ρ}}^2}\frac{{\left(\frac{\∂p}{\∂T}\right)}_{\mathrm{\ρ}}}{{\left(\frac{\∂p}{\∂\mathrm{\ρ}}\right)}_T}-\frac{1}{\mathrm{\ρ}}]---\left(6\right)$

式中c_p为定压比热容，kJ/(kmol·K)。

3)离心式压缩机的控制方程

(1)质量守恒控制

A_dwρ_dwv_dw-A_upρ_upv_up＝M_fuel    (7)

式中：M_fuel为燃料消耗。

(2)增压特性方程

p_dw＝εp_up                      (8)

式中：ε为压比。

(3)升温特性方程

$T_{\mathrm{dw}}=T_{\mathrm{up}}[1+({\mathrm{\ϵ}}^{\frac{m-1}{m}}-1)/\mathrm{\η}]---\left(9\right)$

式中：m为多方指数，η为绝热效率。

${\mathrm{\ϵ}}^2=a_1{\left(\frac{n}{n_0}\right)}^2+b_1\left(\frac{n}{n_0}\right)Q+c_1{Q}^2---\left(10\right)$

4)分支点的控制方程

设分支管为m进n出型。

(1)连续性方程

$\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\mathrm{Inpour}{M}_i=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{OutpourM}}_j$

(2)动量方程

$\∀1\≤i\≤m,$ $\∀1\≤j\≤n,$

Inpourp_i＝Outpourp_j

(3)能量方程

m＝1时，

InpourT＝OutpourT_j

m＞1时，

OutpourT₁＝OutpourT_j

$\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{InpourC}}_{p_i}{M}_i{T}_i=\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{OutpourC}}_{p_j}{M}_j{T}_j$

C_p定压比热容，P为压力，M质量流量，T温度；

3.根据权利要求1所述的天然气管网输送系统动态仿真系统建立方法，其特征是所述天然气管网输送系统动态仿真数学模型离散化是根据管网拓扑顺序对管段进行离散，将非管元件视为边界离散；

1)管段的离散化

根据天然气管网拓扑排序的顺序按照一定的空间步长依次对每个管段进行空间离散，得到一系列计算节点；针对每个计算节点，将其流动的质量守恒、动量守恒、能量守恒三大控制方程分别用数值形式表示即完成管段的离散化；

2)边界的离散化

天然气管网输送系统中的边界包括外部边界和内部边界两种；外部边界主要指气源和分输点，而内部边界主要指阀门、压缩机和分支点；对于所有边界均视为一个计算节点，外部边界和内部边界状态方程离散化不同的地方在于，外部边界节点状态方程只是上游或者下游中的一种情况，而内部边界节点状态方程则包括上下游节点间的约束关系；

将以上每个计算节点所得的控制方程及边界约束关系式联立起来形成一个封闭的方程组；原则上到此一般建模方法中已完成天然气管网的数值离散，管段和边界等节点状态方程排列顺序不影响离散方程组解的存在性；按照以下顺序对联立后方程组进行排序的工作：

(1)管段

按照有向图的深度优先搜索算法对元件进行拓扑排序，按照在拓扑序列中出现的顺序，依次对管段的空间离散计算节点进行排序；

(2)内部边界建立的计算节点间的约束关系式；

(3)除分支点外外部边界建立的计算节点间的约束关系式；

(4)分支点建立的计算节点间的约束关系式；

4.根据权利要求1所述的天然气管网输送系统动态仿真系统建立方法，其特征是所述天然气管网输送系统动态仿真数值模型线性化是：

线性化的处理方法是将数值模型中的非线性部分线性化；

1)管段流动方程线性化处理

首先对离散化后偏微分方程组按自变量重写，都写成以p、v、T为自变量的方程，形式如下：

u_t+Au_x＝F                        (11)

下面将非线性项Au_x进行线性化处理。根据下式

A(v)_w＝A(v₀)w+[D_vA(v₀)；w₀](v-v₀)+O(||v-v₀||²+||w-w₀||²)

其中使用[D_uA(x，v)；z]来表示位于(i，j)的元素为

${[D_{u}A(x,v);z]}_{i,j}=\underset{l=1}{\overset{n}{\mathrm{\Σ}}}{D}_{j+1}{A}_{i,l}(x,v)z_l---\left(12\right)$

的n×n维矩阵，

$D_u{H}_k=(\frac{\∂H_k}{{\∂u}_1},\frac{{\∂H}_k}{{\∂u}_2},L,\frac{{\∂H}_k}{\∂u_n})=(D_2{H}_k,D_3{H}_k,L{,D}_{n+1}{H}_k)---\left(13\right)$

舍掉高阶项便得到管道流动方程线性化处理后的结果；

2)边界特性方程线性化处理

边界特性方程即主要非管元件状态方程；是先将各方程按基本变量重写，然后利用泰勒展开舍掉高阶无穷小量，保留线性项；

最终得到的离散化矩阵如下式所示，最后两行表示边界条件，

表示质量守恒方程中

项的系数，

表示动量守恒方程中

项的系数，

表示边界约束方程中

的系数；

$\left[\begin{array}{cccccccccccc}{C}_{p_0^{j+1}}& C_{v_0^{j+1}}& C_{p_1^{j+1}}& C_{v_1^{j+1}}& 0& L& & & & & & 0\\ M_{p_0^{j+1}}& M_{v_0^{j+1}}& M_{p_1^{j+1}}& M_{v_1^{j+1}}& 0& L& & & & & & 0\\ 0& 0& C_{p_1^{j+1}}& C_{v_1^{j+1}}& C_{p_2^{j+1}}& C_{v_2^{j+1}}& 0& L& & & & 0\\ 0& 0& M_{p_1^{j+1}}& M_{v_1^{j+1}}& M_{p_2^{j+1}}& M_{v_2^{j+1}}& 0& L& & & & M\\ M& M& O& O& O& O& O& O& & & & \\ & & L& 0& C_{p_i^{j+1}}& C_{v_i^{j+1}}& C_{p_{i+1}^{j+1}}& C_{v_{i+1}^{j+1}}& 0& L& & \\ & & L& 0& M_{p_i^{j+1}}& M_{v_i^{j+1}}& M_{p_{i+1}^{j+1}}& M_{v_{i+1}^{j+1}}& 0& L& & \\ & & & & & & O& O& O& O& 0& 0\\ & & & & & & & 0& C_{p_{N-2}^{j+1}}& C_{v_{N-2}^{j+1}}& C_{p_{N-1}^{j+1}}& C_{v_{N-1}^{j+1}}\\ 0& 0& 0& L& & & & 0& M_{p_{N-2}^{j+1}}& M_{v_{N-2}^{j+1}}& M_{p_{N-1}^{j+1}}& M_{v_{N-1}^{j+1}}\\ B_{p_0^{j+1}}& 0& 0& L& & & & & & L& 0& 0\\ 0& 0& 0& L& & & & & & L& B_{p_{N-1}^{j+1}}& 0\end{array}\right]$

经过以上“线性化”处理后，原有的非线性方程组转化成为一个在工程上等价的线性方程组，直接进行求解。

Images