Applications of Symmetry-Enhanced Physics-Informed Neural Networks in High-Pressure Gas Flow Simulations in Pipelines

Abstract

This article presents a detailed examination of the methodology and modeling tools utilized to analyze gas flows in pipelines, rooted in the fundamental principles of gas dynamics. The methodology integrates numerical simulations with modern neural network techniques, particularly focusing on the PINN utilizing the continuous symmetry data inherent in PDEs, which is called the symmetry-enhanced Physics-Informed Neural Network. This innovative approach combines artificial neural networks (ANNs) integrating physical equations, which provide enhanced efficiency and accuracy when modeling various complex processes related to physics with a symmetric and asymmetric nature. The presented mathematical model, based on the system of Euler equations, has been carefully implemented using Python language. Verification with analytical solutions ensures the accuracy and reliability of the computations. In this research, a comparative and comprehensive analysis was carried out comparing the outcomes obtained using the symmetry-enhanced PINN method and those from conventional computational fluid dynamics (CFD) approaches. The analysis highlighted the advantages of the symmetry-enhanced PINN method, which produced smoother pressure and velocity fluctuation profiles while reducing the computation time, demonstrating its capacity as a revolutionary modeling tool. The estimated results derived from this study are of paramount importance for ensuring ongoing energy supply reliability and can also be used to create predictive models related to gas behavior in pipelines. The application of modeling techniques for gas flow simulations has the potential to improve the integrity of our energy infrastructure and utilization of gas resources, contributing to advancing our understanding of symmetry principles in nature. However, it is crucial to emphasize that the effectiveness of such models relies on continuous monitoring and frequent updates to ensure alignment with real-world conditions. This research not only contributes to a deeper understanding of compressible gas flows but also underscores the crucial role of advanced modeling methodologies in the sustainable management of gas resources for both current and future generations. The numerical data covered the physics of the process related to the modeling of high-pressure gas flows in pipelines with regard to density, velocity and pressure, where the PINN model was able to outperform the classical CFD method for velocity by 170% and for pressure by 360%, based on $L_{\infty }$ values.

1. Introduction

The efficient transportation and distribution of natural gas play a pivotal role in the economic development and energy security of nations worldwide [1]. In Kazakhstan, a country abundant in oil and gas resources, the modeling of the gas flow in tubes holds significant importance across various sectors. Understanding and accurately predicting the behavior of gas within pipelines and distribution networks are paramount for optimizing infrastructure, ensuring energy supply reliability and promoting environmental sustainability [2,3]. Thus, modeling gas flow in tubes holds critical significance for Kazakhstan, impacting its oil and gas industry, energy sector development, infrastructure planning, environmental situation, research and innovation endeavors, and safety assurance. Through advanced machine learning and computational fluid dynamics simulations, the country can optimize its energy resources, promote sustainable development, and bolster resilience in its energy infrastructure.

The comprehensive modeling of gas flow in tubes relies fundamentally on solving the partial differential equations (PDEs) governing fluid dynamics. PDEs, such as the Euler equations, form the mathematical framework for describing the conservation laws of mass, momentum, and energy in fluid flows. Solving these equations numerically using advanced computational techniques enables the accurate prediction of gas behavior, including shock wave propagation, turbulence effects, and flow characteristics within pipelines and tubes. In recent years, significant progress has been made in computational fluid dynamics (CFD) [4,5] for numerically solving PDEs, particularly in solving the Navier–Stokes equations [6,7,8], which has revolutionized the ability to simulate and analyze complex fluid flow phenomena using techniques like finite difference (FD) [9,10,11] and finite volume (FV) methods [12,13]. Even though models such as the Navier–Stokes equations [14,15,16,17,18] are able to describe the underlying problem, at this point, modeling still requires a high computation cost to achieve a high precision.

The rapid growth in data, the parallel emergence of computing, and the advent of graphics processors (GPUs), together with the advanced theoretical results in numerical analysis, have resulted in the explosive growth of machine learning (ML) [19,20], particularly in the area of physical problems and mathematical physics. Physics-Informed Neural Networks (PINNs) [21,22,23,24,25] are a type of artificial neural network that can dramatically reduce the computational complexity of modeling physical processes. The further need for developing the potential of these models is still relevant and they are still not a uniform method of CFD, but in the field of forward problems, they have a higher precision and a promising future. In recent advancements, PINNs have been employed to solve nuclear physics problems [26], astrophysics [27], and compressible flows governed by the Navier–Stokes equations [28,29], offering efficient simulations for different aspects of nature. However, despite these advancements, achieving optimal results in compressible flows remains a difficult process and requires skill.

This research considers a transient gas flow in pipelines using symmetry-enhanced Physics-Informed Neural Networks, with manual selection of the necessary hyperparameters using a weighted loss function [30]. The PINN leverages the continuous Lie symmetry information inherent in PDEs, wherein the invariant surface conditions (ISCs) induced by these symmetries are integrated into the PINN’s loss function [31]. When PDEs, along with their initial and boundary conditions, possess a Lie symmetry, the solutions remain invariant under such symmetry and also satisfy the ISCs. These ISCs, being intrinsic properties of the PDEs, impose crucial additional constraints on the solutions. Since the standard loss function in the PINN is the mean square error of the PDE residual, adding the ISCs to the loss function will definitely increase the number of constraints for the objective function during the optimization process and improve the accuracy of the solutions. In addition to this, after an evaluation of the PINN model’s performance, the following problems will be solved to demonstrate the effectiveness of the method: Burgers’ equation [32,33], viscous Burgers’ equation [33,34] and Euler’s equations of gas dynamics [35,36], where the main focus will be on the compressible flow problem.

The article is organized as follows. Section 2 introduces the physical problem in the general form of a PDE with initial and boundary conditions. The following Section 3 illustrates the numerical continuous PINN model with a forward propagation procedure for predictions of true solutions and the multi-objective loss for optimization, in addition to the architecture of the model. In addition to this, the discrete CFD model, which applies the four-step Runge–Kutta method with a general form and space derivative discretization. In Section 4, one-dimensional problems are solved with short discussions of their mathematical form: a PDE with initial and boundary conditions and physical descriptions with graphical comparisons to analytical solutions and L norm error measures. It is also worth noting the modeling approach in comparison with other works; for example, the authors of [37] use the sequence-to-sequence (seq2seq) technique, which splits the domains into parts, while in our case, the model is trained on all time segments, which will be shown to be a valid approach.

2. Physical Problem

The physical problem involves the gas flow in a pipe, which is described by a system of Euler’s equations. Equations govern the motion of adiabatic, inviscid, and compressible fluids, providing a powerful model for analyzing fluid dynamics in many practical scenarios. In order to capture the speed of the rarefaction wave, the discontinuity and the shock contacts, the problem is considered in a chamber that is separated into regions with high pressure and low pressure by interface points at $x=0.5$, as shown in Figure 1. The physical process is observed for the time domain $t\in [0,0.1]$. The space domain is defined as $x\in [0,1]$.

Euler’s system of equations represents the conservation of mass, momentum, and energy. The conservation of mass in the system is described by the continuity equation. It determines how changes in the density affect the fluid flow and ensures that the total mass within a fluid volume remains constant over time:

$$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u\right)}{\partial x}=0,$$

(1)

where the density of gas $\rho \left[\mathrm{kg}·{\mathrm{m}}^{-3}\right]$ changes with respect to time and space as influenced by the flow velocity component $u\left[\mathrm{m}·{\mathrm{s}}^{-1}\right]$.

The conservation of momentum is based on Newton’s laws of motion:

$$\frac{\partial \left(\rho u\right)}{\partial t}+\frac{\partial (\rho u^2+p)}{\partial x}=0.$$

(2)

The momentum equation describes how the velocity of the fluid changes in response to convection and pressure $p\left[\mathrm{kg}·{\mathrm{m}}^{-1}·{\mathrm{s}}^{-2}\right]$. It influences the formation of shocks and determines the acceleration of the gas flow.

The general form of the energy equation is derived from the principles of the conservation of energy, which is expressed as:

$$\frac{\partial E}{\partial t}+\frac{\partial \left[u\right(E+p\left)\right]}{\partial x}=0.$$

(3)

The energy equation describes the conservation of energy in fluid flows; it influences the system by determining how energy is transferred through advection and how work is done by pressure forces, where $E\left[\mathrm{J}·{\mathrm{m}}^{-3}\right]$ is the total energy and is defined as,

$$E=\frac{p}{\gamma -1}+\frac{1}{2}\rho u^2,$$

(4)

with a specific heat ratio $\gamma \left[-\right]$ equal to $1.4$.

Taking into account shock phenomenon generation, the following initial conditions are used:

$$\begin{array}{c}\hfill {(\rho ,u,p)|}_{t=0}=\left\{\begin{array}{c}(1,0,1),0\le x\le 0.5,\hfill \\ (0.125,0,0.1),0.5\le x\le 1.\hfill \end{array}\right.\end{array}$$

(5)

The Dirichlet boundary conditions are applied as boundary conditions, meaning that the quantities take on the values prescribed by the initial conditions at both boundaries, which ensures that the flow at the boundaries remains consistent with the initial conditions.

The transient gas flow in pipelines can be expressed in general form as:

$$\begin{array}{cc}\hfill & \frac{\partial U}{\partial t}+\mathcal{N}[U;\lambda ]=0,x\in \mathsf{\Omega },t\in [0,t_f],\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill & U(x,0)=\mathcal{I}\left(x\right),x\in \mathsf{\Omega },\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill & U(x,t)=\mathcal{B}\left(t\right),x\in \partial \mathsf{\Omega },t\in [0,t_f].\hfill \end{array}$$

(8)

The solution is denoted by $\widehat{U}(x,t)\in \mathbb{R}$ with spatial coordinate $x\in \mathbb{R}$, time t and parameter of the PDE $\lambda$, where $\mathsf{\Omega }$ is the computational domain and $\partial \mathsf{\Omega }$ is the boundary. Here, $\mathcal{N}(·)$ is a linear or nonlinear differential operator, the initial condition operator is $\mathcal{I}(·)$ and the boundary condition operator is $\mathcal{B}(·)$. Note that this form of representation of a parametrized PDE can be used for other problems as well.

3. Methodology

3.1. Continuous PINN Time Model

The continuous symmetry-enhanced PINN time model is used as framework for solving PDEs; the left-hand side of the equation is defined as,

$$f(x,t)=\frac{\partial U}{\partial t}+\mathcal{N}[U;\lambda ],$$

(9)

and the fully connected deep L neural network with inputs $d_{in}$ and outputs $d_{out}$ is defined as,

$$\mathcal{G}:{\mathbb{R}}^{d_{in}}\to {\mathbb{R}}^{d_{out}},$$

(10)

The forward propagation procedure will be used to approximate solution $\widehat{U}(x,t)$ as follows,

$$\begin{array}{cc}\hfill & {\mathcal{G}}^{\left[l\right]}=(x,t)\in {\mathbb{R}}^{d_{in}},l=0\hfill \\ \hfill & {\mathcal{G}}^{\left[l\right]}=\sigma (W^{\left[l\right]}{\mathcal{G}}^{[l-1]}+b^{\left[l\right]})\in {\mathbb{R}}^{n_l},1\le l\le L-1,\hfill \\ \hfill & {\mathcal{G}}^{\left[l\right]}=W^{\left[l\right]}{\mathcal{G}}^{[l-1]}+b^{\left[l\right]}\in {\mathbb{R}}^{d_{out}},l=L,\hfill \end{array}$$

(11)

Here, activation value ${\mathcal{G}}^{\left[l\right]}$ denotes the output of the $l^{th}$ layer with weight matrix $W^{\left[l\right]}\in {\mathbb{R}}^{n_l\times n_{l-1}}$, bias vector $b^{\left[l\right]}\in {\mathbb{R}}^{n_l}$ and nonlinear activation function $\sigma (·)$, where $n_l$ denotes the number of neurons in the $l^{th}$ layer. For the last layer L, a linear activation function with $d_{out}$ outputs is used to approximate the solution, as $\widehat{U}(x,t)\approx U(x,t)=W^{\left[L\right]}{\mathcal{G}}^{[L-1]}+b^{\left[L\right]}$.

The construction that converts the PDE into an optimization problem is the combined multi-objective weighted squared loss $\mathcal{L}$,

$$\mathcal{L}=(1-\alpha )({\mathcal{L}}_{\mathcal{I}}+{\mathcal{L}}_{\mathcal{B}})+\alpha {\mathcal{L}}_{\mathcal{F}},0\le \alpha \le 1,$$

(12)

where the individual loss terms are defined as:

$$\begin{array}{cc}\hfill & {\mathcal{L}}_{\mathcal{I}}=\frac{1}{N_I}\sum _{i=1}^{N_I}{(\mathcal{I}\left(x_I^i\right)-U(x_I^i,t_I^i))}^2,\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill & {\mathcal{L}}_{\mathcal{B}}=\frac{1}{N_B}\sum _{i=1}^{N_B}{(\mathcal{B}\left(t_B^i\right)-U(x_B^i,t_B^i))}^2,\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill & {\mathcal{L}}_{\mathcal{F}}=\frac{1}{N_F}\sum _{i=1}^{N_F}f{(x_F^i,t_F^i)}^2,\hfill \end{array}$$

(15)

Weights and biases denoted by $\theta$ can be trained by minimizing,

$${\theta }^{\ast }=arg\underset{\theta }{min}\mathcal{L}\left(\theta \right),$$

(16)

This minimization process is performed via the gradient-based ADAM Algorithm 1.

Algorithm 1: ADAM Optimization Algorithm

1: Input: $\gamma$ (learning rate lr), ${\beta }_1=0.9,{\beta }_2=0.999$ (coefficients), $ϵ={10}^{-8}$ (numerical stability term), ${\theta }_0$ (parameters), $\mathcal{L}\left(\theta \right)$ (objective) 2: Initialize: $m_0\leftarrow 0$ (first moment), $v_0\leftarrow 0$ (second moment) 3: for $i=1$ to … do 4:     $g_i\leftarrow {\nabla }_{\theta }{\mathcal{L}}_i\left({\theta }_{i-1}\right)$ 5:     $m_i\leftarrow {\beta }_1{m}_{i-1}+(1-{\beta }_1)g_i$ 6:     $v_i\leftarrow {\beta }_2{v}_{i-1}+(1-{\beta }_2)g_i^2$ 7:     ${\widehat{m}}_i\leftarrow m_i/(1-{\beta }_1^i)$ 8:     ${\widehat{v}}_i\leftarrow v_i/(1-{\beta }_2^i)$ 9:     ${\theta }_i\leftarrow {\theta }_{i-1}-\gamma {\widehat{m}}_i/(\sqrt{{\widehat{v}}_i}+ϵ)$ 10: end for

Here, to train the initial conditions, ${\{x_I^i,t_I^i\}}_{i=1}^{i=N_I}$ is used as input and ${\left\{\mathcal{I}\left(x_I^i\right)\right\}}_{i=1}^{i=N_I}$ is used as output; for the boundary conditions, ${\{x_B^i,t_B^i\}}_{i=1}^{i=N_B}$ is the input and ${\left\{\mathcal{B}\left(t_B^i\right)\right\}}_{i=1}^{i=N_B}$ is the output. For training, only the residual inputs are required ${\{x_F^i,t_F^i\}}_{i=1}^{i=N_F}$, where all data points are picked at random either from the entire computational domain or some subset without repetitions and $N_I+N_B\ll N_F$. The multi-objective optimization goal is to train both supervised ${\mathcal{L}}_{\mathcal{I}},{\mathcal{L}}_{\mathcal{B}}$ and unsupervised ${\mathcal{L}}_{\mathcal{F}}$ losses, where the residual loss ${\mathcal{L}}_{\mathcal{F}}$ plays a role as a regularization term that allows the model to generalize information outside of some local region and prevent overfitting by penalizing the physical constraints and $\alpha$ is the balance between loss importances. It can be seen that with a higher number of initial and boundary data, the physics loss becomes less important $\alpha \to 0$ and vice versa, i.e, with less available data, physics loss becomes more important $\alpha \to 1$.

For the calculations of the partial derivatives for ${\mathcal{L}}_{\mathcal{F}}$, automatic differentiation (AD) is applied, which means that outputs of the neural network $u(x,t)$ are directly used to calculate its derivatives with respect to the inputs $(x,t)$, as $[\partial u/\partial x,\partial u/\partial t,...]$.

In Figure 2, a schematic representation of the PINN is shown with key elements (neural network, automatic differentiation AD, loss) and inputs $(x,t)$ and outputs u and v. AD outputs are differentiated with respect to their inputs in order to construct the ${\mathcal{L}}_{\mathcal{F}}$ loss; then, the outputs are directly used to compute ${\mathcal{L}}_{\mathcal{I}}and{\mathcal{L}}_{\mathcal{B}}$ losses. Then, the final loss $\mathcal{L}$ is computed.

3.2. Discrete-Time Model

The Runge–Kutta (RK4) method for time integration is used as a discrete time model in CFD; it is particularly effective in addressing nonlinear problems. This method is characterized by a fourth-order accuracy in time. The RK4 method is applied to the general form of the case problem:

$$\frac{\partial U}{\partial t}=\mathcal{N}\left(U\right).$$

(17)

The Runge–Kutta method’s four steps are written as:

$$\begin{array}{cc}\hfill U^{\left(1\right)}& =U^n+\frac{\Delta t}{2}{\mathcal{N}}^n\hfill \\ \hfill U^{\left(2\right)}& =U^n+\frac{\Delta t}{2}{\mathcal{N}}^{\left(1\right)}\hfill \\ \hfill U^{\left(3\right)}& =U^n+\Delta t{\mathcal{N}}^{\left(2\right)}\hfill \\ \hfill U^{n+1}& =U^n+\frac{\Delta t}{6}({\mathcal{N}}^n+2{\mathcal{N}}^{\left(1\right)}+2{\mathcal{N}}^{\left(2\right)}+{\mathcal{N}}^{\left(3\right)})\hfill \end{array}$$

(18)

Given the fact that the RK4 scheme only affects time, schemes are needed to approximate the space for convective and viscous terms in order to achieve second-order accuracy in space. Thus, the upwind scheme is used for the flux terms:

$$-u\frac{\partial {\xi }_i^t}{\partial x}=-\frac{u_R{\xi }_R-u_L{\xi }_L}{\Delta x},$$

(19)

where $u_R$ and $u_L$ are the average of the nodal point values and are defined as follows:

$$\begin{array}{c}\hfill u_R=\frac{u_i^t+u_{i+1}^t}{2},u_L=\frac{u_i^t+u_{i-1}^t}{2},\end{array}$$

(20)

The functions ${\xi }_R$ and ${\xi }_L$ depend on the signs of $u_R$ and $u_L$ and are defined as follows:

$$\begin{array}{c}\hfill {\xi }_R=\left\{\begin{array}{c}{\xi }_i^t,u_R\ge 0\hfill \\ {\xi }_{i+1}^t,u_R<0\hfill \end{array}\right.,{\xi }_L=\left\{\begin{array}{c}{\xi }_{i-1}^t,u_L\ge 0\hfill \\ {\xi }_i^t,u_L<0\hfill \end{array}\right.\end{array}$$

(21)

This approach, specifically employed for convective terms, calculates the average values of velocities at the cell boundaries of the spatial grid around each node point. It uses the directions of these velocities to establish at which grid node the values of $\xi$ need to be calculated, forming the differences in the flow. The central difference scheme is used for derivatives not involving non-linearity:

$$\frac{\Delta F{\left(u\right)}_i^t}{\Delta x}=\frac{F{\left(u\right)}_{i+1}^t-F{\left(u\right)}_{i-1}^t}{2\Delta x},$$

(22)

The central difference scheme provides an approximation of spatial derivatives for the viscous part of the equations, where F is the flux function.

4. Computational Experiments

First, the viscid and inviscid Burgers equations are considered as a reference method to verify both approaches. The computational experiments of the PINN and RK4 models will be compared to its analytical solutions graphically and with L norm measures. The following are the results of the case study problem.

The input data for the PINN model, such as the NN setup, the number of initial conditions points $N_I\left[-\right]$, the number of boundary conditions points $N_B\left[-\right]$, the number of interior points $N_F\left[-\right]$, the nonlinear activation function $\sigma (·)$, the learning rate $lr$ for the ADAM optimizer, value $\alpha \left[-\right]$, and the number of iterations epochs, are given in

Table 1. Input data for the PINN model.

Model	NN Setup	${\mathit{N}}_{\mathit{x}}\left[-\right]$	${\mathit{N}}_{\mathit{t}}\left[-\right]$	${\mathit{N}}_{\mathit{I}}\left[-\right]$	${\mathit{N}}_{\mathit{B}}\left[-\right]$	${\mathit{N}}_{\mathit{F}}\left[-\right]$	$\mathit{\sigma }(·)$	lr	$\mathit{\alpha }\left[-\right]$	Epochs
Burgers	$[3\times 64]$	128	128	128	256	16,384	Tanh	${10}^{-3}$	0.25	23,000
Viscous Burgers	$[3\times 64]$	128	128	128	256	16,384	Tanh	${10}^{-3}$	0.1	8000
Euler 1D	$[10\times 20]$	256	256	30	50	13,600	Tanh	$3\times {10}^{-4}$	$5.3\times {10}^{-4}$	40,000

, where $N_x\left[-\right],N_t\left[-\right]$ are the grid sizes for the space and time domain, respectively. In addition to this, input data for the CFD model, which are the grid size $N_x$ and the Courant–Friedrichs–Lewy (CFL) $\left[-\right]$ number, are presented in

Table 2. Input data for the CFD model.

Model	${\mathit{N}}_{\mathit{x}}\left[-\right]$	$\mathit{CFL}\left[-\right]$
Burgers	600	0.25
Viscous Burgers	200	0.25
Euler 1D	1000	0.01

.

For the PINN model, all parameters were found either manually or through a grid search [21,30]. The RK4 model usually does not require any sophisticated parameter search; usually, a larger grid size is used to approximate the derivatives and a smaller CFL number is used to ensure the stability is sufficient.

The precision of the proposed models will be measured with the $L_{\infty }$ norm computed for the time $t_f$ and defined as follows:

$$L_{\infty }=\left|\right|\widehat{u}-u{\left|\right|}_{\infty },$$

(23)

where $\widehat{u}$ and u are analytical and numerical solutions, respectively.

4.1. Verification

Inviscid Burgers’ Equation

Despite the simplicity of Burgers’ equation without a viscous part, it plays a significant role in physics due to its ability to capture essential fluid phenomena such as turbulence and wave propagation. This equation describes the behavior of non-linear waves and shock waves:

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=0,x\in [0,1],t\in [0,1],$$

(24)

with the initial condition:

$$u(x,0)=\frac{1}{2\pi t_s}sin\left(2\pi x\right),x\in [0,1],$$

(25)

where the time at which a shock forms is indicated by the $t_s\left[-\right]$ parameter and set equal to 1. Next, periodic boundary conditions are used:

$$u(0,t)=u(1,t)=0,t\in [0,1],$$

(26)

which lead to the next analytical solution $\widehat{u}(x,t)$:

$$\widehat{u}(x,t)=\frac{1}{2\pi t_s}sin\left[2\pi (x-\widehat{u}t)\right],$$

(27)

which can be solved numerically due to implicity by a simple fixed-point iteration for time $t_n$ and location point $x_i$ with the initial condition as ${\widehat{u}}_i^0$:

$${\widehat{u}}_i^{k+1}=\frac{1}{2\pi t_s}sin\left[2\pi (x_i-{\widehat{u}}_i^k{t}_n)\right],k=0,1,...,$$

(28)

with the following convergence criterion:

$$|{\widehat{u}}_i^{k+1}-{\widehat{u}}_i^k|<{10}^{-15}.$$

This solution can be used continuously to measure the numerical accuracy up to shock formation as well as to test the shock-capturing capability after shock formation.

$f(x,t)$ and $\mathcal{N}\left(U\right)$ are given by:

$$\begin{array}{cc}\hfill & f(x,t):=\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x},\hfill \\ \hfill & \mathcal{N}\left(U\right):=-u\frac{\partial u}{\partial x}.\hfill \end{array}$$

The time step stability criterion for the CFD model is defined through the CFL condition and recalculated for each time iteration n:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & s_{max}^{\left(n\right)}=ma{x}_i\left|u_i^n\right|\hfill \\ \hfill & \Delta t^{\left(n\right)}=\frac{CFL\Delta x}{s_{max}^{\left(n\right)}}\hfill \end{array}\end{array}$$

(29)

Figure 3 illustrates the comparison of the numerical solutions of the inviscid Burgers equation with its analytical solution $\widehat{u}$. The graph highlights the accuracy of both numerical approaches in capturing the complex behavior of this fluid dynamics problem. The $L_{\infty }$ norm is computed for the time $t_f=1$ for both models.

Figure 4 shows a contour plot of the u distribution over the entire computational domain $t\in [0,0.1]$ and $x\in [0,1]$. Here, from time $t=0$, the smooth solution turns into a shock at $t=1.0$, which is also in high accordance with the analytical solution. Thus, these methods capture all shock formations during the entire time period.

Viscous Burgers Equation

The viscosity introduces additional mathematical complexity, allowing for the study of phenomena such as the dissipation of kinetic energy, boundary layer separation, and the transition to turbulence. This transition provides valuable insights into the behavior of fluid flows near solid boundaries and in regions with a high velocity gradient. The viscous Burgers equation is described as:

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}=\nu \frac{{\partial }^{2}u}{\partial x^2},x\in [0,1],t\in [0,1],$$

(30)

where $\nu \left[-\right]$ affects the diffusion term, which characterizes the viscosity in the equation, and is set to 0.01. Next, the initial conditions are prescribed:

$$u(x,0)=sin\left(\pi x\right),x\in [0,1],$$

(31)

followed by the boundary conditions:

$$u(0,t)=u(1,t)=0,t\in [0,1].$$

(32)

The analytical solution $\widehat{u}(x,t)$ is defined as follows:

$$\widehat{u}(x,t)=\frac{2\pi \nu {\sum }_{n=1}^{\infty }{a}_{n}exp(-n^2{\pi }^2\nu t)nsin\left(n\pi x\right)}{a_0+{\sum }_{n=1}^{\infty }{a}_{n}exp(-n^2{\pi }^2\nu t)cos\left(n\pi x\right)},$$

(33)

where $a_0$ and $a_n$ are defined as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & a_0={\int }_0^{1}exp\left[-\frac{1}{2\pi \nu }(1-cos\left(\pi x\right))\right]dx\hfill \\ \hfill & a_n=2{\int }_0^{1}exp\left[-\frac{1}{2\pi \nu }(1-cos\left(\pi x\right))\right]cos\left(n\pi x\right)dx,n=1,2,...\hfill \end{array}\end{array}$$

(34)

$f(x,t)$ and $\mathcal{N}\left(U\right)$ are given by:

$$\begin{array}{cc}\hfill & f(x,t):=\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}-\nu \frac{{\partial }^{2}u}{\partial x^2},\hfill \\ \hfill & \mathcal{N}\left(U\right):=-u\frac{\partial u}{\partial x}+\nu \frac{{\partial }^{2}u}{\partial x^2},\hfill \end{array}$$

The time step stability criterion for time iteration n for the CFD model through the CFL condition is defined as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & s_{max}^{\left(n\right)}=ma{x}_i\left|u_i^n\right|\hfill \\ \hfill & \Delta t^{\left(n\right)}=min\left[\frac{CFL\Delta x}{s_{max}^{\left(n\right)}},\frac{\Delta x^2}{2\nu }\right]\hfill \end{array}\end{array}$$

(35)

Figure 5 demonstrates the comparison of the numerical solutions and the analytical solution of the Burgers equation with viscosity with the $L_{\infty }$ norm measures at time $t_f=1$. The presence of a viscous term in the equation has a distinct impact on the solution. Both methods show a high accuracy in capturing shock waves.

Figure 6 demonstrates the behavior of the fluid flow for the entire computational domain. The flow exhibits the same pattern as in the previous case; from time $t=0$, the smooth solution turns into a shock at $t=1.0$. The shock formations are captured over the entire time period as well.

4.2. Case Study

Since both methods have been verified, a case study is presented for a further analysis of their ability to capture the unique characteristics of gas flows in pipes, which involves contact discontinuity and shock waves. The analytical solutions $(\widehat{\rho },\widehat{u},\widehat{p})$ for the case study are taken from Section 2.4.1 in [35]. In Case 1, the rarefaction shock, the whole tube region/computational domain is divided into five regions and in each part, the solution is found through thermodynamical equations.

$f(x,t)$ and $\mathcal{N}\left(U\right)$ are given by,

$$\begin{array}{cc}\hfill & f(x,t):=\frac{\partial }{\partial t}\left[\begin{array}{c}\rho \\ \rho u\\ E\end{array}\right]+\frac{\partial }{\partial x}\left[\begin{array}{c}\rho u\\ \rho u^2+p\\ u(E+p)\end{array}\right],\hfill \\ \hfill & \mathcal{N}\left(U\right):=-\frac{\partial }{\partial x}\left[\begin{array}{c}\rho u\\ \rho u^2+p\\ u(E+p)\end{array}\right],\hfill \end{array}$$

In this case, the time stability condition for the CFD model for time iteration n is found through the CFL condition as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & s_{max}^{\left(n\right)}=ma{x}_i\left[|u_i^n|+\sqrt{\frac{\gamma p_i^n}{{\rho }_i^n}}\right]\hfill \\ \hfill & \Delta t^{\left(n\right)}=\frac{CFL\Delta x}{s_{max}^{\left(n\right)}}\hfill \end{array}\end{array}$$

(36)

Figure 7 illustrates plots of the density $\rho$, velocity u and pressure p with comparisons to the analytical solutions.

The significance of these plots lies in their ability to demonstrate the accuracy of the numerical methods in capturing the complex behavior of compressible flows. It can be seen that the RK4 model has some oscillations in velocity u and pressure p profiles, while the PINN model is more smooth. The $L_{\infty }$ norms were computed at time $t_f=0.1$. Given that $L_{\infty }$ norms for the PINN and CFD models for the velocity u are 0.1633 and 0.4472 and 0.0261 and 0.1213 for the pressure p, the PINN model significantly outperformed the CFD model by 170 % for the velocity u and by 360 % for the pressure p relatively, which is the main advantage of the PINN model.

In Figure 8, a contour plot for the density $\rho$ over the entire computational domain with the interface, a left nonlinear wave, a right nonlinear wave, rarefaction, and shock formations are shown, where a high pressure is present in the left zone from the discontinuity point $x=0.5$, shown in red, and the right zone with a low pressure is shown in blue. The distribution looks similar at the initial time $t=0$, as shown in the physical domain in Figure 1.

In order to understand the physics behind the estimated results, let us consider the case before any wave has reached the left or right boundary, for instance, at time $t=0.15$. Figure 9 illustrates density, velocity and pressure formation in this case.

Regions I and V show values of the physical quantities that are the same as the initial conditions. Region II denotes the position where the head and tail of the rarefaction wave move leftward. While the solution remains continuous within this region, certain derivatives of the fluid quantities might not exhibit continuity. Region III shows the element position of the fluid that is initially at interface point $x=0.5$ reached at time $t=0.15$. The right end points of Region III and Region IV are called a contact discontinuity. It is observed that the pressure and the normal component of velocity remain continuous across the contact discontinuity. However, the density is not continuous across a contact discontinuity. The end point of Region IV is the location of the shock wave moving to the right. During the shock, all of the quantities will in general be discontinuous.

The pressure results reveal the locations and strengths of the shock waves and rarefaction waves as they propagate through the tube. The density results provide information about the compression and expansion of the fluid as the waves pass through it. The velocity results show how the fluid velocities change as waves propagate through the tube, including its acceleration due to compression and deceleration due to expansion. All the results obtained from computational experiments are sufficient for both models in terms of comparing them with analytical plots and using L norms. The results show that the CFD model outperforms the PINN model, but this was achieved due to the higher number of points for space discretization. The CFD model starts to exhibit oscillations for velocity and pressure, while the PINN model gives a smoother solution. This effect could be explained by the PINN’s nature, where it integrates physical laws into the numerical solution. Thus, it creates smooth and continuous solutions. Although the CFD model gives better computational and optimal results for the case of one-dimensional problems, as the dimensionality of the problem and the number of equations increase, the CFD model will require a large grid size and a long computational process and the obtained solutions may exhibit oscillations. Meanwhile, the PINN model can generalize the information using a subset of the entire computational domain in a relatively small number of iterations compared to the CFD model, resulting in a smoother and more optimal solution.

5. Discussion

In the realm of numerical techniques for solving differential equations, both classical CFD and PINN methods are unique and offer distinct advantages. Understanding the comparative strengths and limitations of these methods is crucial for selecting the most appropriate approach for a given problem.

Accuracy and Precision:

Classical CFD methods are renowned for their ability to achieve a high accuracy, especially when employing fine grid resolutions. These methods excel in problems where the underlying physics is well understood and can be accurately represented by explicit equations. However, CFD methods may encounter challenges near discontinuities or sharp gradients, leading to accuracy degradations in such regions.

Conversely, PINNs adopt a data-driven approach, leveraging neural networks to approximate solutions based on available data. While PINNs may sacrifice some accuracy compared to FDMs, particularly in problems with well-defined physics, they offer advantages in handling complex or data-driven scenarios. PINNs can provide smoother and more continuous solutions, making them suitable for problems with irregular geometries or unstructured data.

Computational Efficiency:

Classical CFD methods typically involve explicit grid generation and solving equations at each grid point, which can be computationally intensive for large-scale problems. Conversely, PINNs offer computational efficiency once trained, as they can rapidly generate solutions without the need for explicit grid structures. This efficiency is particularly advantageous for problems with complex geometries or time-dependent constraints.

Flexibility and Adaptability:

Classical CFD methods are constrained by grids shapes and may struggle with problems involving complex or irregular geometries. In contrast, PINNs offer greater flexibility and adaptability, as they can handle unstructured data and irregular domains more effectively. PINNs have the potential to generalize well to new conditions or scenarios not explicitly included in the training data, enhancing their applicability to a wide range of problems.

In conclusion, the choice between classical CFD and PINN methods depends on the specific requirements of the problem, including the nature of the differential equations, the availability of data, the complexity of the geometry, and the computational resources. CFD methods excel in problems with well-understood physics and structured domains, offering a high accuracy with a fine grid resolution. On the other hand, PINNs provide flexibility, generalization, and efficiency advantages, especially for complex or data-driven problems with irregular geometries. By understanding the comparative strengths and limitations of these methods, researchers and engineers can make informed decisions to select the most suitable approach for their specific application.

6. Conclusions

This article investigates gas flow in a pipeline to capture the propagation speed of the rarefaction wave, the contact discontinuity and the shock discontinuity. First, the physical model is defined to provide gas characteristics. Initial and boundary conditions were specified according to the case problem. Then, two numerical models were proposed: a continuous symmetry-enhanced PINN model and a discrete CFD model (Runge–Kutta 4). The methods were verified through the one-dimensional inviscid and viscous Burgers equations. This was followed by case study problem estimation, which is described by a system of one-dimensional Euler equations. The results show a sufficient accuracy in the graphical comparison with the analytical solution and in the L norm value. According to the results of numerical studies, it can be seen that for relatively simple problems, the classical CFD method is still superior to PINNs, but with an increasing complexity of the problem, CFD suffers from increasing grid sizes and as a result may exhibit oscillations. In contrast, the PINN presents a smoother solution, which in turn leads to the expectation of a better result with an increasing dimensionality of the problem and an increasing number of input equations. Representing the main numerical simulation result, the PINN model was able to outperform the CFD model in terms of $L_{\infty }$ norms, improving the velocity u by 170% and the pressure p by 360% in a relative comparison. Moreover, graphical results that are indistinguishable from the analytical solution for both models were obtained; in the case of the Burgers equations, plots at time 0.3, 0.7, and 1 and also contour plots for the whole computational domain were obtained, and for the Euler equations, plots at time 0.08, 0.09, and 0.1 were obtained, along with contour plots for the density over the whole computational domain. In conclusion, this study highlights the effectiveness of both Physics-Informed Neural Networks and classical CFD methods in modeling gas flows in pipelines. By accurately capturing the dynamic flow characteristics such as rarefaction waves and shock discontinuities, these numerical techniques provide valuable insights into the behavior of gas transportation systems. The comparison between PINN and CFD approaches underscores the importance of selecting the most appropriate method based on the specific requirements of the problem, ensuring reliable and efficient solutions for modeling gas flows in pipelines.

Future research could be aimed at improving both models: For the continuous PINN model, a way to find the most appropriate hyperparameters using, for instance, Bayesian hyperparameter optimization, should be found. For the discrete CFD model, the approximation of derivatives using the finite difference method should be replaced by more accurate methods with a higher order of approximation over space and a flexibility in complex domains, such as the finite volume method. In addition to this, the physical problem case study (Euler’s equations) should be extended to two-dimensional forms and even further to two-dimensional Navier–Stokes equations.