ANALYTICAL METHOD FOR DYNAMIC ANALYSIS OF NATURAL GAS NETWORK

Abstract

Disclosed is an analytical method for dynamic analysis of a natural gas network in the field of energy system modeling and operational analysis, which includes establishing adynamic model of natural gas transmission according to the conservation equations, and reconstructing the dynamic model into the equations in a heat conduction equation form. The present disclosure directly constructs an analytical method for dynamic analysis of a natural gas network, avoiding approximation errors, numerical dispersion, and dissipation compared with the traditional numerical methods. The discretization process is avoided during the solution, greatly improving the computational efficiency and solution accuracy of dynamic analysis of the natural gas network.

Description

TECHNICAL FIELD

The present disclosure relates to the field of energy system modeling and operational analysis, and in particular to an analytical method for dynamic analysis of a natural gas network.

BACKGROUND

The growing energy consumption and environmental pressures are driving the transformation of low-carbon green energy network technologies. Cities, as the main subjects of energy consumption and transformation, are transitioning towards complex energy networks with multiple energy flows and dynamics. As an important component of the urban energy network, a natural gas system can significantly improve the overall utilization efficiency of energy system and enhance the renewable consumption through the integration with the power system and heating system, as well as through optimized energy management. The time scales for the operation and management of different energy networks vary significantly, making it necessary to obtain real-time, reliable, and consistent network information based on accurate simulation models and technologies. However, since multi-energy networks are managed and operated by different companies, the exchange of information is highly limited, requiring special attention to information protection during co-simulation and operational optimization.

The simulation of the natural gas system essentially involves defining a set of state variables to describe the key characteristics of the system, followed by mechanism analysis to acquire the evolution of all state variables under a given boundary. Since the natural gas system model is composed of a set of partial differential-algebraic equations, one of the mainstream methods is to discretize the pipeline model using space-time segmentation and recursively calculate the state distribution within the system based on boundary and initial conditions. However, to ensure calculation efficiency, a large number of segments is typically required for each pipeline, which results in a low computational efficiency for the recursive process. Additionally, it is difficult to directly and quantitatively characterize the response of the state variables in the natural gas system to external boundaries.

In light of the above shortcomings, it is necessary to develop an analytical model for the transmission dynamics of the natural gas system that balances both accuracy and complexity from a modeling perspective. This model would allow for the rapid and precise description of the system's dynamic processes, making it more conducive to multi-agent simulation and optimization operation, and better suited to the practical application conditions in engineering.

SUMMARY

In view of the deficiencies of the prior art, the present disclosure proposes an analytical method for dynamic analysis of a natural gas network. The disclosure utilizes a constant variation method to derive an equivalent form of a dynamic transmission model of a natural gas system. Furthermore, based on the superposition, time-invariant, and causal properties of natural gas transmission, practical analytical forms of the dynamic model under different initial conditions and boundary conditions, as well as the equivalent network model, are constructed, providing a modeling basis for the operational analysis of an integrated energy system.

An object of the present disclosure may be achieved by the following technical solutions:

An analytical method for dynamic analysis of a natural gas network is provided, including:

• • establishing a dynamic model of natural gas transmission according to the conservation equations, and reconstructing the dynamic model into the equations in a heat conduction equation form;
  • transforming the heat conduction equation into a partial differential equation set with homogeneous boundary conditions using a method of variation of parameters, and constructing a general analytical expression for the partial differential equation set under constant pressure boundaries;
  • constructing, based on superposition principle, practical analytical expressions for the dynamic model of natural gas transmission;
  • constructing, according to system topologies and time-invariant, and causal properties of a natural gas system, a natural gas network equivalent model based on the practical analytical expressions; and
  • determining pressures at load nodes and mass flows at source nodes based on the natural gas network equivalent model.

In some embodiments, the establishing a dynamic model of natural gas transmission according to the conservation equations includes the following steps:

• • establishing a dynamic model of natural gas transmission in a pipeline, including a mass conservation equation and a momentum conservation equation as below:

$\begin{array}{cc}\frac{\partial p_{x,t}}{\partial t}+\frac{c^2}{S}\frac{\partial q_{x,t}}{\partial x}=0,\frac{\partial p_{x,t}}{\partial x}+\frac{\lambda {\mathrm{wq}}_{x,t}}{2\mathrm{SD}}=0,& \left(1\right)\end{array}$

• • where x and t represent space and time, respectively, and c represents a speed of sound, taken as 340 m/s; p_x,t represents a pressure of natural gas, Pa; q_x,t represents mass flow of natural gas, kg/s; w represents an average velocity of natural gas, m/s; D represents an inner diameter of the pipeline, m; S represents a cross-sectional area of the pipeline, m²; and λ represents a resistance coefficient of the pipeline; and
  • setting initial conditions and boundary conditions in the dynamic analysis of the natural gas system as follows:

$\begin{array}{cc}{p}_{x,0}={\phi }_{p,x},q_{x,0}={\phi }_{q,x},p_{0,t},q_{L,t},& \left(2\right)\end{array}$

• • where p_x,0 and q_x,0 represent distributions of pressure and mass flow in the pipeline at moment 0 respectively, namely, the initial conditions; φ_p,x and φ_q,x represent initial condition functions of pressure and mass flow, respectively; L represents a length of the pipeline, m; p_0,t represents a time-series of pressure at an inlet of the pipeline, namely, pressure boundary; and q_L,t represents a time-series of mass flow at an outlet of the pipeline, namely, flow boundary.

In some embodiments, the reconstructing the dynamic gas flow model into the equations in a heat conduction equation form includes the following steps:

deriving a partial derivative for x from a second equation in Eq. before substituting into a first equation in Eq. (1) to obtain a heat conduction form equation as follows:

$\begin{array}{cc}\frac{\partial q_{x,t}}{\partial x}=\frac{-2\mathrm{DS}}{\lambda w}\frac{{\partial }^2{p}_{x,t}}{\partial x^2}\Rightarrow \{\begin{array}{c}\frac{\partial p_{x,t}}{\partial t}=\alpha \frac{{\partial }^2{p}_{x,t}}{\partial x^2}\\ \alpha =\frac{2{\mathrm{Dc}}^2}{\lambda w}\end{array},& \left(3\right)\end{array}$

• • where α represents a coefficient term in the heat conduction form equation; accordingly, the partial derivative of pipeline outlet pressure for x can be obtained through the following equation and is defined as m_t.

$\begin{array}{cc}\frac{\partial p}{\partial x}{❘}_{x=L}=-\frac{\lambda w}{2\mathrm{DS}}{q}_{L,\iota }=m_{\iota },& \left(4\right)\end{array}$

Eq. (1) is converted to:

$\begin{array}{cc}\{\begin{array}{c}\frac{\partial p_{x,t}}{\partial t}=\alpha \frac{{\partial }^2{p}_{x,t}}{\partial x^2},\frac{\partial p_{x,t}}{\partial x}{❘}_{x=L}=m_t\\ p\left(0,t\right)=p_{0,t},p\left(x,0\right)={\phi }_{p,x}\end{array}.& \left(5\right)\end{array}$

In some embodiments, the transforming the heat conduction equation into a partial differential equation set with homogeneous boundary conditions using a method of variation of parameters includes the following steps:

• • taking p_0,t as constant with the value of ψ_p,1, introducing an intermediate variable u_x,t and reconstructing p_x,t as follows:

$\begin{array}{cc}{p}_{x,t}=u_{x,t}+x\left(\frac{\partial p}{\partial x}{❘}_{x=L}\right)+p_{0,t};& \left(6\right)\end{array}$

• • substituting x=0 and t=0 into Eq. (6) to obtain initial conditions and a set of boundary conditions for u_x,t, respectively:

$\begin{array}{cc}\begin{array}{c}\{\begin{array}{c}{u}_{x,0}=p_{x,0}-{x\left({\frac{\partial p}{\partial x}|}_{x=L}\right)|}_{t=0}-p_{0,0}\\ p_{0,t}=u_{0,t}+p_{0,t}\Rightarrow u_{0,t}=0\end{array}\end{array};& \left(7\right)\end{array}$

• • deriving a partial derivative for x on both sides of Eq. (6) at x=L to obtain a further set of boundary conditions for u_x,t, expressed as follows:

$\begin{array}{cc}{\frac{\partial p_{x,t}}{\partial x}=\frac{\partial u_{x,t}}{\partial x}+{\frac{\partial p_{x,t}}{\partial x}|}_{x=L}\Rightarrow \frac{\partial u_{x,t}}{\partial x}❘}_{x=L}=0;& \left(8\right)\end{array}$

• • deriving a second-order partial derivative for x and a first-order partial derivative for t on both sides of Eq. (6), respectively:

$\begin{array}{cc}{\frac{{\partial }^2{p}_{x,t}}{\partial x^2}=\frac{{\partial }^2{p}_{x,t}}{\partial x^2}+\frac{\partial p_{x,t}}{\partial t}=\frac{\partial u_{x,t}}{\partial t}+x\frac{\partial }{\partial t}(\frac{\partial p_{x,t}}{\partial x}❘}_{x=L});& \left(9\right)\end{array}$

and

• • substituting Eq. (9) into Eq. (5), with Eq. (5) re-expressing a partial differential equation set with homogeneous boundary conditions as follows:

$\begin{array}{cc}\{\begin{array}{c}\frac{\partial u_{x,t}}{\partial t}=\alpha \frac{{\partial }^2{u}_{x,t}}{\partial x^2}+f_{x,t}\\ {{u_{x,0}=p_{x,0}-x(\frac{\partial p}{\partial x}❘}_{x=L})❘}_{t=0}-p_{0,0}\frac{\partial u_{x,t}}{\partial x}{|}_{x=L}=0.\\ u_{0,t}=0,f_{x,t}=-x\frac{\partial }{\partial t}\left({\frac{\partial p}{\partial x}|}_{x=L}\right)=-x\frac{\partial m_t}{\partial t}\end{array}& \left(10\right)\end{array}$

In some embodiments, the constructing a general analytical expression for the partial differential equation set under constant pressure boundaries includes the following steps:

• • defining u_x,t* as a solution of a homogeneous problem of Eq. (10), and solving the homogeneous problem as follows:

$\begin{array}{cc}{u}_{x,t}^{*}=T_0\left(t\right){\sum }_{n=1}^{\infty }{k}_2\sin \frac{\left(2n-1\right)\pi x}{2L},& \left(11\right)\end{array}$

• • where k₂ represents a coefficient term; T₀(t) represents a function component about time in u_x,t*; n represents the number of Fourier components; analogy to the form of u_x,t*, a form of u_x,t in Eq. (10) is assumed as follows:

$\begin{array}{cc}{u}_{x,t}={\sum }_{n=1}^{\infty }{u}_{n,t}\sin \frac{\left(2n-1\right)\pi x}{2L},& \left(12\right)\end{array}$

• • where u_n,t represents a function component about time in u_x,t; f_x,t in Eq. (10) is expanded according to a base in Eq. (11):

$\begin{array}{cc}{f}_{x,t}={\sum }_{n=1}^{\infty }{f}_{n,t}\sin \frac{\left(2n-1\right)\pi x}{2L},& \left(13\right)\end{array}$ $\begin{array}{cc}{f}_{n,t}=\frac{2}{L}{\int }_0^L{f}_{x,t}\sin \frac{\left(2n-1\right)\pi x}{2L}\mathrm{dx},& \left(14\right)\end{array}$

• • where f_n,t represents a function component about time in f_x,t, obtained by inverse Fourier transformation shown in Eq. (14); Eq. (12) to (14) are substituted into a first equation in Eq. (10):

$\begin{array}{cc}\{\begin{array}{c}\frac{\partial u_{n,t}}{\partial t}+g_n{u}_{n,t}=f_{n,t}\\ g_n=\frac{{\left(2n-1\right)}^2{\pi }^2\alpha }{4L^2}\\ u_{n,0}=\frac{2}{L}{\int }_0^L{u}_{x,0}\sin \frac{\left(2n-1\right)\pi x}{2L}\mathrm{dx}\end{array},& \left(15\right)\end{array}$

• • where g_n represents a coefficient term in an n-th nonhomogeneous partial differential, u_n,0 represents a value of u_n,t when t=0; Eq. (15) is solved to obtain analytical expressions for u_n,t and u_x,t, respectively as follows:

$\begin{array}{cc}{u}_{n,t}=e^{-g_{n}t}\left\{u_{n,0}+{\int }_0^t{f}_{n,\tau }{e}^{g_n\tau }d\tau \right\},& \left(16\right)\end{array}$ $\begin{array}{cc}{u}_{x,t}={\sum }_{n=1}^{\infty }{e}^{-g_{n}t}\left\{u_{n,0}+{\int }_0^t{f}_{n,\tau }{e}^{g_n\tau }d\tau \right\}\sin \frac{\left(2n-1\right)\pi x}{2L};& \left(17\right)\end{array}$

• • substituting Eq. (17) into Eq. (6) to obtain an analytical expression for p_x,t at constant pressure:

$\begin{array}{cc}{p}_{x,t}={\sum }_{n=1}^{\infty }{e}^{-g_{n}t}\left\{u_{n,0}+{\int }_0^t{f}_{n,\tau }{e}^{g_n\tau }d\tau \right\}\sin \frac{\left(2n-1\right)\pi x}{2L}-\frac{\lambda \mathrm{wx}}{2\mathrm{Ds}}{q}_{L,t}+p_{0,t};& \left(18\right)\end{array}$

• • substituting Eq. (18) into a second equation in Eq. (1) to obtain an analytical expression for q_x,t at constant pressure:

$\begin{array}{cc}{q}_{x,t}=-\frac{2\mathrm{DS}}{\lambda w}\frac{\partial p_{x,t}}{\partial x}=q_{L,t}+\left\{{\sum }_{n=1}^{\infty }\frac{\left(1-2n\right)\pi \mathrm{DS}}{\lambda \mathrm{wL}{e}^{g_{n}t}}\left\{u_{n,0}+{\int }_0^t{f}_{n,\tau }{e}^{g_n\tau }d\tau \right\}\cos \frac{\left(2n-1\right)\pi x}{2L}\right\},& \left(19\right)\end{array}$

• • where p_0,t represents a constant pressure boundary; Eq. (18) and Eq. (19) are general analytical expressions for the partial differential equation set of the natural gas dynamic model under constant pressure boundaries.

In some embodiments, the constructing, based on superposition principle, practical analytical expressions for the dynamic model of natural gas transmission includes the following steps:

• • reconstructing the pressure boundary, expressed as follows:

$\begin{array}{cc}\{\begin{array}{c}{\psi }_p=\left[{\psi }_{p,1},{\psi }_{p,2},\dots ,{\psi }_{p,\mathrm{Nt}}\right]\\ {\psi }_{p,i}=\{\begin{array}{cc}{p}_{0,i}& i=1\\ p_{0,i}-p_{0,i-1}& i>1\end{array},t_i=i\Delta t\end{array},& \left(20\right)\end{array}$

• • where Nt represents a length of the boundary condition, and Δt is a time step; ψ_p represents a new boundary condition corresponding to a pressure increment, and an element ψ_p,i represents an increment of a pressure at moment t_i for a pressure at moment t_i−1;
  • decomposing, according to superposition principle of natural gas transmission, the pressures p_x,t as follows:

$\begin{array}{cc}{p}_{x,t}={\sum }_{k=1}^i{\chi }_{k,x,t-t_{k-1}}1\le i\le \mathrm{Nt},t_{k-1}\le t\le t_i,& \left(21\right)\end{array}$

• • where each component χ_k,x,t satisfies:

$\begin{array}{cc}\{\begin{array}{c}\frac{\partial {\chi }_{k,x,t-t_{k-1}}}{\partial t}=\alpha \frac{{\partial }^2{\chi }_{k,x,t-t_{k-1}}}{\partial x^2},\text{}t>t_{k-1},0<x<L\\ \frac{\partial {\chi }_{k,x,t-t_{k-1}}}{\partial t}{|}_{x=L}=\{\begin{array}{cc}{m}_t& k=1\\ 0& k>1\end{array},t>t_{k-1}\\ {\chi }_{k,x,t-t_{k-1}}{|}_{x=0}=\begin{array}{cc}{\psi }_{p,k}\text{},& t>t_{k-1}\end{array}\\ {\chi }_{k,x,t-t_{k-1}}{|}_{t=t_{k-1}}=\{\begin{array}{cc}{\phi }_{p,x}& k=1,\\ 0& k>1\end{array},\text{}0\le x\le L\end{array},& \left(22\right)\end{array}$

• • when k=1, a form of χ_k,x,t is consistent with Eq. (18), namely:

$\begin{array}{cc}{\chi }_{k,x,t-t_{k-1}}={\sum }_{n=1}^{\infty }{e}^{-g_{n}t}\left\{u_{n,0}+{\int }_0^t{f}_{n,\tau }{e}^{g_n\tau }d\tau \right\}\sin \frac{\left(2n-1\right)\pi x}{2L}\begin{array}{cc}+{\mathrm{xm}}_t+{\psi }_{p,1}& k=1,t>t_{k-1}\end{array},& \left(23\right)\end{array}$

• • when k>1, a form of χ_k,x,t is as follows:

$\begin{array}{cc}{\chi }_{k,x,t-t_{k-1}}={\psi }_{p,k}k>1,t>t_{k-1};& \left(24\right)\end{array}$

• • calculating the initial conditions according to steady-state energy flow in practice, expressed as follows:

$\begin{array}{cc}{\phi }_{p,x}=K_{1}x+K_2,{\phi }_{q,x}=K_3,m_0=\frac{-\lambda w{\phi }_{q,L}}{2\mathrm{DS}},& \left(25\right)\end{array}$

• • where K₁, K₂, and K₃ represent constant coefficients in the initial conditions; m₀ represents a value of m_t at t=0; and φ_q,L represents a value of the initial conditions of mass flow at x=L, namely, the mass flow at an outlet of the pipeline at t=0;
  • expressing m_t as a segmented step function as follows:

$\begin{array}{cc}{m}_t={\sum }_{i=1}^{\mathrm{Nt}}{m}_i\left({\delta }_{t-t_i}-{\delta }_{t-t_{i-1}}\right)t>0;& \left(26\right)\end{array}$

• • substituting f_x,t in Eq. (10) and f_n,t in Eq. (14) into Eq. (16) as follows:

$\begin{array}{cc}\begin{array}{c}={\int }_0^t\left(\frac{-2}{L}{\int }_0^L\times \frac{\partial m_{\tau }}{\partial \tau }\sin \frac{\left(2n-1\right)\pi x}{2L}\mathrm{dx}\right)e^{g_n^{\tau }}d\tau \\ {\int }_0^t{f}_{n,\pi }{e}^{g_n^{\tau }}d\tau =\frac{-2}{L}{\int }_0^{L}x\sin \frac{\left(2n-1\right)\pi x}{2L}\mathrm{dx}{\int }_0^t\frac{\partial m_{\tau }}{\partial \tau }{e}^{g_n^{\tau }}d\tau \\ =\frac{8{L\left(-1\right)}^n}{{\left(2n-1\right)}^2{\pi }^2}{\int }_0^t\frac{\partial m_{\tau }}{\partial \tau }{e}^{g_n^{\tau }}d\tau \end{array};& \left(27\right)\end{array}$

• • modifying an integral formula in Eq. (27) as follows:

$\begin{array}{cc}\frac{\partial m_{\tau }}{\partial \tau }{e}^{g_n\tau }=\frac{\partial \left(m_{\tau }{e}^{g_n\tau }\right)}{\partial \tau }-g_n{m}_{\tau }{e}^{g_n\tau },& \left(28\right)\end{array}$ $\begin{array}{cc}{g}_n{\int }_0^{t_k}{m}_{\tau }{e}^{g_n\tau }d\tau ={\sum }_{i=1}^k{m}_i\left(e^{g_n{t}_i}-e^{g_n{t}_{i-1}}\right);& \left(29\right)\end{array}$

• • substituting Eq. (28) and Eq. (29) into Eq. (27) as follows:

$\begin{array}{cc}\frac{1}{e^{g_n{t}_k}}{\int }_0^{t_k}{f}_{n,\tau }{e}^{g_n\tau }d\tau =-E_{1n}\left(m_k-\frac{m_0}{e^{g_n{t}_k}}-{\sum }_{i=1}^k{m}_i\frac{e^{g_n\Delta t}-1}{e^{g_n\left(t_k-t_{i-1}\right)}}\right),& \left(30\right)\end{array}$ $\begin{array}{cc}{E}_{1n}=\frac{8{L\left(-1\right)}^{n-1}}{{\left(2n-1\right)}^2{\pi }^2};& \left(31\right)\end{array}$

and

• • substituting Eq. (30) into Eq. (18) and Eq. (19), to obtain the practical analytical expressions for the dynamic model of natural gas transmission as follows:

$\begin{array}{cc}{p}_{x,k}={\sum }_{n=1}^{\infty }\left\{\frac{u_{n,0}}{e^{g_n{t}_k}}-E_{1n}\left(m_k-\frac{m_0}{e^{g_n{t}_k}}-{\sum }_{i=1}^k{m}_i\frac{e^{g_n\Delta t}-1}{e^{g_n\left(t_k-t_{i-1}\right)}}\right)\right\}\sin \frac{\left(2n-1\right)\pi x}{2L}-\frac{\lambda \mathrm{wx}}{2\mathrm{DS}}{q}_{L,k}+p_{0,k},& \left(32\right)\end{array}$ $\begin{array}{cc}{q}_{x,k}=q_{L,k}+\left\{{\sum }_{n=1}^{\infty }\frac{\left(1-2n\right)\pi \mathrm{DS}}{\lambda \mathrm{wL}}\left\{\frac{u_{n,0}}{e^{g_n{t}_k}}-E_{1n}\left(m_k-\frac{m_0}{e^{g_n{t}_k}}-{\sum }_{i=1}^k{m}_i\frac{e^{g_n\Delta t}-1}{e^{g_n\left(t_k-t_{i-1}\right)}}\right)\right\}\cos \frac{\left(2n-1\right)\pi x}{2L}\right\}.& \left(33\right)\end{array}$

In some embodiments, the constructing, according to system topologies and time-invariant, and causal properties of a natural gas system, a natural gas network equivalent model based on the practical analytical expressions includes the following steps:

• • substituting x=L and x=0 into Eq. (32) and Eq. (33), respectively, to obtain pressure at the outlet and flow at the inlet of the pipeline at arbitrary t_j as follows:

$\begin{array}{cc}{p}_{L,j}={\sum }_{n=1}^{\infty }{e}^{-g_n{t}_j}{u_{n,0}\left(-1\right)}^{n-1}-\frac{\lambda \mathrm{wL}}{2\mathrm{DS}}{q}_{L,j}+p_{0,j}+{\sum }_{n=1}^{\infty }\frac{\lambda {\mathrm{wE}}_{1n}}{2\mathrm{DS}}{\left(-1\right)}^{n-1}\left(q_{L,j}-\frac{q_{L,0}}{e^{g_n{t}_j}}-{\sum }_{i=1}^j{q}_{L,i}\frac{e^{g_n\Delta t}-1}{e^{g_n\left(j-i+1\right)\Delta t}}\right),& \left(34\right)\end{array}$ $\begin{array}{cc}{q}_{0,j}=q_{L,j}+\frac{\mathrm{DS}\pi }{\lambda \mathrm{wL}}{\sum }_{n-1}^{\infty }\left(1-2n\right)u_{n,0}{e}^{-g_n{t}_j}+\frac{\pi }{2L}{\sum }_{n-1}^{\infty }\left(1-2n\right)E_{1n}\left(q_{L,j}-\frac{q_{L,0}}{e^{g_n{t}_j}}-{\sum }_{i=1}^j{q}_{L,i}\frac{e^{g_n\Delta t}-1}{e^{g_n\left(j-i+1\right)\Delta t}}\right);& \left(35\right)\end{array}$

• • reorganizing Eq. (34) and Eq. (35) as follows:

$\begin{array}{cc}\left[\begin{array}{c}{p}_L\\ q_0\end{array}\right]=\left[\begin{array}{cc}{H}_1& H_2\\ H_3& H_4\end{array}\right]\left[\begin{array}{c}{p}_0\\ q_L\end{array}\right]+\left[\begin{array}{c}{\gamma }_p\\ {\gamma }_q\end{array}\right],& \left(36\right)\end{array}$

• • where p_L and p₀ represent a pressure vector at the outlet and inlet of the pipeline, respectively; q_L and q₀ represent a flow vector at the outlet and inlet of the pipeline, respectively; H₁ to H₄ represent constant coefficient matrices with the dimension of N_t×N_t; and γ_p and γ_q represent a constant coefficient vector of pressure and mass flow with the dimension of N_t×1, respectively, determined by the initial conditions; the elements may be represented as follows:

$\begin{array}{cc}\{\begin{array}{c}{p}_0=\left[p_{0,1},\dots ,p_{0,\mathrm{Nt}}\right],p_L=\left[p_{L,1},\dots ,p_{L,\mathrm{Nt}}\right]\\ q_L=\left[q_{L,1},\dots ,q_{L,\mathrm{Nt}}\right],q_0=\left[q_{0,1},\dots ,q_{0,\mathrm{Nt}}\right]\end{array},& \left(37\right)\end{array}$ $\begin{array}{cc}{\gamma }_{p,j}={\sum }_{n1}^{\infty }\left(u_{n,0}+E_{1n}{m}_0\right){\left(-1\right)}^{n-1}{e}^{-g_n{t}_j},& \left(38\right)\end{array}$ $\begin{array}{cc}{\gamma }_{q,j}=\frac{\mathrm{DS}\pi }{\lambda \mathrm{wL}}{\sum }_{n-1}^{\infty }\left(1-2n\right)\left(u_{n,0}+E_{1n}{m}_0\right)e^{-g_n{t}_j},& \left(39\right)\end{array}$ $\begin{array}{cc}{H}_1=1,H_3=0,\underset{i>j}{H_{2.\mathrm{ji}}}=\underset{i>j}{H_{4.\mathrm{ji}}}=0,& \left(40\right)\end{array}$ $\begin{array}{cc}\underset{i=j}{H_{2,\mathrm{ji}}}=\frac{\lambda w}{2\mathrm{DS}}\left({\sum }_{n=1}^{\infty }{E_{1n}\left(-1\right)}^{n-1}{e}^{-g_n\Delta t}-L\right),& \left(41\right)\end{array}$ $\begin{array}{cc}\underset{\left(i<j\right)}{H_{2,\mathrm{ji}}}=\frac{\lambda w}{2\mathrm{DS}}\left({\sum }_{n=1}^{\infty }{E_{1n}\left(-1\right)}^{n-1}\frac{1-e^{g_n\Delta t}}{e^{g_n\left(j-i+1\right)\Delta t}}\right),& \left(42\right)\end{array}$ $\begin{array}{cc}\underset{\left(i=j\right)}{H_{4,\mathrm{ij}}}=\frac{\pi }{2L}{\sum }_{n=1}^{\infty }{E}_{1n}\left(1-2n\right)e^{-g_n\Delta t}+1,& \left(43\right)\end{array}$ $\begin{array}{cc}\underset{\left(i<j\right)}{H_{4,\mathrm{ji}}}=\frac{\pi }{2L}{\sum }_{n=1}^{\infty }{E}_{1n}\left(1-2n\right)\frac{1-e^{g_n\Delta t}}{e^{g_n\left(j-i+1\right)\Delta t}};& \left(44\right)\end{array}$

• • defining the number of nodes and pipelines in the natural gas system as N_g and N_b, respectively, and defining an identity matrix and a zero matrix as 1 and 0, respectively; introducing incidence matrices A_in, A_out, A_nb1, and A_nb2 to describe topological properties in the natural gas system,
  • where A_in includes N_g×N_b sub-matrices for associating flow at end of a branch with flow at a node, and in response to a sub-matrix of an i-th row and j-th column being 1, the flow flows from the end of the branch j into the node i, otherwise, the sub-matrix is 0; A_out includes N_g×N_b sub-matrices for associating flow at head end of a branch with flow at a node, and in response to a sub-matrix of an i-th row and j-th column being 1, the flow flows from the node i into the head end of the branch j, otherwise, the sub-matrix is 0; A_nb1 includes N_b×N_g sub-matrices for associating pressure at end of a branch with pressure at a node, and in response to a sub-matrix of an i-th row and j-th column being 1, pressure at the node j is equal to pressure at the end of the branch i, otherwise, the sub-matrix is 0; A_nb2 includes N_b×N_g sub-matrices for associating pressure at a head end of a branch with pressure at a node, and in response to a sub-matrix of an i-th row and j-th column being 1, pressure at the node j is equal to pressure at the head end of the branch i, otherwise, the sub-matrix is 0; particularly, in response to the node j being a compressor node with a pressure ratio of K_cp and the node j being connected to the head end of the branch i, the sub-matrix of the i-th row and jth column is K_cp×1;
  • constructing, based on the incidence matrices, equations of mass conservation at nodes and pressure continuity between nodes and pipelines, expressed as follows:

$\begin{array}{cc}{A}_{\mathrm{in}}{q}_L-A_{\mathrm{out}}{q}_0=\mathrm{qn};& \left(45\right)\end{array}$ $\begin{array}{cc}{p}_L=A_{\mathrm{nb}1}\mathrm{pn},p_0=A_{\mathrm{nb}2}\mathrm{pn};& \left(46\right)\end{array}$

• • where pn and qn represent a pressure and a vector at the node, respectively; p_L, p₀, q_L, and q₀ are all vectors of Nt×1; the sub-matrices 1 and 0 in coefficient matrices A_in, A_out, A_nb1, and A_out2 are all Nt×Nt matrices;
  • substituting Eq. (46) into a first row of Eq. (36) as follows:

$\begin{array}{cc}{q}_L=H_2^{-1}\left(\left(A_{\mathrm{nb}1}-A_{{\mathrm{nb}}_2}\right)\mathrm{pn}-{\gamma }_p\right)=H_{a1}\mathrm{pn}+H_{a2},& \left(47\right)\end{array}$

• • where H_a1 and H_a2 represent a constant coefficient matrix and a vector describing mass flow at an outlet of a pipeline, respectively; Eq. (46) and Eq. (47) are substituted into a second row of Eq. (36) as follows:

$\begin{array}{cc}{q}_0=H_4\left(H_{a1}\mathrm{pn}+H_{a2}\right)+{\gamma }_q=H_{a_3}\mathrm{pn}+H_{a_4},& \left(48\right)\end{array}$

• • where H_a3 and H_a4 represent a constant coefficient matrix and a vector describing mass flow at an inlet of the pipeline, respectively; Eq. (47) and Eq. (48) are substituted into Eq. (45) to associate pressure and mass flow at the node as follows:

$\begin{array}{cc}\mathrm{qn}\begin{array}{c}=\left(A_{\mathrm{in}}{H}_{a1}-A_{\mathrm{out}}{H}_{a3}\right)\mathrm{pn}+\left(A_{\mathrm{in}}{H}_{a2}-A_{\mathrm{out}}{H}_{a4}\right)\\ =H_{a5}\mathrm{pn}+H_{a6}\end{array},& \left(49\right)\end{array}$

• • where H_a5 and H_a6 are a constant coefficient matrix and a vector describing a relationship between pressure and mass flow at a node, respectively; subscripts sr and ns are defined as a variable identification of a gas source node and a gas load node, respectively, and Eq. (49) is expanded to obtain the natural gas network equivalent model based on the practical analytical expressions:

$\begin{array}{cc}\left[\begin{array}{c}{\mathrm{qn}}_{\mathrm{sr}}\\ {\mathrm{qn}}_{\mathrm{ns}}\end{array}\right]=\left[\begin{array}{cc}{H}_{a51}& H_{a52}\\ H_{a53}& H_{a54}\end{array}\right]\left[\begin{array}{c}{\mathrm{pn}}_{\mathrm{sr}}\\ {\mathrm{pn}}_{\mathrm{ns}}\end{array}\right]+\left[\begin{array}{c}{H}_{a61}\\ H_{a62}\end{array}\right],& \left(50\right)\end{array}$

• • where H_a51, H_a52, H_a53, and H_a54 are sub-matrices in H_a5, and H_a61 and H_a62 are sub-vectors in H_a6; pn_sr and pn_ns are a pressure vector at a source node and a load node, respectively; qn_sr and qn_ns are a mass flow vector at the source node and the load node, respectively; the mass flow at the source node and the pressure at the load node in the natural gas system are fully expressed as functions of the boundary conditions using matrix transformation, namely, a simplified method of natural gas network equivalent model is as follows:

$\begin{array}{cc}{\mathrm{pn}}_{\mathrm{ns}}=H_{a54}^{-1}\left({\mathrm{qn}}_{\mathrm{ns}}-H_{a53}{\mathrm{pn}}_{\mathrm{sr}}-H_{a62}\right),& \left(51\right)\end{array}$ $\begin{array}{cc}{\mathrm{qn}}_{\mathrm{sr}}=H_{a51}{\mathrm{pn}}_{\mathrm{sr}}+H_{a52}{\mathrm{pn}}_{\mathrm{ns}}+H_{a61}.& \left(52\right)\end{array}$

The beneficial effects of the present disclosure are as follows:

This method directly establishes an analytical solution for the dynamic model of natural gas system. Compared to traditional numerical methods based on discretization, this method completely avoids approximation errors, numerical dispersion, and dissipation. Additionally, by avoiding the discretization process during the solution, this method significantly improves the computational efficiency and solution accuracy of dynamic analysis of the natural gas system.

BRIEF DESCRIPTION OF THE DRAWINGS

The disclosure will now be further described with reference to the drawings.

FIG. 1 is a specific flow chart of the present disclosure;

FIG. 2 is a structural diagram of a natural gas system in an exemplary embodiment of the present disclosure; and

FIG. 3 is a calculation result graph of mass flow at a node in a natural gas system of the present disclosure.

DETAILED DESCRIPTION

The technical solutions in the embodiments of the present disclosure will be described clearly and completely below in combination with the drawings in the embodiments of the present disclosure. Clearly, the described embodiment is not all but only part of embodiments of the present disclosure. All other embodiments obtained by those ordinarily skilled in the art based on the embodiments in the present disclosure without creative work shall fall within the scope of protection of the present disclosure.

In the description of this specification, references to the terms “one embodiment”, “an example”, “a specific example”, and the like indicate that a specific feature, structure, material, or characteristic described in connection with the embodiment or example is included in at least one embodiment or example of the present disclosure. In this specification, schematic representations of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific feature, structure, material, or characteristic described may be combined in any suitable manner in any one or more embodiments or examples.

Embodiment: Taking the natural gas system shown in FIG. 2 as an example, the state distribution is studied under given initial conditions and boundary conditions.

As shown in FIG. 1, an embodiment of the present disclosure provides an analytical method for dynamic analysis of a natural gas system, including the following steps:

At step 10, a dynamic model of natural gas transmission is established according to the conservation equations and reconstructed into the equations in a heat conduction equation form.

At step 101, a dynamic model of natural gas transmission is established in a pipeline, including a mass conservation equation and a momentum conservation equation as below:

$\begin{array}{cc}\frac{\partial p_{x,t}}{\partial t}+\frac{c^2}{s}\frac{\partial q_{x,t}}{\partial x}=0,\frac{\partial p_{x,t}}{\partial t}+\frac{\lambda {\mathrm{wq}}_{x,t}}{2\mathrm{SD}}=0,& \left(1\right)\end{array}$

• • where x and t represent space and time, respectively, and c represents a speed of sound, taken as 340 m/s; p_x,t represents a pressure of natural gas, Pa; q_x,t represents mass flow of natural gas, kg/s; w represents an average velocity of natural gas, m/s; D represents an inner diameter of the pipeline, m; S represents a cross-sectional area of the pipeline, m²; and λ represents a resistance coefficient of the pipeline. In addition, initial conditions and boundary conditions are set in the dynamic analysis of the natural gas system as follows:

$\begin{array}{cc}{p}_{x,0}={\phi }_{p,x},q_{x,0}={\phi }_{q,x},p_{0,t},q_{L,t},& \left(2\right)\end{array}$

• • where p_x,0 and q_x,0 represent distributions of pressure and mass flow in the pipeline at moment 0 respectively, namely, the initial conditions; φ_p,x and φ_q,x represent initial condition functions of pressure and mass flow, respectively; L represents a length of the pipeline, m; p_0,t represents a time-series of pressure at an inlet of the pipeline, namely, pressure boundary; and q_L,t represents a time-series function of mass flow at an outlet of the pipeline, namely, flow boundary.

At step 102, a partial derivative for x is derived from a second equation in Eq. (1) before substituting into a first equation in Eq. (1) to obtain a heat conduction form equation as follows:

$\begin{array}{cc}\frac{\partial q_{x,t}}{\partial x}=\frac{-2\mathrm{DS}}{\lambda w}\frac{{\partial }^2{p}_{x,t}}{\partial x^2}\Rightarrow \{\begin{array}{c}\frac{\partial p_{x,t}}{\partial t}=\alpha \frac{{\partial }^2{p}_{x,t}}{\partial x^2}\\ \alpha =\frac{2{\mathrm{Dc}}^2}{\lambda w}\end{array},& \left(3\right)\end{array}$

• • where α is a coefficient term in the heat conduction form equation; accordingly, the partial derivative of pipeline outlet pressure for x can be obtained through the following equation and is defined as m_t:

$\begin{array}{cc}\frac{\partial p}{\partial x}{❘}_{x=L}=-\frac{\lambda w}{2\mathrm{DS}}{q}_{L,t}=m_t.& \left(4\right)\end{array}$

Eq. (1) may be converted to as follows:

$\begin{array}{cc}\{\begin{array}{c}\frac{\partial p_{x,t}}{\partial t}=\alpha \frac{{\partial }^2{p}_{x,t}}{\partial x^2},\frac{\partial p_{x,t}}{\partial t}{❘}_{x=L}=m_t\\ p\left(0,t\right)=p_{0,t},p\left(x,0\right)={\phi }_{p,x}\end{array}.& \left(5\right)\end{array}$

At step 20, the reconstructed dynamic model of natural gas transmission is transformed into a partial differential equation set with homogeneous boundary conditions using a method of variation of parameters, followed by deriving a general analytical expression for the dynamic model for natural gas transmission when the pressure boundary is constant.

At step 201, p_0,t is taken constant with the value of ψ_p,1 in response to the pressure boundary being constant, followed by introducing an intermediate variable u_x,t and reconstructing p_x,t as follows:

$\begin{array}{cc}{p}_{x,t}=u_{x,t}+x\left(\frac{\partial p}{\partial x}{❘}_{x=L}\right)+p_{0,t}.& \left(6\right)\end{array}$

• • x=0 and t=0 are substituted into Eq. (6) to obtain initial conditions and a set of boundary conditions for u_x,t, respectively:

$\begin{array}{cc}\{\begin{array}{c}{u}_{x,0}=p_{x,0}-x\left(\frac{\partial p}{\partial x}{❘}_{x=L}\right){❘}_{t=0}-p_{0,0}\\ p_{0,t}=u_{0,t}+p_{0,t}\Rightarrow u_{0,t}=0\end{array}.& \left(7\right)\end{array}$

A partial derivative for x is derived on both sides of Eq. (6) at x=L to obtain a further set of boundary conditions for u_x,t, expressed as follows:

$\begin{array}{cc}\frac{\partial p_{x,t}}{\partial x}=\frac{\partial u_{x,t}}{\partial x}+\frac{\partial p_{x,t}}{\partial x}{❘}_{x=L}\Rightarrow \frac{\partial u_{x,t}}{\partial x}{❘}_{x=L}=0.& \left(8\right)\end{array}$

A second-order partial derivative for x and a first-order partial derivative for t are derived on both sides of Eq. (6), respectively:

$\begin{array}{cc}\frac{{\partial }^2{p}_{x,t}}{\partial x^2}=\frac{\partial u_{x,t}^2}{\partial x^2},\frac{\partial p_{x,t}}{\partial t}=\frac{\partial u_{x,t}}{\partial t}+x\frac{\partial }{\partial t}\left(\frac{\partial p}{\partial x}{❘}_{x=L}\right).& \left(9\right)\end{array}$

Eq. (9) is substituted into Eq. (5), with Eq. (5) re-expressing a partial differential equation set with homogeneous boundary conditions as follows:

$\begin{array}{cc}\{\begin{array}{c}\frac{\partial u_{x,t}}{\partial t}=\alpha \frac{{\partial }^2{u}_{x,t}}{\partial x^2}+f_{x,t}\\ u_{x,0}=p_{x,0}-x\left(\frac{\partial p}{\partial x}{❘}_{x=L}\right){❘}_{t=0}-p_{,0,0},\frac{\partial u_{x,t}}{\partial x}{❘}_{x=L}=0\\ u_{0,t}=0,f_{x,t}=-x\frac{\partial }{\partial t}\left(\frac{\partial p}{\partial x}{❘}_{x=L}\right)=-x\frac{\partial m_t}{\partial t}\end{array}.& \left(10\right)\end{array}$

At step 202, u_x,t* is defined as a solution of a homogeneous problem of Eq. (10), and the homogeneous problem is solved as follows:

$\begin{array}{cc}{u}_{x,t}^{*}=T_0\left(t\right){\sum }_{n=1}^{\infty }{k}_2\sin \frac{\left(2n-1\right)\pi x}{2L},& \left(11\right)\end{array}$

• • where k₂ represents a coefficient term; T₀(t) represents a function component about time in u_x,t*; n represents the number of Fourier components; analogy to a form of u_x,t*, a form of u_x,t in Eq. (10) is as follows:

$\begin{array}{cc}{u}_{x,t}={\sum }_{n=1}^{\infty }{u}_{n,t}\sin \frac{\left(2n-1\right)\pi x}{2L},& \left(12\right)\end{array}$

• • where u_n,t represents a function component about time in u_x,t, f_x,t in Eq. (10) is expanded according to a base in Eq. (11):

$\begin{array}{cc}{f}_{x,t}={\sum }_{n=1}^{\infty }{f}_{n,t}\sin \frac{\left(2n-1\right)\pi x}{2L},& \left(13\right)\end{array}$ $\begin{array}{cc}{f}_{n,t}=\frac{2}{L}{\int }_0^L{f}_{x,t}\sin \frac{\left(2n-1\right)\pi x}{2L}\mathrm{dx},& \left(14\right)\end{array}$

• • where f_n,t represents a function component about time in f_x,t, obtained by inverse Fourier transformation shown in Eq. (14); Eq. (12) to (14) are substituted into a first equation in Eq. (10):

$\begin{array}{cc}\{\begin{array}{c}\frac{\partial u_{n,t}}{\partial t}+g_n{u}_{n,t}=f_{n,t}\\ g_n=\frac{{\left(2n-1\right)}^2{\pi }^2\alpha }{4L^2}\\ u_{n,0}=\frac{2}{L}{\int }_0^L{u}_{x,0}\sin \frac{\left(2n-1\right)\pi x}{2L}\mathrm{dx}\end{array},& \left(15\right)\end{array}$

• • where g_n represents a coefficient term in an n-th nonhomogeneous partial differential, u_n,0 represents a value of u_n,t when t=0; Eq. (15) is solved to obtain analytical expressions for u_n,t and u_x,t, respectively as follows:

$\begin{array}{cc}{u}_{n,t}=e^{-g_{n}t}\left\{u_{n,0}+{\int }_0^t{f}_{n,\tau }{e}^{g_n\tau }d\tau \right\},& \left(16\right)\end{array}$ $\begin{array}{cc}{u}_{x,t}={\sum }_{n=1}^{\infty }{e}^{-g_{n}t}\left\{u_{n,0}+{\int }_0^t{f}_{n,\tau }{e}^{g_n\tau }d\tau \right\}\sin \frac{\left(2n-1\right)\pi x}{2L}.& \left(17\right)\end{array}$

Eq. (17) is substituted into Eq. (6) to obtain an analytical expression for p_x,t at constant pressure:

$\begin{array}{cc}{p}_{x,t}={\sum }_{n=1}^{\infty }{e}^{-g_{n}t}\left\{u_{n,0}+{\int }_0^t{f}_{n,\tau }{e}^{g_n\tau }d\tau \right\}\sin \frac{\left(2n-1\right)\pi x}{2L}-\frac{\lambda wx}{2Ds}{q}_{L,t}+p_{0,t}.& \left(18\right)\end{array}$

Eq. (18) is substituted into a second equation in Eq. to obtain an analytical expression for q_x,t at constant pressure:

$\begin{array}{cc}{q}_{x,t}=-\frac{2DS}{\lambda w}\frac{\partial p_{x,t}}{\partial x}=q_{L,t}+\left\{{\sum }_{n=1}^{\infty }\frac{\left(1-2n\right)\pi DS}{\lambda wL{e}^{g_{n}t}}\left\{u_{n,0}+{\int }_0^t{f}_{n,\tau }{e}^{g_n\tau }d\tau \right\}\cos \frac{\left(2n-1\right)\pi x}{2L}\right\},& \left(19\right)\end{array}$

• • where p_0,t represents a constant-voltage boundary, taking a constant value.

Eq. (18) and Eq. (19) are general analytical expressions for the partial differential equation set of the dynamic model for natural gas transmission under constant pressure boundaries.

At step 30, a practical analytical expression for the dynamic model for natural gas transmission is constructed combining the superposition principle of natural gas transmission and according to the discrete form of the initial conditions and boundary conditions.

At step 301, first, the pressure boundary is reconstructed, expressed as:

$\begin{array}{cc}\{\begin{array}{c}{\psi }_p=\left[{\psi }_{p,1}{\psi }_{p,2},\dots ,{\psi }_{p,\mathrm{Nt}}\right]\\ {\psi }_{p,i}=\{\begin{array}{cc}{p}_{0,i}& i=1\\ p_{0,i}-p_{0,i-1}& i>1\end{array}\end{array},t_i=i\Delta t,& \left(20\right)\end{array}$

• • where Nt represents a length of the boundary condition, and Δt is a time step; ψ_p represents a new boundary condition corresponding to a pressure increment, and an element ψ_p,i represents an increment of a pressure at moment t_i for a pressure at moment t_i−1.

At step 302, according to superposition principle of natural gas transmission, the pressures p_x,t is decomposed as follows:

$\begin{array}{cc}{p}_{x,t}={\sum }_{k=1}^i{\chi }_{k,x,t-t_{k-1}}1\le i\le Nt,t_{k-1}\le t\le t_i,& \left(21\right)\end{array}$

• • where each component χ_k,x,t satisfies:

$\begin{array}{cc}\{\begin{array}{cc}\frac{\partial {\chi }_{k,x,t-t_{k-1}}}{\partial t}=\alpha \frac{{\partial }^2{\chi }_{k,x,t-t_{k-1}}}{\partial x^2},& t>t_{k-1},0<x<L\\ {\frac{\partial {\chi }_{k,x,t-t_{k-1}}}{\partial t}❘}_{x=L}=\{\begin{array}{cc}{m}_t& k=1\\ 0& k>1\end{array},& t>t_{k-1}\\ {{\chi }_{k,x,t-t_{k-1}}❘}_{x=0}={\psi }_{p,k},& t>t_{k-1}\\ {{\chi }_{k,x,t-t_{k-1}}❘}_{t=t_{k-1}}=\{\begin{array}{cc}{\phi }_{p,x}& k=1\\ 0& k>1\end{array},& 0\le x\le L\end{array}.& \left(22\right)\end{array}$

When k=1, a form of χ_k,x,t is consistent with Eq. (18), namely:

$\begin{array}{cc}{\chi }_{k,x,t-t_{k-1}}={\sum }_{n=1}^{\infty }{e}^{-g_{n}t}\left\{u_{n,0}+{\int }_0^t{f}_{n,\tau }{e}^{g_n\tau }d\tau \right\}\sin \frac{\left(2n-1\right)\pi x}{2L}+x{m}_t+{\psi }_{p,1}k=1,t>t_{k-1}.& \left(23\right)\end{array}$

When k>1, a form of χ_k,x,t is as follows:

$\begin{array}{cc}{\chi }_{k,x,t-t_{k-1}}={\psi }_{p,k}k>1,t>t_{k-1}.& \left(24\right)\end{array}$

Eq. (23) and Eq. (24) are substituted into Eq. (21), it may be found that at any pressure boundary, the analytical expressions for p_x,t and q_x,t are still as shown in Eq. (18) and Eq. (19) with the difference that p_0,t is transformed from a constant to a set of discrete sequences.

At step 303, the initial conditions are calculated according to steady-state energy flow in practice, expressed as follows:

$\begin{array}{cc}{\phi }_{p,x}=K_{1}x+K_2,{\phi }_{q,x}=K_3,m_0=\frac{-\lambda w{\phi }_{q,L}}{2Ds},& \left(25\right)\end{array}$

• • where K₁, K₂, and K₃ represent constant coefficients in the initial conditions; m₀ represents a value of m_t at t=0; and φ_q,L represents a value of the initial conditions of mass flow at x=L, namely, the mass flow at an outlet of the pipeline at t=0. m_t is expressed as a segmented step function as follows:

$\begin{array}{cc}{m}_t={\sum }_{i=1}^{\mathrm{Nt}}{m}_i\left({\delta }_{t-t_i}-{\delta }_{t-t_{i-1}}\right)t>0.& \left(26\right)\end{array}$

f_x,t in Eq. (10) and f_n,t in Eq. (14) are substituted into Eq. (16) as follows:

$={\int }_0^t\left(\frac{-2}{L}{\int }_0^{L}x\frac{\partial m_{\tau }}{\partial \tau }\sin \frac{\left(2n-1\right)\pi x}{2L}dx\right)e^{g_n^{\tau }}d\tau$ $\begin{array}{cc}{\int }_0^t{f}_{n,\pi }{e}^{g_n^{\tau }}d\tau =\frac{-2}{L}{\int }_0^{L}x\sin \frac{\left(2n-1\right)\pi x}{2L}dx{\int }_0^t\frac{\partial m_{\tau }}{\partial \tau }{e}^{g_n^{\tau }}d\tau .& \left(27\right)\end{array}$ $=\frac{8{L\left(-1\right)}^n}{{\left(2n-1\right)}^2{\pi }^2}{\int }_0^t\frac{\partial m_{\tau }}{\partial \tau }{e}^{g_n^{\tau }}d\tau$

An integral formula in Eq. (27) may be modified as follows:

$\begin{array}{cc}\frac{\partial m_{\tau }}{\partial \tau }{e}^{g_n\tau }=\frac{\partial \left(m_{\tau }{e}^{g_n\tau }\right)}{\partial \tau }-g_n{m}_{\tau }{e}^{g_n\tau },& \left(28\right)\end{array}$ $\begin{array}{cc}{g}_n{\int }_0^{t_k}{m}_{\tau }{e}^{g_n\tau }d\tau ={\sum }_{i=1}^k{m}_i\left(e^{g_n{t}_i}-e^{g_n{t}_{i-1}}\right).& \left(29\right)\end{array}$

Eq. (28) and Eq. (29) are substituted into Eq. (27) as follows:

$\begin{array}{cc}\frac{1}{e^{g_n{t}_k}}{\int }_0^{t_k}{f}_{n,\tau }{e}^{g_n\tau }d\tau =-E_{1n}\left(m_k-\frac{m_0}{e^{g_n{t}_k}}-{\sum }_{i=1}^k{m}_i\frac{e^{g_n\Delta t}-1}{e^{g_n\left(t_k-t_{i-1}\right)}}\right),& \left(30\right)\end{array}$ $\begin{array}{cc}{E}_{1n}=\frac{8{L\left(-1\right)}^{n-1}}{{\left(2n-1\right)}^2{\pi }^2}.& \left(31\right)\end{array}$

Eq. (30) is substituted into Eq. (18) and Eq. (19), to obtain the practical analytical expressions for the dynamic model for natural gas transmission as follows:

$\begin{array}{cc}{p}_{x,k}={\sum }_{n=1}^{\infty }\left\{\frac{u_{n,0}}{e^{g_n{t}_k}}-E_{1n}\left(m_k-\frac{m_0}{e^{g_n{t}_k}}-{\sum }_{i=1}^k{m}_i\frac{e^{g_n\Delta t}-1}{e^{g_n\left(t_k-t_{i-1}\right)}}\right)\right\}\sin \frac{\left(2n-1\right)\pi x}{2L}-\frac{\lambda wx}{2DS}{q}_{L,k}+p_{0,k},& \left(32\right)\end{array}$ $\begin{array}{cc}{q}_{x,k}=q_{L,k}+\left\{{\sum }_{n=1}^{\infty }\frac{\left(1-2n\right)\pi DS}{\lambda wL}\left\{\frac{u_{n,0}}{e^{g_n{t}_k}}-E_{1n}\left(m_k-\frac{m_0}{e^{g_n{t}_k}}-{\sum }_{i=1}^k{m}_i\frac{e^{g_n\Delta t}-1}{e^{g_n\left(t_k-t_{i-1}\right)}}\right)\right\}\cos \frac{\left(2n-1\right)\pi x}{2L}\right\}.& \left(33\right)\end{array}$

At step 40, the dynamic model of natural gas transmission is reconstructed; and the natural gas network equivalent model is derived based on the analytical method, as well as the simplified analysis method thereof combining topological characteristics, the time-invariant and causal properties of natural gas transmission.

At step 401, x=L and x=0 are substituted into Eq. (32) and Eq. (33), respectively, to obtain pressure at the outlet and flow at the inlet of the pipeline at arbitrary t_j as follows:

$\begin{array}{cc}{p}_{L,j}={\sum }_{n=1}^{\infty }{e}^{-g_n{t}_j}{u_{n,0}\left(-1\right)}^{n-1}-\frac{\lambda wL}{2Ds}{q}_{L,j}+p_{0,j}+{\sum }_{n=1}^{\infty }\frac{\lambda w{E}_{1n}}{2Ds}{\left(-1\right)}^{n-1}\left(q_{L,j}-\frac{q_{L,0}}{e^{g_n{t}_j}}-{\sum }_{i=1}^j{q}_{L,i}\frac{e^{g_n\Delta t}-1}{e^{g_n\left(j-i+1\right)\Delta t}}\right),& \left(34\right)\end{array}$ $\begin{array}{cc}{q}_{0,j}=q_{L,j}+\frac{Ds\pi }{\lambda wL}{\sum }_{n-1}^{\infty }\left(1-2n\right)u_{n,0}{e}^{-g_n{t}_j}+\frac{\pi }{2L}{\sum }_{n-1}^{\infty }\left(1-2n\right)E_{1n}\left(q_{L,j}-\frac{q_{L,0}}{e^{g_n{t}_j}}-{\sum }_{i=1}^j{q}_{L,i}\frac{e^{g_n\Delta t}-1}{e^{g_n\left(j-i+1\right)\Delta t}}\right).& \left(35\right)\end{array}$

Eq. (34) and Eq. (35) are reorganized as follows:

$\begin{array}{cc}\left[\begin{array}{c}{p}_L\\ q_0\end{array}\right]=\left[\begin{array}{cc}{H}_1& H_2\\ H_3& H_4\end{array}\right]\left[\begin{array}{c}{p}_0\\ q_L\end{array}\right]+\left[\begin{array}{c}{\gamma }_p\\ {\gamma }_q\end{array}\right],& \left(36\right)\end{array}$

• • where p_L and p₀ represent a pressure vector at the outlet and inlet of the pipeline, respectively; q_L and q₀ represent a flow vector at the outlet and inlet of the pipeline, respectively; H_1-4 represent constant coefficient matrices with the dimension of N_t×N_t; and γ_p and γ_q represent a constant coefficient vector of pressure and mass flow with the dimension of N_t×1, respectively, determined by the initial conditions; the elements may be represented as follows:

$\begin{array}{cc}\{\begin{array}{cc}{p}_0=\left[p_{0,1},\cdots ,p_{0,\mathrm{Nt}}\right],& p_L=\left[p_{L,1},\cdots ,p_{L,\mathrm{Nt}}\right]\\ q_L=\left[q_{L,1},\cdots ,q_{L,\mathrm{Nt}}\right],& q_0=\left[q_{0,1},\cdots ,q_{0,\mathrm{Nt}}\right]\end{array},& \left(37\right)\end{array}$ $\begin{array}{cc}{\gamma }_{p,j}={\sum }_{n=1}^{\infty }\left(u_{n,0}+E_{1n}{m}_0\right){\left(-1\right)}^{n-1}{e}^{-g_n{t}_j},& \left(38\right)\end{array}$ $\begin{array}{cc}{\gamma }_{q,j}=\frac{\mathrm{DS}\pi }{\lambda \mathrm{wL}}{\sum }_{n-1}^{\infty }\left(1-2n\right)\left(u_{n,0}+E_{1n}{m}_0\right)e^{-g_n{t}_j},& \left(39\right)\end{array}$ $\begin{array}{cc}{H}_1=1,H_3=0,\underset{i>j}{H_{2.\mathrm{ji}}}=\underset{i>j}{H_{4.\mathrm{ji}}}=0,& \left(40\right)\end{array}$ $\begin{array}{cc}\underset{i=j}{H_{2,\mathrm{ji}}}=\frac{\lambda w}{2\mathrm{DS}}\left({\sum }_{n=1}^{\infty }{E}_{1n}{\left(-1\right)}^{n-1}{e}^{-g_n\Delta t}-L\right),& \left(41\right)\end{array}$ $\begin{array}{cc}\underset{\left(i<j\right)}{H_{2,\mathrm{ji}}}=\frac{\lambda w}{2\mathrm{DS}}\left({\sum }_{n=1}^{\infty }{E}_{1n}{\left(-1\right)}^{n-1}\frac{1-e^{g_n\Delta t}}{e^{g_n\left(j-i+1\right)\Delta t}}\right),& \left(42\right)\end{array}$ $\begin{array}{cc}\underset{\left(i=j\right)}{H_{4,\mathrm{ij}}}=\frac{\pi }{2L}{\sum }_{n=1}^{\infty }{E}_{1n}\left(1-2n\right)e^{-g_n\Delta t}+1,& \left(43\right)\end{array}$ $\begin{array}{cc}\underset{\left(i<j\right)}{H_{4,\mathrm{ji}}}=\frac{\pi }{2L}{\sum }_{n=1}^{\infty }{E}_{1n}\left(1-2n\right)\frac{I-e^{g_n\Delta t}}{e^{g_n\left(j-i+1\right)\Delta t}}.& \left(44\right)\end{array}$

At step 402, the number of nodes and pipelines in the natural gas system are defined as N_g and N_b, respectively, and an identity matrix and a zero matrix are defined as 1 and 0, respectively; incidence matrices are introduced to describe topological properties in the natural gas system.

(1) A_in includes N_g×N_b sub-matrices for associating flow at end of a branch with flow at a node, and in response to a sub-matrix of an i-th row and j-th column being 1, the flow flows from the end of the branch j into the node i, otherwise, the sub-matrix is 0.

(2) A_out includes N_g×N_b sub-matrices for associating flow at head end of a branch with flow at a node, and in response to a sub-matrix of an i-th row and j-th column being 1, the flow flows from the node i into the head end of the branch j, otherwise, the sub-matrix is 0.

(3) A_nb1 includes N_b×N_g sub-matrices for associating pressure at end of a branch with pressure at a node, and in response to a sub-matrix of an i-th row and j-th column being 1, pressure at the node j is equal to pressure at the end of the branch i, otherwise, the sub-matrix is 0.

(4) A_nb2 includes N_b×N_g sub-matrices for associating pressure at a head end of a branch with pressure at a node, and in response to a sub-matrix of an i-th row and j-th column being 1, pressure at the node j is equal to pressure at the head end of the branch i, otherwise, the sub-matrix is 0; particularly, in response to the node j being a compressor node with a pressure ratio of K_cp and the node j being connected to the head end of the branch i, the sub-matrix of the i-th row and jth column is K_cp×1.

At step 403, equations of mass conservation at nodes and pressure continuity between nodes and pipelines are constructed based on the incidence matrices, expressed as follows:

$\begin{array}{cc}{A}_{\mathrm{in}}{q}_L-A_{\mathrm{out}}{q}_0=\mathrm{qn},& \left(45\right)\end{array}$ $\begin{array}{cc}{p}_L=A_{\mathrm{nb}1}\mathrm{pn},p_0=A_{\mathrm{nb}2}\mathrm{pn},& \left(46\right)\end{array}$

• • where pn and qn represent a pressure and a vector at the node, respectively; p_L, p₀, q_L, and q₀ are all vectors of Nt×1; the sub-matrices 1 and 0 in coefficient matrices A_in, A_out, A_nb1, and A_out2 are all Nt×Nt matrices.

Eq. (46) is substituted into a first row of Eq. (36) as follows:

$\begin{array}{cc}{q}_L=H_2^{-1}\left(\left(A_{\mathrm{nb}1}-A_{{\mathrm{nb}}_2}\right)\mathrm{pn}-{\gamma }_p\right)=H_{a1}\mathrm{pn}+H_{a2},& \left(47\right)\end{array}$

• • where H_a1 and H_a2 represent a constant coefficient matrix and a vector describing mass flow at an outlet of a pipeline, respectively; Eq. (46) and Eq. (47) are substituted into a second row of Eq. (36) as follows:

$\begin{array}{cc}{q}_0=H_4\left(H_{a1}\mathrm{pn}+H_{a2}\right)+{\gamma }_q=H_{a_3}\mathrm{pn}+H_{a_4},& \left(48\right)\end{array}$

• • where H_a3 and H_a4 represent a constant coefficient matrix and a vector describing mass flow at an inlet of the pipeline, respectively; Eq. (47) and Eq. (48) are substituted into Eq. (45) to associate pressure and mass flow at the node as follows:

$\begin{array}{cc}\mathrm{qn}\begin{array}{c}=\left(A_{\mathrm{in}}{H}_{a1}-A_{\mathrm{out}}{H}_{a3}\right)pn+\left(A_{\mathrm{in}}{H}_{a2}-A_{\mathrm{out}}{H}_{a4}\right)\\ =H_{a5}\mathrm{pn}+H_{a6}\end{array},& \left(49\right)\end{array}$

• • where H_a5 and H_a6 are a constant coefficient matrix and a vector describing a relationship between pressure and mass flow at a node, respectively; subscripts sr and ns are defined as a variable identification of a gas source node and a gas load node, respectively, and Eq. (49) is expanded to obtain the natural gas network equivalent model based on the practical analytical expressions:

$\begin{array}{cc}\left[\begin{array}{c}{\mathrm{qn}}_{\mathrm{sr}}\\ {\mathrm{qn}}_{\mathrm{ns}}\end{array}\right]=\left[\begin{array}{cc}{H}_{a51}& H_{a52}\\ H_{a53}& H_{a54}\end{array}\right]\left[\begin{array}{c}{\mathrm{pn}}_{\mathrm{sr}}\\ {\mathrm{pn}}_{\mathrm{ns}}\end{array}\right]+\left[\begin{array}{c}{H}_{a61}\\ H_{a62}\end{array}\right],& \left(50\right)\end{array}$

• • where H_a51, H_a52, H_a53, and H_a54 are sub-matrices in H_a5, and H_a61 and H_a62 are sub-vectors in H_a6; pn_sr and pn_ns are a pressure vector at a source node and a load node, respectively; qn_sr and qn_ns are a mass flow vector at the source node and the load node, respectively; the mass flow at the source node and the pressure at the load node in the natural gas system are fully expressed as functions of the boundary conditions using matrix transformation, namely, a simplified method of natural gas network equivalent model is as follows:

$\begin{array}{cc}{\mathrm{pn}}_{\mathrm{ns}}=H_{a54}^{-1}\left({\mathrm{qn}}_{\mathrm{ns}}-H_{a53}{\mathrm{pn}}_{\mathrm{sr}}-H_{a62}\right),& \left(51\right)\end{array}$ $\begin{array}{cc}{\mathrm{qn}}_{\mathrm{sr}}=H_{a51}{\mathrm{pn}}_{\mathrm{sr}}+H_{a52}{\mathrm{pn}}_{\mathrm{ns}}+H_{a61}.& \left(52\right)\end{array}$

The pressures at load nodes and mass flows at source nodes are calculated using the analytical functions as shown in FIG. 3.

The present disclosure proposes an analytical method for dynamic analysis of a natural gas system, which is based on a dynamic model of the natural gas system. Furthermore, based on the superposition principle, time-invariant and causal properties of natural gas transmission, practical analytical forms of the dynamic model under different initial conditions and boundary conditions, as well as the equivalent network model, are constructed, reducing the computational complexity and greatly improving the computational accuracy.

The foregoing has shown and described the basic principles, principal features, and advantages of the present disclosure. It is to be understood by those skilled in the art that the present disclosure is not limited to the above embodiments. The above embodiments and description in the present disclosure are merely illustrative of the principles of the present disclosure. Without deviating from the spirit and scope of the present disclosure, the present disclosure is subject to various changes and improvements, and these changes and improvements fall within the scope of the present disclosure that requires protection.

Claims

1. An analytical method for dynamic analysis of a natural gas network, comprising: establishing a dynamic model of natural gas transmission according to the conservation equations, and reconstructing the dynamic model into the equations in a heat conduction equation form;transforming the heat conduction equation into a partial differential equation set with homogeneous boundary conditions using a method of variation of parameters, and constructing a general analytical expression for the partial differential equation set under constant pressure boundaries;constructing, based on superposition principle, practical analytical expressions for the dynamic model of natural gas transmission;constructing, according to system topologies and time-invariant, and causal properties of a natural gas system, a natural gas network equivalent model based on the practical analytical expressions; anddetermining pressure at load nodes and mass flows at source nodes based on the natural gas network equivalent model.

2. The analytical method for dynamic analysis of a natural gas network according to claim 1, wherein the establishing a dynamic model of natural gas transmission according to the conservation equations comprises the following steps: establishing a dynamic model of natural gas transmission in a pipeline, comprising a mass conservation equation and a momentum conservation equation as below:

3. The analytical method for dynamic analysis of a natural gas network according to claim 2, wherein the reconstructing the dynamic gas flow model into the equations in a heat conduction equation form comprises the following steps: deriving a partial derivative for x from a second equation in Eq. before substituting into a first equation in Eq. (1) to obtain a heat conduction form equation as follows:

4. The analytical method for dynamic analysis of a natural gas network according to claim 3, wherein the transforming the heat conduction equation into a partial differential equation set with homogeneous boundary conditions using a method of variation of parameters comprises the following steps:taking, by p0,t, a constant value ψp,1 in response to the pressure boundary being constant; introducing an intermediate variable ux,t and reconstructing px,t as follows:

5. The analytical method for dynamic analysis of a natural gas network according to claim 4, wherein the constructing a general analytical expression for the partial differential equation set under constant pressure boundaries comprises the following steps: defining ux,t* as a solution of a homogeneous problem of Eq. (10), and solving the homogeneous problem as follows:

6. The analytical method for dynamic analysis of a natural gas network according to claim 1, wherein the constructing, based on superposition principle, practical analytical expressions for the dynamic model of natural gas transmission comprises the following steps: reconstructing the pressure boundary, expressed as follows:

7. The analytical method for dynamic analysis of a natural gas network according to claim 1, wherein the constructing, according to system topologies, and time-invariant, and causal properties of a natural gas system, a natural gas network equivalent model based on the practical analytical expressions comprises the following steps:substituting x=L and x=0 into Eq. (32) and Eq. (33), respectively, to obtain pressure at the outlet and flow at the inlet of the pipeline at arbitrary tj as follows: