基于流动参数与组分解耦的掺氢天然气管网动态仿真方法

Transient Simulation Method Based on Decoupling the Hydrothermal Parameters and Component Tracking of Hydrogen Blended Natural Gas Pipeline Network

摘要

氢能具备清洁、高效、零碳等优点，对推动能源转型具有重要意义。目前，利用现有的长输天然气管道进行掺氢输送是整合氢能大规模、长距离、高经济性输送的有效方式。在天然气管网公平开放运行模式下，天然气计量模式逐渐从体积计量转变为能量计量，而通过运行仿真进行气体组分追踪是实现天然气管网能量计量的一种有效的方法。本文分别提出了流动参数与气体组分耦合和解耦的瞬态仿真模型，所提模型中考虑了管流控制方程、对流扩散方程和节点关系方程，并采用隐式中心差分法和阻尼牛顿法对模型进行求解。耦合模型中，将组分追踪部分和流动参数计算部分作为一个整体，而解耦模型中，则将整个过程分为气体状态方程参数计算部分、流动参数计算部分和气体组分追踪部分。根据当前时步初始条件给定的不同节点气体组分来计算气体状态方程参数，并通过其计算管道流动参数，再将流动参数代入组分追踪部分来计算组分，如果迭代前后各节点组分偏差均满足设定的精度要求，则当前时步的仿真迭代结束。对某枝状管网案例进行验证，案例分析表明：耦合模型与TGNET的模拟结果中压力、温度、甲烷和氢气的摩尔分数（物质的量分数）的最大平均绝对偏差分别为0.086 MPa、0.44 K、0.10%和0.10%；在所选定的时空步长下，解耦模型计算时间比耦合模型计算时间缩短了76.04%；对于不同边界条件组合，解耦模型与TGNET的模拟结果中压力、温度、甲烷和氢气的摩尔分数的最大平均绝对偏差分别为0.078 MPa、0.44 K、0.21%和0.21%。所建立的两个模型可以达到商业软件的精度要求，能够有效地模拟水力热力参数和气体组分的时空变化，为天然气管网实际运行提供了参考；相比于耦合模型，解耦模型可以有效提高仿真效率且具有较好的边界条件适应性，能应用于实际的天然气管网运行。

Abstract

Objective Considering the cleanliness and high efficiency of hydrogen energy, the development of hydrogen energy is of great significance for raising the energy transition. Because the gas energy of hydrogen-blended natural gas varies due to differences in gas components, the measurement approach for natural gas is gradually shifting from the traditional volume-based measurement mode to the energy-based measurement mode to ensure fair natural gas transfer. Transient operation simulation that considers gas component tracking represents an effective method for achieving virtual real-time metering of natural gas components and gas energy, particularly for stations without gas chromatographs. The transient simulation results of gas components and other operating parameters provide valuable references for pipeline scheduling and risk assessment. Methods Firstly, a coupling method that considered the spatiotemporal evolution of state variables and gas components during the transient flow of gas networks was proposed to simulate the time-varying operating variables and mole fraction of each component of natural gas at each discretized node in the gas pipeline network. The pipeline governing equations consisted of the continuity equation, the momentum equation, the energy equation, and the convection-diffusion equation. The BWRS EoS and enthalpy equation were applied at each discretized node. Equations of variable relations were added at the connection nodes of multiple devices to ensure the united operation of the entire gas pipeline network within the model. The coupling model considered the component tracking part and the flow parameter calculation part as an integrated whole, whereas the decoupling model divided the entire process into the gas equation of state parameter calculation part, the flow parameter calculation part, and the gas component tracking part. The parameters of the gas equation of state were calculated based on the gas components at different nodes, which were given by the initial conditions at the current time step based on the parameter calculation equations. Based on the calculated parameters, the flow parameters were obtained by solving the transient simulation model, which consisted of the pipeline governing model, equations of variable relations at the connection nodes, boundary conditions, and initial conditions. Then, the flow parameters were substituted into the component tracking model to calculate the gas components, and this model consisted of the convection-diffusion equation, boundary conditions, and initial conditions. If the pre- and post-iteration component deviations at each node met the specified accuracy requirements, the iteration ended at the current time step; otherwise, the components were updated using the calculated values, which were then used as the initial conditions to recalculate the parameters of the gas equation of state, and the above process was repeated until the component deviation between two consecutive iterations satisfied the error criteria. Results and Discussions The simulation results of the proposed coupling model and the software TGNET were compared, and the average absolute deviations of pressure and temperature at delivery Node 3# were 0.086 MPa and 0.45 K, respectively, while the average absolute deviations of methane and hydrogen mole fractions were both 0.098%. For the decoupling model, to obtain simulation results that were independent of the discrete grid, 30 s and 60 s were selected as temporal steps, and 1 km, 2.5 km, and 5 km were selected as spatial steps for the simulations. The simulation results were compared under different combinations of temporal and spatial steps. The average absolute deviations of pressure and temperature at Node 1# were 0.072 MPa and 0.44 K, respectively, while the average absolute deviations of methane and hydrogen mole fractions were both 0.10%. The average relative deviation between the methane mole fraction calculated by the decoupling model under different combinations of spatiotemporal step sizes and the methane mole fraction simulated by commercial software was analyzed, together with the computation time required under each combination. When the temporal step size was constant, the average relative deviation decreased as the spatial step size decreased. Similarly, when the spatial step size was constant, the average relative deviation decreased as the temporal step size decreased. Considering both the computation time of the decoupling model and the average relative deviation relative to commercial software, 2.5 km and 30 s were selected as the spatial and temporal steps to simulate subsequent transient scenarios. For a simulation duration of 1 h at this spatiotemporal step size, the computation times of the coupling model and the decoupling model were 4 383 s and 1 050 s, respectively, and the computation time of the decoupling model was reduced by 76.04%. Then, the model results under three different boundary condition combinations were investigated. Under boundary condition combination A, the average absolute deviations of pressure, temperature, methane mole fraction, and hydrogen mole fraction at delivery Node 3# were 0.078 MPa, 0.39 K, 0.21%, and 0.21%, respectively, while the corresponding values at Node 1# were 0.067 MPa, 0.43 K, 0.14%, and 0.14%. Under boundary condition combination B, the average absolute deviations of pressure, temperature, methane mole fraction, and hydrogen mole fraction at delivery Node 3# were 0.078 MPa, 0.41 K, 0.21%, and 0.21%, respectively, while the corresponding values at Node 1# were 0.068 MPa, 0.44 K, 0.13%, and 0.13%. Boundary condition combination C cannot be supported by the commercial software TGNET. Conclusions The proposed coupling and decoupling models demonstrate high accuracy in flow parameter calculation and gas component tracking, reaching a performance level comparable to that of commercial software. These models are not only applicable to linear natural gas pipelines but also exhibit high accuracy in simulating multi-branch natural gas pipeline systems. At a specific spatiotemporal step size, the computation time of the decoupling model is 76.04% shorter than that of the coupling model, indicating superior simulation efficiency. In addition, the decoupling model can accurately capture the spatiotemporal evolution patterns of state variables and gas components under various combinations of boundary conditions, demonstrating strong adaptability to different boundary condition scenarios.

Graphical abstract

关键词

掺氢天然气管网 / 瞬态模拟 / 组分追踪 / 解耦模型 / 耦合模型

Key words

hydrogen blended natural gas pipeline networks / transient simulation / components tracking / decoupling model / coupling model

本刊网刊

随着中国“双碳”战略目标的提出，绿色低碳能源的地位愈加显著^[1]。氢能作为清洁、高效、零碳的能源载体，被视为21世纪最具发展潜力的清洁能源，大力发展氢能对推动能源转型具有重要意义^[2]。目前，利用现有的长输天然气管道进行掺氢输送^[3]，可以极大地降低氢气长输管道建设成本，是整合氢能大规模、长距离、高经济性输送的有效方式^[4]。考虑到一般天然气的性质差异不大，传统的天然气管网仿真针对水力和热力计算进行研究^[5‒7]，得到管网系统的压力、温度和流量等水热力参数的时空演化规律^[8‒9]，以确保管网的安全运行^[10‒11]。然而，氢气的加入会改变天然气的热值^[12]，为了保证天然气贸易的公平，天然气的计量方式正逐渐由传统的体积计量模式转变为以天然气气质为基础的能量计量模式^[13‒14]，因此，还有必要对掺氢后的天然气进行组分追踪^[15‒17]。掺氢天然气瞬态组分追踪模拟是实现掺氢天然气组分和气体能量虚拟实时计量的有效方法，尤其适用于没有气相色谱仪的站场^[18]，其瞬态仿真结果和其他运行参数可为管道调度和风险评估提供参考。

在掺氢天然气管道的实际运行中，运行状态受到各种因素的影响，如交付节点或注入节点的气体流量波动、装置运行方案的故障或调整等。管道的运行参数，包括压力、温度、流量、密度和气体组分等，随时间和空间的变化而变化，所以，天然气管道的流动状态通常是瞬态的^[19‒20]。许多学者提出了瞬态流模拟模型^[21‒23]，所提出的瞬态仿真模型主要包括管流控制方程、压缩机和阀门的运行特性、初始条件和一些边界条件。特征线和中心差分法是求解天然气管网瞬态仿真模型的主要方法^[24‒25]。针对瞬态气体组分追踪的需要，考虑气体组分与运行参数之间的耦合关系，学者们提出了不同的基于组分追踪的瞬态运行仿真模型。对于单条管道，Chaczykowski等^[26]提出了一种瞬态气体组分追踪模型。对于复杂管网，Zhang等^[27]提出了基于气体组分追踪的等温瞬态流仿真模型。考虑到实际管输过程中的大部分工况为非等温工况，在构建组分追踪仿真模型时有必要考虑管道温度的时空变化情况。Bermúdez等^[28]假设热力参数为已知变量，提出了一种水力方程和组分方程相解耦的方法，分别求解每个时间步的流动方程和对流扩散方程。Fan等^[29]考虑热力参数为未知变量，提出了一种解耦瞬态气体组分追踪模型，先对水力参数进行迭代修正，再根据修正后的水力参数直接计算热力参数和组分参数。由于气体流动参数和组分参数的数量多且不同参数之间非线性耦合，所以现有的基于组分追踪的管网运行仿真求解方法存在求解精度低或求解时间长的问题。

为此，本文分别提出了流动参数与气体组分耦合和解耦的瞬态仿真模型。在所构建的基于流动参数和气体组分耦合模型的基础上，提出了流动参数和气体组分相解耦的仿真方法，该方法在保证计算精度的同时，能够显著缩短仿真计算时间。同时，为了适应实际管网运行过程中边界条件组合的变化，本文在模型验证部分采用了不同边界条件组合来分析模型的边界条件适应性。

1 流动参数和组分耦合模型

构建了一种考虑管网瞬态流动过程中流动参数和气体组分时空演变规律的耦合模型，用于预测掺氢天然气管网中每个离散节点处的流动参数和掺氢天然气组分的时变规律。天然气管网瞬态仿真耦合模型框架如图1所示，管流控制方程由连续性方程、动量方程和能量方程组成，使用对流扩散方程来实现对管网中气体的追踪，将BWRS（Benedict-Webb-Rubin-Starling）气体状态方程和焓方程用于每个离散化节点。同时，在多个管段的连接节点处考虑变量关系方程，保证了整个管网是统一运行的。通过上述方程、边界条件和初始条件可以建立耦合模型，隐式中心差分法被用于对管段微分方程进行离散化，采用阻尼牛顿法求解模型，得到每个时刻的压力、温度、质量流量、密度、焓和组分的分布。

1.1 1　管道流动偏微分方程组

在流动参数计算中，气体的压力、温度和流量是由连续性方程、动量方程和能量方程连接的。以上方程被用于描述流动参数的时空变化规律，连续性方程、动量方程和能量方程的分别如下：

\begin{array}{c}\frac{\partial \rho (x,\tau )}{\partial \tau }+\frac{\mathrm{1}}{A}\cdot \frac{\partial M(x,\tau )}{\partial x}=\mathrm{0}\end{array}（1）

\begin{array}{c}\frac{\mathrm{1}}{A}\cdot \frac{\partial M(x,\tau )}{\partial \tau }+\frac{\partial \left(p(x,\tau )+\frac{{\left(M(x,\tau )\right)}^{\mathrm{2}}}{A^{\mathrm{2}}\rho (x,\tau )}\right)}{\partial x}=\\ -\frac{\lambda }{\mathrm{2}d}\cdot \frac{{\left(M(x,\tau )\right)}^{\mathrm{2}}}{A^{\mathrm{2}}}\cdot \frac{\mathrm{1}}{\rho (x,\tau )}-\rho (x,\tau )g\frac{\mathrm{d}s}{\mathrm{d}x}\end{array}（2）

\begin{array}{c}\frac{\partial }{\partial \tau }\left[\left(h(x,\tau )-\frac{p(x,\tau )}{\rho (x,\tau )}+\frac{{\left(M(x,\tau )\right)}^{\mathrm{2}}}{\mathrm{2}{A}^{\mathrm{2}}{\left(\rho (x,\tau )\right)}^{\mathrm{2}}}\right)\rho (x,\tau )\right]+\\ \frac{\partial }{\partial x}\left[\left(h(x,\tau )+\frac{M{(x,\tau )}^{\mathrm{2}}}{\mathrm{2}{A}^{\mathrm{2}}{\left(\rho (x,\tau )\right)}^{\mathrm{2}}}\right)\frac{M(x,\tau )}{A}\right]+\\ \frac{\mathrm{4}K\left(T(x,\tau )-T_{\mathrm{0}}\right)}{D}+\frac{M(x,\tau )}{A}g\frac{\mathrm{d}s}{\mathrm{d}x}=\mathrm{0}\end{array}（3）

式（1）～（3）中：\rho为气体密度，kg/m³；M为气体质量流量，kg/s；A为管段横截面积，m²；p为气体压力，Pa；\lambda为水力摩阻系数；d为管段内径，m；g为重力加速度，m/s²；s为管段的高程；h为气体的焓，J/kg；D为管段外径，m；K为总传热系数，W/(m²\cdotK)；T为气体温度，K；T₀为环境温度，K；\tau为模拟时间，s；x为管段长度，m。

为了准确计算水力摩阻系数，本文采用了Colebrook-White方程：

\begin{array}{c}\frac{\mathrm{1}}{\sqrt[]{\lambda }}=-\mathrm{2}\mathrm{l}\mathrm{g}\left(\frac{\epsilon }{\mathrm{3}\mathrm{700}d}+\frac{\mathrm{2.51}}{Re\sqrt[]{\lambda }}\right)\end{array}（4）

式中：\epsilon为管段粗糙度，mm；\mu为气体黏度，\mathrm{P}\mathrm{a}\cdot \mathrm{s}；Re为雷诺数，计算公式为：

\begin{array}{c}Re=\frac{Md}{\mu A}\end{array}（5）

同时，为了预测天然气组分的变化规律，本文引入了对流扩散方程来描述天然气各组分的质量守恒：

\begin{array}{c}\frac{\partial {\rho }_n(x,\tau )}{\partial \tau }+\frac{\mathrm{1}}{A}\cdot \frac{\partial M_n(x,\tau )}{\partial x}=\mathrm{0}\end{array}（6）

式中：{\rho }_n为天然气第n个组分的密度，kg/m³；M_n 为天然气第n个组分的质量流量，kg/s。

1.1 2　控制方程的离散

考虑到求解以上偏微分方程组的解析解较为困难，本文采用了具有高稳定性的隐式中心差分法来对偏微分方程进行离散，如图2所示，从而获得对应的数值解。

通过对连续性方程、动量方程、能量方程和对流扩散方程进行离散，可以得到以下离散后的连续性方程、动量方程、能量方程和对流扩散方程：

\begin{array}{l}{\rho }_{i+\mathrm{1}}^{k+\mathrm{1}}-{\rho }_{i+\mathrm{1}}^k+{\rho }_i^{k+\mathrm{1}}-{\rho }_i^k+\\ \frac{\mathrm{1}}{A}\frac{\mathrm{\Delta }\tau }{\mathrm{\Delta }x}\left(M_{i+\mathrm{1}}^{k+\mathrm{1}}-M_i^{k+\mathrm{1}}+M_{i+\mathrm{1}}^k-M_i^k\right)=\mathrm{0}\end{array}（7）

\begin{array}{l}\frac{\mathrm{1}}{A}(M_{i+\mathrm{1}}^{k+\mathrm{1}}+M_i^{k+\mathrm{1}}-M_{i+\mathrm{1}}^k-M_i^k)+\\ \frac{\mathrm{\Delta }\tau }{\mathrm{\Delta }x}\left[{\left(\frac{M_{i+\mathrm{1}}^{k+\mathrm{1}}}{A}\right)}^{\mathrm{2}}\cdot \frac{\mathrm{1}}{{\rho }_{i+\mathrm{1}}^{k+\mathrm{1}}}+p_{i+\mathrm{1}}^{k+\mathrm{1}}-{\left(\frac{M_i^{k+\mathrm{1}}}{A}\right)}^{\mathrm{2}}\cdot \frac{\mathrm{1}}{{\rho }_i^{k+\mathrm{1}}}-\right.\\ \left.p_i^{k+\mathrm{1}}+{\left(\frac{M_{i+\mathrm{1}}^k}{A}\right)}^{\mathrm{2}}\cdot \frac{\mathrm{1}}{{\rho }_{i+\mathrm{1}}^k}+p_{i+\mathrm{1}}^k-{\left(\frac{M_i^k}{A}\right)}^{\mathrm{2}}\cdot \frac{\mathrm{1}}{{\rho }_i^k}-p_i^k\right]+\\ \frac{\mathrm{\Delta }\tau \lambda }{\mathrm{4}d}\left[{\left(\frac{M_{i+\mathrm{1}}^{k+\mathrm{1}}}{A}\right)}^{\mathrm{2}}\cdot \frac{\mathrm{1}}{{\rho }_{i+\mathrm{1}}^{k+\mathrm{1}}}+{\left(\frac{M_i^{k+\mathrm{1}}}{A}\right)}^{\mathrm{2}}\cdot \frac{\mathrm{1}}{{\rho }_i^{k+\mathrm{1}}}+\right.\\ \left.{\left(\frac{M_{i+\mathrm{1}}^k}{A}\right)}^{\mathrm{2}}\cdot \frac{\mathrm{1}}{{\rho }_{i+\mathrm{1}}^k}+{\left(\frac{M_i^k}{A}\right)}^{\mathrm{2}}\cdot \frac{\mathrm{1}}{{\rho }_i^k}\right]+\\ \mathrm{\Delta }\tau g\frac{\mathrm{d}s}{\mathrm{d}x}\cdot \frac{\left({\rho }_{i+\mathrm{1}}^{k+\mathrm{1}}+{\rho }_{i+\mathrm{1}}^k+{\rho }_i^{k+\mathrm{1}}+{\rho }_i^k\right)}{\mathrm{2}}=\mathrm{0}\end{array}（8）

\begin{array}{c}{h}_i^{k+\mathrm{1}}{\rho }_i^{k+\mathrm{1}}-p_i^{k+\mathrm{1}}+{\left(\frac{M_i^{k+\mathrm{1}}}{A}\right)}^{\mathrm{2}}\cdot \frac{\mathrm{1}}{\mathrm{2}{\rho }_i^{k+\mathrm{1}}}-h_i^k{\rho }_i^k+\\ p_i^k-{\left(\frac{M_i^k}{A}\right)}^{\mathrm{2}}\cdot \frac{\mathrm{1}}{\mathrm{2}{\rho }_i^{k+\mathrm{1}}}+h_{i+\mathrm{1}}^{k+\mathrm{1}}{\rho }_{i+\mathrm{1}}^{k+\mathrm{1}}-\\ p_{i+\mathrm{1}}^{k+\mathrm{1}}+{\left(\frac{M_{i+\mathrm{1}}^{k+\mathrm{1}}}{A}\right)}^{\mathrm{2}}\cdot \frac{\mathrm{1}}{\mathrm{2}{\rho }_{i+\mathrm{1}}^{k+\mathrm{1}}}-h_{i+\mathrm{1}}^k{\rho }_{i+\mathrm{1}}^k+\\ p_{i+\mathrm{1}}^k-{\left(\frac{M_{i+\mathrm{1}}^k}{A}\right)}^{\mathrm{2}}\cdot \frac{\mathrm{1}}{\mathrm{2}{\rho }_{i+\mathrm{1}}^k}+\\ \frac{\mathrm{\Delta }\tau }{\mathrm{\Delta }x}\left[h_{i+\mathrm{1}}^{k+\mathrm{1}}\frac{M_{i+\mathrm{1}}^{k+\mathrm{1}}}{A}+{\left(\frac{M_{i+\mathrm{1}}^{k+\mathrm{1}}}{A}\right)}^{\mathrm{3}}\cdot \frac{\mathrm{1}}{\mathrm{2}({\rho }_{i+\mathrm{1}}^{k+\mathrm{1}}{)}^{\mathrm{2}}}-\right.\\ h_i^{k+\mathrm{1}}\frac{M_i^{k+\mathrm{1}}}{A}-{\left(\frac{M_i^{k+\mathrm{1}}}{A}\right)}^{\mathrm{3}}\cdot \frac{\mathrm{1}}{\mathrm{2}({\rho }_i^{k+\mathrm{1}}{)}^{\mathrm{2}}}+\\ h_{i+\mathrm{1}}^k\frac{M_{i+\mathrm{1}}^k}{A}+{\left(\frac{M_{i+\mathrm{1}}^k}{A}\right)}^{\mathrm{3}}\cdot \frac{\mathrm{1}}{\mathrm{2}({\rho }_{i+\mathrm{1}}^k{)}^{\mathrm{2}}}-\\ \left.h_i^k\frac{M_i^k}{A}-{\left(\frac{M_i^k}{A}\right)}^{\mathrm{3}}\cdot \frac{\mathrm{1}}{\mathrm{2}({\rho }_i^k{)}^{\mathrm{2}}}\right]+\\ \mathrm{2}\frac{K\mathrm{\Delta }\tau }{D}(T_{i+\mathrm{1}}^{k+\mathrm{1}}+T_{i+\mathrm{1}}^k+T_i^{k+\mathrm{1}}+T_i^k-\mathrm{4}{T}_{\mathrm{0}})+\\ \frac{g\mathrm{\Delta }\tau }{\mathrm{2}}\cdot \frac{\mathrm{d}s}{\mathrm{d}x}\left(\frac{M_{i+\mathrm{1}}^{k+\mathrm{1}}}{A}+\frac{M_{i+\mathrm{1}}^k}{A}+\frac{M_i^{k+\mathrm{1}}}{A}+\frac{M_i^k}{A}\right)=\mathrm{0}\end{array}（9）

\begin{array}{l}{\rho }_{n,i+\mathrm{1}}^{k+\mathrm{1}}-{\rho }_{n,i+\mathrm{1}}^k+{\rho }_{n,i}^{k+\mathrm{1}}-{\rho }_{n,i}^k+\\ \frac{\mathrm{1}}{A}\cdot \frac{\mathrm{\Delta }\tau }{\mathrm{\Delta }x}\left(M_{n,i+\mathrm{1}}^{k+\mathrm{1}}+M_{n,i+\mathrm{1}}^k-M_{n,i}^{k+\mathrm{1}}-M_{n,i}^k\right)=\mathrm{0}\end{array}（10）

式（7）～（10）中，每个变量的上标k和k+1表示两个相邻的时间步长，每个变量的下标i和i+1表示两个相邻的空间步长，\mathrm{\Delta }\tau为时间步长，\mathrm{\Delta }x为空间步长。

1.2 气体物性参数计算

在管流控制方程中，气体物性参数是重要的组成部分。因此，为了满足数学模型的封闭性，在模型中需要加入对应的方程。气体压力、温度、组分和气体状态方程是计算气体物性参数的基础。

1.2 1　气体状态方程

为了准确计算控制方程中涉及的运行参数的变化，本文采用了被认为是最精确的天然气状态方程之一的BWRS状态方程（式（11））^[30‒32]来反映气体压力、温度、密度和组分之间的关系。

\begin{array}{c}p=\left[{\rho }_{m}RT+\left(B_{\mathrm{0}}RT-A_{\mathrm{0}}-\frac{C_{\mathrm{0}}}{T^{\mathrm{2}}}+\frac{D_{\mathrm{0}}}{T^{\mathrm{3}}}-\frac{E_{\mathrm{0}}}{T^{\mathrm{4}}}\right){\rho }_m^{\mathrm{2}}+\right.\\ \left(bRT-a-\frac{d}{T}\right){\rho }_{\mathrm{m}}^{\mathrm{3}}+\alpha \left(a+\frac{d}{T}\right){\rho }_{\mathrm{m}}^{\mathrm{6}}+\\ \left.\frac{c{\rho }_{\mathrm{m}}^{\mathrm{3}}}{T^{\mathrm{2}}}\left(\mathrm{1}+\gamma {\rho }^{\mathrm{2}}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(-\gamma {\rho }_{\mathrm{m}}^{\mathrm{2}}\right)\right]\cdot \mathrm{1}\mathrm{000}\end{array}（11）

式中：{\rho }_{\mathrm{m}}为气体摩尔密度，kmol/m³；A₀、B₀、C₀、D₀、E₀、a、b、c、d、\gamma和\alpha是气体混合物的11个特征参数。

对于纯物质，A_0,i 、B_0,i 、C_0,i 、D_0,i 、E_0,i 、a_i 、b_i 、c_i 、d_i 、{\gamma }_i和{\alpha }_i的数值可以从数据库中得到。对于由多个纯物质组成的气体混合气体，A₀、B₀、C₀、D₀、E₀、a、b、c、d、\gamma和\alpha可以通过各纯物质的特征参数和特定的混合规则计算得到，具体的混合规则分别如下：

\begin{array}{c}\begin{array}{c}{A}_{\mathrm{0}}=\end{array}\end{array}{\displaystyle \sum _{i=\mathrm{1}}^k}{\displaystyle \sum _{j=\mathrm{1}}^k}{y}_i{y}_j{A}_{\mathrm{0}i}^{\frac{\mathrm{1}}{\mathrm{2}}}{A}_{\mathrm{0}j}^{\frac{\mathrm{1}}{\mathrm{2}}}(\mathrm{1}-k_{ij})（12）

\begin{array}{c}\begin{array}{c}{B}_{\mathrm{0}}={\displaystyle \sum _{i=\mathrm{1}}^k}{y}_i{B}_{\mathrm{0}i}\end{array}\end{array}（13）

\begin{array}{c}\begin{array}{c}{C}_{\mathrm{0}}={\displaystyle \sum _{i=\mathrm{1}}^k}{\displaystyle \sum _{j=\mathrm{1}}^k}{y}_i{y}_j{C}_{\mathrm{0}i}^{\frac{\mathrm{1}}{\mathrm{2}}}{C}_{\mathrm{0}j}^{\frac{\mathrm{1}}{\mathrm{2}}}(\mathrm{1}-k_{ij}{)}^{\mathrm{3}}\end{array}\end{array}（14）

\begin{array}{c}\begin{array}{c}{D}_{\mathrm{0}}={\displaystyle \sum _{i=\mathrm{1}}^k}{\displaystyle \sum _{j=\mathrm{1}}^k}{y}_i{y}_j{D}_{\mathrm{0}i}^{\frac{\mathrm{1}}{\mathrm{2}}}{D}_{\mathrm{0}j}^{\frac{\mathrm{1}}{\mathrm{2}}}(\mathrm{1}-k_{ij}{)}^{\mathrm{4}}\end{array}\end{array}（15）

\begin{array}{c}\begin{array}{c}{E}_{\mathrm{0}}={\displaystyle \sum _{i=\mathrm{1}}^k}{\displaystyle \sum _{j=\mathrm{1}}^k}{y}_i{y}_j{E}_{\mathrm{0}i}^{\frac{\mathrm{1}}{\mathrm{2}}}{E}_{\mathrm{0}j}^{\frac{\mathrm{1}}{\mathrm{2}}}(\mathrm{1}-k_{ij}{)}^{\mathrm{5}}\end{array}\end{array}（16）

\begin{array}{c}\begin{array}{c}a={\left({\displaystyle \sum _{i=\mathrm{1}}^k}{y}_i{a}_i^{\frac{\mathrm{1}}{\mathrm{3}}}\right)}^{\mathrm{3}}\end{array}\end{array}（17）

{\begin{array}{c}\begin{array}{c}b=\left({\displaystyle \sum _{i=\mathrm{1}}^k}{y}_i{b}_i^{\frac{\mathrm{1}}{\mathrm{3}}}\right)\end{array}\end{array}}^{\mathrm{3}}（18）

\begin{array}{c}\begin{array}{c}c{\begin{array}{c}\begin{array}{c}=\left({\displaystyle \sum _{i=\mathrm{1}}^k}{y}_i{c}_i^{\frac{\mathrm{1}}{\mathrm{3}}}\right)\end{array}\end{array}}^{\mathrm{3}}{}^{}\end{array}\end{array}（19）

d{\begin{array}{c}\begin{array}{c}=\left({\displaystyle \sum _{i=\mathrm{1}}^k}{y}_i{d}_i^{\frac{\mathrm{1}}{\mathrm{3}}}\right)\end{array}\end{array}}^{\mathrm{3}}（20）

\begin{array}{c}\begin{array}{c}\gamma =\end{array}\end{array}{\left(y_i{y}_i^{\frac{\mathrm{1}}{\mathrm{2}}}\right)}^{\mathrm{2}}（21）

\begin{array}{c}\begin{array}{c}\alpha =\end{array}\end{array}{\left({{\displaystyle \sum _{i=\mathrm{1}}^k}{y}_i\alpha }_i^{\frac{\mathrm{1}}{\mathrm{3}}}\right)}^{\mathrm{3}}（22）

式（12）～（22）中：y_i 、y_j 分别为气体的第i和j的组分的摩尔分数（物质的量分数），k_ij 为气体的第i和j的组分的交互系数，k为混合物气体的组分数。

1.2 2　相关气体物性参数

在热力计算中，焓是能量方程中非常重要的一个参数。在掺氢天然气管网组分追踪过程中，气体的焓总是受气体组分、压力和温度影响。通过BWRS气体状态方程和热力学微分方程可以推导出实际的混合气体焓方程，将理想混合气体的焓值代入到以上推导的方程中可以得到实际混合气体的焓值：

\begin{array}{l}h=\mathrm{1}\mathrm{000}\cdot H^{\mathrm{0}}+\\ \left(\left(B_{\mathrm{0}}RT-\mathrm{2}{A}_{\mathrm{0}}-\frac{\mathrm{4}{C}_{\mathrm{0}}}{T^{\mathrm{2}}}+\frac{\mathrm{5}{D}_{\mathrm{0}}}{T^{\mathrm{3}}}-\frac{\mathrm{6}{E}_{\mathrm{0}}}{T^{\mathrm{4}}}\right)\right.{\rho }_{\mathrm{m}}+\\ \frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{2}bRT-\mathrm{3}a-\frac{\mathrm{4}d}{T}\right){\rho }_{\mathrm{m}}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{5}}\alpha \left(\mathrm{6}a+\frac{\mathrm{7}d}{T}\right){\rho }_{\mathrm{m}}^{\mathrm{5}}+\\ \left.\frac{c}{\gamma T^{\mathrm{2}}}\left(\mathrm{3}-\left(\mathrm{3}+\frac{\gamma {\rho }_{\mathrm{m}}^{\mathrm{2}}}{\mathrm{2}}-{\gamma }^{\mathrm{2}}{\rho }_{\mathrm{m}}^{\mathrm{4}}\right)\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(-\gamma {\rho }_{\mathrm{m}}^{\mathrm{2}}\right)\right)\cdot \frac{\mathrm{1}\mathrm{000}}{M_{\mathrm{g}}}\end{array}（23）

式中：H⁰为理想混合气体的焓，kJ/kg；M_g为气体相对分子质量，g/mol。

1.3 连接节点处变量方程

在天然气管网仿真中，存在许多连接节点，如管段与压缩机之间的连接节点、管段与阀门之间的连接节点、多个管段之间的连接节点。每个连接节点上都有变量关系，为了准确模拟节点连接处的流动状态，应在连接节点周围设置一些虚拟节点以反映这些特殊节点的流入和流出参数。对于每个多管段连接节点（图3），有压力关系方程（式（24））、流量平衡方程（式（25））、能量平衡方程（式（26））和天然气各组分的流量平衡方程（式（27））。

\frac{\partial p_{\mathrm{J}\mathrm{u}\mathrm{n}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}}{\partial \tau }=\frac{\partial p_j^{\mathrm{o}\mathrm{u}{\mathrm{t}}^{\text{'}}}}{\partial \tau },j=\mathrm{1,2},\cdots ,n_{\mathrm{1}}（24）

\frac{\partial M_{\mathrm{J}\mathrm{u}\mathrm{n}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}}{\partial \tau }={\displaystyle \sum _{j=\mathrm{1}}^{n_{\mathrm{1}}}}\frac{\partial M_j^{\mathrm{o}\mathrm{u}{\mathrm{t}}^{\text{'}}}}{\partial \tau }（25）

\frac{\partial {(TM)}_{\mathrm{J}\mathrm{u}\mathrm{n}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}}{\partial \tau }={\displaystyle \sum _{j=\mathrm{1}}^{n_{\mathrm{1}}}}\frac{{(TM)}_j^{\mathrm{o}\mathrm{u}{\mathrm{t}}^{\text{'}}}}{\partial \tau }（26）

\frac{\partial (x_n{M)}_{\mathrm{J}\mathrm{u}\mathrm{n}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}}}{\partial \tau }={\displaystyle \sum _{j=\mathrm{1}}^{n_{\mathrm{1}}}}\frac{\partial (x_n{M)}_j^{\mathrm{o}\mathrm{u}{\mathrm{t}}^{\text{'}}}}{\partial \tau },n=\mathrm{1,2},\cdots ,k（27）

式（24）～（27）中，下标Junction表示节点处于连接节点处，上标\mathrm{o}\mathrm{u}{\mathrm{t}}^{\text{'}}表示节点处于出口节点处，下标j表示节点处于第j个出口节点处，x_n 为气体第n个组分的质量分数。

1.4 边界条件

边界条件是确定仿真模型唯一解的重要条件，表示边界节点在整个时域的运行状态。天然气管网的瞬态仿真分析中，对各供气节点选择3个边界条件，给出气源压力、温度、组分的时变函数，或给出气源的温度、流量和组分的时变函数作为边界条件组合。在每个交付节点上，选择质量流量或压力随时间变化的函数作为边界条件。

1.5 初始条件

在建立瞬态仿真模型时，初始条件是必不可少的。在天然气管网运行中，初始条件可以通过数据采集与监视控制（SCADA）系统、瞬态预仿真或稳态仿真获得。本文通过初始时间的边界条件进行稳态模拟，得到了初始条件。从离散方程可以看出，时间τ+Δτ的运行参数仅取决于时间τ的运行参数，因此，采用时间递推方法进行仿真。对于第1个时间步后的仿真，初始条件可以由前1个时间步的仿真结果提供。

1.6 模型求解

通过上述分析可以看出，耦合模型考虑了运行过程中水力变量、热力变量和组分参数的非线性耦合关系，这意味着需要对大规模非线性方程组进行求解才可以得到管网的运行参数。求解非线性方程组的方法一般为传统的牛顿法，但不能保证牛顿法的迭代方向是一定沿着函数值下降的方向，因此，本文采用阻尼牛顿法作为求解上述非线性方程组的方法。该方法能够保证每一次的迭代是朝着函数值下降的方向移动，同时，可以保持传统牛顿法的优点，即快速收敛和高度近似到最优。阻尼牛顿法模拟流程如图4所示。图4中，\mathit{\delta }为迭代步长，\mathit{\beta }为搜索方向。在求解非线性方程组 F ( x )=0之前，应先设定迭代初值 x₀和迭代次数r为0；再计算梯度{\mathit{g}}_r和Hessian矩阵 H_r；判断迭代收敛要求‖{\mathit{g}}_r‖\le \xi是否被满足，如果不满足，则应更新 x_r+1和r，程序在此时间步长进入下一次迭代；否则，当前时间步长的模拟将停止，程序进入下一个时间步长的模拟。

需要注意的是，该方法同样被用于求解解耦模型中的非线性方程组，以上求解时梯度计算方法、Hessian矩阵计算方法如式（28）和（29）所示：

\mathit{g}(\mathit{x})=\nabla \mathit{F}(\mathit{x})={\left(\frac{\partial \mathit{F}}{\partial x_{\mathrm{1}}},\frac{\partial \mathit{F}}{\partial x_{\mathrm{2}}},\cdots ,\frac{\partial \mathit{F}}{\partial x_r}\right)}^{\mathrm{T}}（28）

\mathit{H}(\mathit{x})=\left[\begin{array}{cccc}\frac{{\partial }^{\mathrm{2}}\mathit{F}}{\partial x_{\mathrm{1}}\partial x_{\mathrm{1}}}& \frac{{\partial }^{\mathrm{2}}\mathit{F}}{\partial x_{\mathrm{1}}\partial x_{\mathrm{2}}}& \cdots & \frac{{\partial }^{\mathrm{2}}\mathit{F}}{\partial x_{\mathrm{1}}\partial x_r}\\ \frac{{\partial }^{\mathrm{2}}\mathit{F}}{\partial x_{\mathrm{2}}\partial x_{\mathrm{1}}}& \frac{{\partial }^{\mathrm{2}}\mathit{F}}{\partial x_{\mathrm{2}}\partial x_{\mathrm{2}}}& \cdots & \frac{{\partial }^{\mathrm{2}}\mathit{F}}{\partial x_{\mathrm{2}}\partial x_r}\\ ⋮& ⋮& & ⋮\\ \frac{{\partial }^{\mathrm{2}}\mathit{F}}{\partial x_r\partial x_{\mathrm{1}}}& \frac{{\partial }^{\mathrm{2}}\mathit{F}}{\partial x_r\partial x_{\mathrm{2}}}& \cdots & \frac{{\partial }^{\mathrm{2}}\mathit{F}}{\partial x_r\partial x_r}\end{array}\right]（29）

2 流动参数和组分解耦模型

考虑到上述瞬态组分追踪耦合模型在管网规模较大时的仿真效率相对较低，本文通过对组分追踪部分和水力热力参数计算部分分别解耦，建立了一种新的仿真模型。所提出的解耦模型在保证仿真结果计算精度的同时，可以有效缩短仿真模型的求解时间。

2.1 解耦模型的框架

图5为耦合模型与解耦模型的关系。如图5所示，解耦模型由3个部分组成。

首先，给定初始条件，应用气体状态方程混合规则建立气体状态方程参数计算子模型。然后，结合连续性方程、动量方程、能量方程、连接节点方程、焓方程、BWRS气体状态方程、边界条件和初始条件可以建立水力热力子模型。求解水力热力子模型后，可以得到一个时步内每个节点的密度和质量流量。最后，结合水力热力子模型所提供参数，对于对流扩散方程、连接节点方程、边界条件和初始条件，可以建立组分追踪子模型。

2.2 求解框架

分析解耦仿真模型的计算步骤如图6所示。具体步骤如下：

1）在每一个时步内，给定各节点气体组分的质量分数，可计算出A₀、B₀、C₀、D₀、E₀、a、b、c、d、\gamma和\alpha等参数。

2）通过式（7）～（9）、（11）、（23）、（24）～（27），边界条件和初始条件可以建立水力热力子模型，对该非线性方程组求解后可以得到各节点的气体压力、温度、质量流量、密度和焓。

3）将经由求解水力热力子模型得到的密度和质量流量提供给组分追踪子模型，通过式（10）、（24）～（27）、结合边界条件和初始条件，可以建立并求解组分追踪模型，得到各节点的气体组分。

4）如果两次迭代前后各节点的组分质量分数偏差值大于误差偏差，则每个节点的组分质量分数使用组分追踪过程计算的质量分数进行更新；否则，水力热力参数计算过程和组分追踪过程的迭代将停止，并进入下一个时间步的模拟。

5）当当前仿真时间大于或等于结束仿真时间时，仿真结束。

本文解耦仿真模型与TGNET软件算法相比，有以下不同之处：1）TGNET通过线性化简化了水力参数、热力参数、组分参数和设备运行参数之间的耦合关系；在所提解耦仿真模型中，保留了该非线性耦合关系。2）在每一个时步的仿真过程中，所提解耦仿真模型根据连续两次迭代的组分偏差判断迭代过程是否满足收敛条件，若满足，则该时步仿真结束；否则，以新计算的组分进行下一次迭代。TGNET直接计算组分参数，未作回代修正。3）在动态仿真时，TGNET的供气节点和交付节点的边界条件选择存在限制，必须至少包含一个压力边界才能使模型被正确求解；所提解耦仿真模型不存在该限制。

所提解耦仿真模型有如下优点：1）考虑了水力变量、热力变量、组分变量和设备运行参数的非线性耦合关系，提高了模拟结果的准确性。2）通过对连续两次迭代的组分偏差进行误差判定，使用回代新组分进行修正的方法提高了仿真结果的准确性。3）对边界条件有较好的适应性，可适用于不同的运行条件；在一些极端运行条件下，该模型仍然有效。

3 结果与讨论

通过以下的案例验证了所提出掺氢天然气管网瞬态组分追踪模型的有效性。

图7为掺氢天然气管网系统的简单拓扑图。管段基本参数，如表1所示。

管道系统的整体传热系数为1.15 W/(m²∙K)。供气节点1#的边界条件为气体压力、温度和气体组分，供气节点2#、3#的边界条件为气体压力、质量流量和气体组分。选择随时间变化的气体质量流量函数作为所有交付节点的边界条件。供气节点的边界条件、供气节点和交付节点的瞬态边界条件及随时间的变化如表2和图8所示。以下瞬态仿真时间为1 h。

3.1 封闭性分析

耦合模型的封闭性分析如下：假设管段i被划分成\Psi_i 个子管段，即含\Psi_i +1个节点，则所有管段被分为\Omega（即{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{9}}}{\Psi }_i）个子管段，包含\Omega+9个节点。该仿真模型中气体压力、温度、质量流量、密度、焓和每个组分的质量分数为变量（假设气体混合物由k个纯物质组成），那么，在每个时步中该模型有\left(\mathrm{5}+k\right)\times \left(\Omega +\mathrm{9}\right)个变量。每个管段处，有3\Omega个基础管流控制方程和k\Omega个气体组分质量守恒方程，即共\left(\mathrm{3}+k\right)\Omega个方程。每个节点处都有气体状态方程和焓方程，即\mathrm{2}\left(\Omega +\mathrm{9}\right)个方程。每个多管段的混合流入节点处有k+4个方程，每个多管段的混合流出节点处有5+2k个方程。分析该管网可知，共计2个多管段的混合流入节点和2个多管段混合的流出节点，因此，多管段混合节点处有\mathrm{2}\left(\mathrm{9}+\mathrm{3}k\right)个方程。此外，每个供气节点都有k+2个边界条件，每个交付节点都有1个边界条件，总共\mathrm{3}\left(k+\mathrm{2}\right)+\mathrm{3}个边界条件。上述情况下，仿真模型中的方程数等于每个时间步长下的变量数，因此，耦合模型是封闭的。

解耦模型的封闭性分析如下：假设管网子管段和节点划分与耦合模型一样，所有管段被分为\Omega（即{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{9}}}{\Psi }_i）个子管段、\Omega+9个节点。在水力热力子模型中，气体压力、温度、密度、质量流量和焓均为变量，那么在每个时步中该子模型有\mathrm{5}\left(\Omega +\mathrm{9}\right)个变量。每个管段处有3个基础管流控制方程，每个节点处有气体状态方程和焓方程，则共计\mathrm{3}\Omega +\mathrm{2}\left(\Omega +\mathrm{9}\right)个方程。每个多管段的混合流入节点处有4个方程，每个多管段的混合流出节点处有5个方程。分析该管网可知，共有2个多管段的混合流入节点和2个多管段的混合流出节点，所以，多管段混合节点处共18个方程。此外，供气节点1#的边界条件为气体压力和温度，供气节点2#、3#的边界条件为气体质量流量和温度，交付节点采用气体质量流量作为边界条件，共有9个方程。在水力热力子模型中，每个时步的方程数等于变量数，因此，水力热力子模型是封闭的。

在组分追踪子模型中，气体组分的质量分数为变量，那么，每个时步中该子模型有k\left(\Omega +\mathrm{9}\right)个变量。每个管段处有k个组分的质量守恒方程，则共有k\Omega个方程。每个多管段的混合流入节点处有k个方程，每个多管段的混合流出节点处有2k个方程。分析管网可知，共2个多管段混合流入节点和2个多管段混合流出节点，所以，多管段的混合节点处共6k个方程。此外，给出各供气节点的气体组分，设为供气节点的边界条件，即有3k个方程。因此，组分追踪子模型中每个时步的方程数等于变量数，即组分追踪子模型也是封闭的。综上，解耦模型是封闭的。

3.2 结果分析

在构建耦合模型时，选择1 km和60 s分别作为仿真的空间步长和时间步长。节点#1和交付节点#3处本文耦合模型和商业软件TGNET的压力、温度、甲烷摩尔分数、氢气摩尔分数仿真结果，如图9所示。

由图9可知：在交付节点3#位置处，压力和温度的平均绝对偏差分别为0.086 MPa和0.45 K，甲烷和氢气摩尔分数的平均绝对偏差分别为0.098%和0.098%。在节点1#位置处，压力和温度的平均绝对偏差分别为0.072 MPa和0.44 K，甲烷和氢气摩尔分数的平均绝对偏差分别为0.10%和0.10%。这表明本文耦合模型的结果与TGNET仿真的结果非常接近，所建立的耦合模型可以为解耦模型的后续分析奠定基础。

为了得到与离散网格无关的仿真结果，本文解耦模型中分别选择30和60 s为时间步长及1.0、2.5和5.0 km为空间步长进行仿真。在上述时空间步长下，TGNET模型和本文解耦模型在交付节点3#和节点1#处的压力、温度、氢气摩尔分数和甲烷摩尔分数如图10所示。图10中，未写模型名称仅注明时空步长参数的为本文解耦模型结果。

由图10可以看出，交付节点3#和节点1#处本文解耦模型和TGNET计算的气体压力、温度、甲烷和氢气的摩尔分数的模拟结果非常接近。

此外，不同时空步长组合下本文解耦模型计算得出与TGNET模拟得出的甲烷摩尔分数的平均相对偏差及本文解耦模型计算所需时间，如图11所示。

由图11可见：当时间步长恒定时，甲烷摩尔分数的平均相对偏差随空间步长的减小而减小；当空间步长恒定时，甲烷摩尔分数的平均相对偏差随时间步长的减小而减小。

综合考虑到本文解耦模型的计算时间小及本文解耦模型与TGNET模拟相比的甲烷摩尔分数的平均相对偏差不大，选择2.5 km和30 s作为空间步长和时间步长，模拟以下瞬态场景。在2.5 km时空步长下模拟1 h，本文的耦合模型和解耦模型的计算时间分别为4 383 s和1 050 s，本文解耦模型在所选时空步长下的计算时间比耦合模型的计算时间缩短了76.04%。

在天然气管网实际运行中，供气节点和交付节点的边界条件组合可能多种多样。本文所提解耦模型不仅对上述边界条件组合有效，而且适用于所有类型的边界条件组合。

3.2.1 边界条件组合A

供气节点1#的压力控制边界条件改为质量流量控制，该质量流量在整个时间范围内为100 kg/s，其余边界条件与第3.1节中相同；供气节点2#的边界条件由质量流量控制改为压力控制，整个时间范围的控制压力为9.72 MPa，其余边界条件也与第3.1节中相同。在上述边界条件组合A下，分析交付节点3#和节点1#处的本文解耦模型和TGNET模拟的压力、温度、甲烷摩尔分数、氢气摩尔分数结果，如图12所示。

由图12可知，交付节点3#位置处的压力、温度、甲烷摩尔分数和氢气摩尔分数的平均绝对偏差分别为0.078 MPa、0.39 K、0.21%、0.21%，节点1#位置处的压力、温度、甲烷摩尔分数和氢气摩尔分数的平均绝对偏差分别为0.067 MPa、0.43 K、0.14%和0.14%。

3.2.2 边界条件组合B

供气节点2#的质量流量控制边界条件改为压力控制，在整个时间范围内控制压力为9.72 MPa，其余边界条件与第3.1节中相同。在上述边界条件组合B下，分析交付节点3#和节点1#处的本文解耦模型和TGNET的压力、温度、甲烷摩尔分数、氢气摩尔分数模拟结果，如图13所示。

由图13可以知道，交付节点3#位置处的压力、温度、甲烷摩尔分数和氢气摩尔分数的平均绝对偏差分别为0.078 MPa、0.41 K、0.21%、0.21%，节点#1位置处的平均绝对偏差分别为0.068 MPa、0.44 K、0.13%和0.13%。

3.2.3 边界条件组合C

在天然气管网调峰的实际运行中，供应节点和交付节点通常会使用质量流量控制。该情况中，供气节点1#的控制方式由压力控制改为质量流量控制，其余边界条件与第3.1节中相同。在整个时间范围内，供气节点1#的控制质量流量为100 kg/s。在上述边界条件组合C下，分析交付节点3#和节点1#处的本文解耦模型的压力、温度、甲烷摩尔分数、氢气摩尔分数模拟结果，如图14所示。需要注意的是，TGNET暂时还不能设定边界组合模式C。

4 结 论

本文分别建立了将流动参数和组分追踪耦合、解耦的仿真模型。通过对复杂管网进行仿真和分析，可以得到以下结论：

1）所提耦合模型和解耦模型在流动参数计算和气体组分追踪方面均具有较好的精度，达到了商业软件的水平。

2）所提耦合模型和解耦模型不仅适用于线性天然气管网，而且在多支路天然气管网仿真中也具有较高的精度。

3）在特定时空步长下，所提解耦模型计算时间比耦合模型计算时间缩短了76.04%，具有较高的仿真效率。所提解耦模型能够提供不同边界条件组合下的状态变量和气体组分的精确时空演化规律，表明该模型具有较强的边界条件适应性。

参考文献

[1] Hou Zhengmeng, Luo Jiashun, Cao Cheng,et al.Development and contribution of natural gas industry under the goal of carbon neutrality in China[J].Advanced Engineering Sciences,2023,55(1):243-252.

[2] 侯正猛,罗佳顺,曹成,等.中国碳中和目标下的天然气产业发展与贡献[J].工程科学与技术,2023,55(1):243-252.

[3] Fang Yanli, Hou Zhengmeng, Yue Ye,et al.A new concept of multifunctional salt cavern hydrogen storage applied to the integration of hydrogen energy industry[J].Advanced Engineering Sciences,2022,54(1):128-135.

[4] 方琰藜,侯正猛,岳也,等.一种应用于氢能产业一体化的新型多功能盐穴储氢库[J].工程科学与技术,2022,54(1):128-135.

[5] Li Jingfa, Su Yue, Zhang Heng,et al.Research progresses on pipeline transportation of hydrogen-blended natural gas[J].Natural Gas Industry,2021,41(4):137-152.

[6] 李敬法,苏越,张衡,等.掺氢天然气管道输送研究进展[J].天然气工业,2021,41(4):137-152.

[7] Xie Ping, Wu Yi, Li Changjun,et al.Research progress on pipeline transportation technology of hydrogenmixed natural gas[J].Oil & Gas Storage and Transportation,2021,40(4):361-370.

[8] 谢萍,伍奕,李长俊,等.混氢天然气管道输送技术研究进展[J].油气储运,2021,40(4):361-370.

[9] Zhang Duihong, Yang Yi.Development of a sectionalizing method for simulation of large-scale complicated natural gas pipeline networks[J].Journal of Pipeline Science and Engineering,2024,4(4):100209. doi:10.1016/j.jpse.2024.100209

[10] Zuo Lili, Jin Antai, Chen Qian,et al.Simulation on venting process and valve opening control method for gas trunk pipelines[J].Petroleum Science,2022,19(6):3016-3028. doi:10.1016/j.petsci.2022.09.027

[11] Zhou Dengji, Jia Xingyun, Ma Shixi,et al.Dynamic simulation of natural gas pipeline network based on interpretable machine learning model[J].Energy,2022,253:124068. doi:10.1016/j.energy.2022.124068

[12] Yu Weichao, Huang Weihe, Wen Kai,et al.Subset simulation-based reliability analysis of the corroding natural gas pipeline[J].Reliability Engineering & System Safety,2021,213:107661. doi:10.1016/j.ress.2021.107661

[13] Wang Peng, Yu Bo, Han Dongxu,et al.Fast method for the hydraulic simulation of natural gas pipeline networks based on the divide-and-conquer approach[J].Journal of Natural Gas Science and Engineering,2018,50:55-63. doi:10.1016/j.jngse.2017.11.017

[14] Madoliat R, Khanmirza E, Moetamedzadeh H R.Transient simulation of gas pipeline networks using intelligent methods[J].Journal of Natural Gas Science and Engineering,2016,29:517-529. doi:10.1016/j.jngse.2016.01.018

[15] Chaczykowski M, Zarodkiewicz P.Simulation of natural gas quality distribution for pipeline systems[J].Energy,2017,134:681-698. doi:10.1016/j.energy.2017.06.020

[16] Hong Bingyuan, Li Xiaoping, Song Shangfei,et al.Optimal planning and modular infrastructure dynamic allocation for shale gas production[J].Applied Energy,2020,261:114439. doi:10.1016/j.apenergy.2019.114439

[17] Zhang Jiaqi, Li Zheng, Zheng Xiaoying,et al.Long-term planning and coupling optimization of multi-regional natural gas and hydrogen supply systems:A case study of China[J].Computers & Chemical Engineering,2024,183:108593. doi:10.1016/j.compchemeng.2024.108593

[18] Pinto G M, da Costa R B R, de Souza T A Z,et al.Experimental investigation of performance and emissions of a CI engine operating with HVO and farnesane in dual-fuel mode with natural gas and biogas[J].Energy,2023,277:127648. doi:10.1016/j.energy.2023.127648

[19] Chen Qian, Guan Aocheng, Chen Feng,et al.A transient gas pipeline network simulation model for decoupling the hydraulic-thermal process and the component tracking process[J].Energy,2024,301:131613. doi:10.1016/j.energy.2024.131613

[20] Gong Jing, Yin Xiong, Li Weijia,et al.Exploration on the function and operation mode of natural gas pipeline networks in energy Internet[J].Oil & Gas Storage and Transportation,2022,41(6):702-711.

[21] 宫敬,殷雄,李维嘉,等.能源互联网中的天然气管网作用及其运行模式探讨[J].油气储运,2022,41(6):702-711.

[22] Bleschke T, Chaczykowski M.Composition tracking of natural gas-hydrogen mixtures in pipeline flow using high-resolution schemes[J].International Journal of Hydrogen Energy,2024,79:756-770. doi:10.1016/j.ijhydene.2024.06.402

[23] Yang Zhaoming, Liu Zhe, Zhou Jing,et al.A graph neural network (GNN) method for assigning gas calorific values to natural gas pipeline networks[J].Energy,2023,278:127875. doi:10.1016/j.energy.2023.127875

[24] Zhou Weiguo, Liu Xiaojing, Wang Hai.Development of gas pipeline networks simulation technologies[J].World Sci-Tech R & D,2013,35(1):99-100.

[25] 周伟国,刘晓婧,王海.燃气管网仿真技术的发展状况[J].世界科技研究与发展,2013,35(1):99-100.

[26] Dorao C A, Fernandino M.Simulation of transients in natural gas pipelines[J].Journal of Natural Gas Science and Engineering,2011,3(1):349-355. doi:10.1016/j.jngse.2011.01.004

[27] Chaczykowski M.Transient flow in natural gas pipeline—The effect of pipeline thermal model[J].Applied Mathematical Modelling,2010,34(4):1051-1067. doi:10.1016/j.apm.2009.07.017

[28] Abbaspour M, Chapman K S, Glasgow L A.Transient modeling of non-isothermal,dispersed two-phase flow in natural gas pipelines[J].Applied Mathematical Modelling,2010,34(2):495-507. doi:10.1016/j.apm.2009.06.023

[29] Osiadacz A J, Chaczykowski M.Comparison of isothermal and non-isothermal pipeline gas flow models[J].Chemical Engineering Journal,2001,81(1/2/3):41-51. doi:10.1016/s1385-8947(00)00194-7

[30] Wang Peng, Yu Bo, Han Dongxu,et al.Adaptive implicit finite difference method for natural gas pipeline transient flow[J].Oil & Gas Sciences and Technology,2018,73:21. doi:10.2516/ogst/2018013

[31] Ding Yanpeng, Li Yuxing, Fu Juntao,et al.Effect of selection of basic parameters on simulation of fast transients in natural gas pipelines[J].Journal of Petrochemical Universities,2011,24(3):82-86.

[32] 丁延鹏,李玉星,付俊涛,等.基础变量选取对输气管道快瞬变流数值模拟的影响[J].石油化工高等学校学报,2011,24(3):82-86.

[33] Chaczykowski M, Sund F, Zarodkiewicz P,et al.Gas composition tracking in transient pipeline flow[J].Journal of Natural Gas Science and Engineering,2018,55:321-330. doi:10.1016/j.jngse.2018.03.014

[34] Zhang Zihang, Saedi I, Mhanna S,et al.Modelling of gas network transient flows with multiple hydrogen injections and gas composition tracking[J].International Journal of Hydrogen Energy,2022,47(4):2220-2233. doi:10.1016/j.ijhydene.2021.10.165

[35] Bermúdez A, Shabani M.Numerical simulation of gas composition tracking in a gas transportation network[J].Energy,2022,247:123459. doi:10.1016/j.energy.2022.123459

[36] Fan Di, Gong Jing, Zhang Shengnan,et al.A transient composition tracking method for natural gas pipe networks[J].Energy,2021,215:119131. doi:10.1016/j.energy.2020.119131

[37] 李玉星,姚光镇.输气管道设计与管理[M].2版.东营:中国石油大学出版社,2009.

[38] Yang Guangfeng, Yang Fan, Yang Tinghao,et al.Predication of five state equations on saturated gas-liquid phase density[J].Oil-Gasfield Surface Engineering,2018,37(8):9-12.

[39] 杨广峰,杨帆,杨霆浩,等.五种状态方程对饱和气液相密度的预测[J].油气田地面工程,2018,37(8):9-12.

[40] Yuan Weimin.Use of equations of state in predicting thermophysical properties of LNG[J].Petroleum Engineering Construction,2018,44(1):23-26.

[41] 苑伟民.应用状态方程预测液化天然气的热物理性质[J].石油工程建设,2018,44(1):23-26.

基金资助

国家自然科学基金面上项目(52174064)

湖北省教育厅科研计划项目青年人才项目(Q20241308)