Modeling and Simulation of Real Gas Flow in a Pipeline

In this paper, a mathematical model that describes the flow of gas in a pipe is formulated. The model is simplified by making some assumptions. It is considered that the natural gas flowing in a long horizontal pipe, no heat source occurs inside the volume, transfer of heat due to heat conduction is dominated by heat exchange with the surrounding. The flow equations are coupled with equation of state. Different types of equations of state, ranging from the simple Ideal gas law to the more complex equation of state Benedict Webb Rubin Starling (BWRS), are considered. The flow equations are solved numerically using the Godunov scheme with Roe solver. Some numerical results are also presented.


Introduction
The purpose of this paper is to describe the flow of natural gas in a pipeline by employing the full set of differential equations along with different types of equations of states(EOS), ranging from the simple Ideal gas law to the more complex equation of state, Benedict Webb Rubin Starling (BWRS).The flow equations are derived from the physical principles of conservation of mass, momentum, and energy.More detailed discussion of conservation laws can be found in [1]- [4].The natural gas is inviscid and compressible.The gas flows in along a horizontal pipe, and then can be considered as one-dimensional flow.It is assumed no heat source occurs inside the pipe and transfer of heat due to the heat conduction is much less than the heat exchange with the surrounding.
In this paper, the results obtained by solving the flow equations along with different types of EOS are compared [5].The ideal gas equation works reasonably well over limited temperature and pressure ranges for many substances.However, pipelines commonly operate outside these ranges and may move substances that are not ideal under any conditions.The more complicated EOS will approximate the real gas behavior for a wide range of pressure and temperature conditions.The Godunov scheme with Roe solver [3] is used to solve the Euler equations numerically.The Godunov scheme for conservation laws is known for its shock-capturing capability.
The rest of the article is organized as follows.In Section (2) we review the set of partial differential equations which describe the flow of gas in a pipe.Several equations of states are discussed in this section.In Section (3) a thermodynamical relationships among the physical quantities are presented.One can refer [6] for more thermodynamical relationships.Section (4) contains the discussion of the numerical method used to solve the flow equations together with different types equation of states.Some numerical results are given in this section.Conclusions are deferred to Section (5).

Governing Equations of Real Gas Flow in a Pipe
Let us consider a gas occupying a sub domain 0 ∈ Ω describes the posi- tion of the particle x at time t.Then at time t the gas occupies the domain ∈ Ω .The velocity of the gas at position x and time t is given by ( ) ( )

Transport Theorem
Let [ ] : 0, →  be some physical quantity transported by the fluid.The total amount ( ) F t of the quantity f contained in t Ω a time t is given by ( ) ( ) The rate of change of ( ) F t is given by: ( Then we get the transport theorem: . The transport theorem is useful in the derivation of the governing equations.

Conservation of Mass (The Continuity Equation)
The total mass m in a volume t Ω is given by ( ) Since the above integral holds true for arbitrary region t Ω

Conservation of Momentum (Equation of the Motion)
The total momentum M of particles contained in t Ω is given by According to Newton's second law: The rate of change of momentum equals the action of all the forces F applied on t Ω We have two types of forces acting on t Ω : 1) Volume forces v f , for example gravitation, which is given by where g is the gravitational acceleration.
and n is the outer normal.The total force By applying divergence theorem, the second term on the right side of the above equation can be transformed to integral over the domain t Ω and then we get: For Newtonian fluid, the stress tensor depends linearly on the deformation velocity, u ∇ , i.e.

(
) where τ is the viscous part of σ , p is pressure, I is the identity matrix, λ and µ are friction coefficients, and D is the strain tensor given by ( ) ( ) For inviscid fluid, friction is neglected and then pI σ = − Therefore, the equation of motion for inviscid fluid becomes ( ) ( )

Conservation of Energy
Conservation of energy accounts for effects of temperature variations on the flow or the transfer of heat with in the flow.The 1 st Law of Thermodynamics states that: The total energy of a system and its surroundings remains constant.
Let  be the total energy of the fluid in t Ω and Q be the amount of heat transfered to t Ω .The rate of change of the total energy of the fluid occupying t Ω is the sum of powers of the volume force acting on the volume t Ω , powers of the surface force acting on the surface where q is the density of heat sources (per unit mass), and q is the heat flux (transfer of heat by conduction).The transfer of heat by conduction is given by Fourier's law: where T is the absolute temperature and 0 κ ≥ is the coefficient of thermal conductivity of the fluid.
q is the density of heat transfered from the surrounding and is given by: ( ) where L k is the total heat transfer coefficient and out T is the temperature of the surrounding.
Then the energy equation for inviscid gas flow becomes: x t E x t x x t g x t u x t x pI u x t s t x t q x t x q x t s q x t x By applying the transport and divergence theorems to the above equation we obtain the following equation: surface area volume There fore, from the equations (1), ( 2), (3) we get the following system of equations. (

Simplifications
In practice the form of mathematical model varies with the assumptions made as regards of operation conditions of the pipeline.Simplified models are obtained by neglecting some terms in the basic equations.In our case, we consider natural gas (Methane) flowing in a long horizontal pipeline.Hence we can consider the flow as a one dimensional flow.By assuming the pipe is horizontal, we can neglect the contribution of the gravitational force.
Assume also no heat source occurs inside the volume.For a cylindrical pipe, surface area 4 volume where D is the diameter of the pipe.By applying the assumptions we made, ( 4) is reduced to Furthermore, Methane gas has the following properties.The specific heat capacity [ ] Prandtl number (Pr), defined as Another dimensionless constant we can use to simplify our system of equations is the Nusselt number (Nu).The Nusselt number is defined as , where D is a characteristic width of a flow, for example the diameter of the pipe.The Nusselt number compares convection heat transfer to fluid conduction heat transfer.For Methane gas flowing through an insulated pipe of diameter 0.5 m buried underground, the value of Nu is approximately 10.If the pipe is exposed to air, it will be around 300.Therefore, the term included in the energy equation due to heat conduction can be neglected in favor of the term due to heat exchange with the surrounding ( ) . Incorporating these assumptions to Equation (5) yields: where ( ) is an equation of state used to complete the system of conservation laws.In the next chapter we will solve Equation (6) with different equation of state numerically.

Equation of State (EOS)
An equation of state is a relationship between state variables, such that specification of two state variables permits the calculation of the other state variables.For an ideal gas, the equation of state is the ideal gas law.More complicated EOS have been formulated by several workers to try to model the behavior of real gases over a range of pressures and temperatures.This includes Van der Waals (VDW), Sovae Redlich Kwong (SRK), Peng Robinson (PR), and Benedict Webb Rubin Starling (BWRS).

Ideal Gas law
The ideal gas law is given by where p is the pressure, ρ is the density, R is the gas constant, and T is the absolute temperature.
The ideal gas law is derived based on two assumptions:  The gas molecules occupy a negligible fraction of the total volume of the gas. The force of attraction between gas molecules is zero.
The ideal gas equation works reasonably well over limited temperature and pressure ranges for many substances.However, pipelines commonly operate outside these ranges and may move substances that are not ideal under any conditions.Hence, we need to look for equation of state with wider validity.

Van der Waals (VDW) EOS
It was observed that the ideal gas law didn't quite work for higher pressures and temperatures.The first assumption works at low pressures.But this assumption is not valid as the gas is compressed.Imagine for the moment that the molecules in a gas were all clustered in one corner of a cylinder, as shown in the figure below.At normal pressures, the volume occupied by these particles is a negligibly small fraction of the total volume of the gas.But at high pressures (when the gas is compressed), this is no longer true.As a result, real gases are not as compressible at high pressures as an ideal gas.The volume of real gas is therefore larger than expected from the ideal gas equation at high pressures.Van der Waals proposed that we correct for the fact that the volume of real gas is too large at high pressures by subtracting a term from the volume of the real gas before we substitute it in to the ideal gas equation.He therefore introduced a constant b in to the ideal gas equation that was equal to the volume actually occupied by the gas particles.When the pressure is small, and the volume is reasonably large, the subtracted term is too small to make any difference in the calculation.But at high pressures, when the volume of the gas is small, the subtracted term corrects for the fact that the volume of a real gas is larger than expected from the ideal gas equation.
The assumption that there is no force of attraction between the gas particles cannot be true.If it was, gases would never condense to form liquids.In reality, there is a small force of attraction between gas molecules that tends to hold the molecules together.This force of attraction has two consequences: (1) gases condense to form liquids at low temperatures and (2) the pressure of a real gas is sometimes smaller than expected for an ideal gas.To correct for the fact that the pressure of a real gas is smaller than expected from the ideal gas equation, Van der Waals added a term to the pressure in the ideal gas equation.This term contains a second constant a.The complete Van der Waals equation is written as follows: Or in terms of molar volume ( ) P critical pressure, and c T critical temperature Note that the values of the constants a and b differ from gas to gas.Even though, VDW EOS is better than Ideal gas law, still it is inadequate to describe real gas behavior.
We will consider three widely used equations of state that do work reasonably well near the dew point: Sovae-Redlich-Kwong (SRK), Peng-Robinson (PR), and Benedict-Webb-Rubin-Starling (BWRS).In addition to covering a wide range of conditions, these equations also can be expressed in generalized forms with mixing rules that permit the calculation of the coefficients for different compositions.
SRK and PR, along with VDW are called cubic equation of state, because expansion of the equations into a polynomial results in the highest order terms in density (or specific volume) being cubic.BWRS adds fifth and sixth power and exponential density terms.The cubic equation are all of the form 2 2 1 1

The Sovae-Redlich-Kwong (SRK) EOS
The SRK EOS of state is given by where is the reduced temperature.

The Peng-Robinson (PR) EOS
The PR EOS is defined as where

Benedict-Webb-Rubin-Starling (BWRS) EOS
Probably because of its ability to cover both liquids and gases and the availability of coefficients and mixing rules for many hydrocarbons in one place, BWRS is the most widely used equation of state for simulation of pipelines with high density hydrocarbons, or with condensation.
Simplicity is not among the good qualities of the BWRS equation of state.The form of the equation is: where the eleven coefficients , , , , , , , , , A B C D E a b c d α , and γ are determined from , , , and ω of the gas of interest and the universal constants i A and i B as follows.
( ) ( ) ∑∑ where i x is the mole fraction of the pure component i of the mixture, and ij k are the binary interaction coef- ficients.

The Universal Gas Law
The universal gas law is p Z RT ρ = where Z is called the compressibility factor (Real gas factor).It is a measure of how far the gas is from ideality.At atmospheric conditions, the value of Z is typically around 0.99.Under pipeline conditions, the value is typically around 0.9.A good equation of state can be selected by its ability to approximate the compressibility factor at critical conditions c Z .For example the experimental value of c Z for Methane is 0.288.But its approximate value by VDW is 0.3025, by SRK is 0.2904, by PR it is 0.2894, and by BWRS it is 0.2890.

Thermodynamical Relations
In this section we will briefly discuss thermodynamical relations that exist among different physical quantities.First law of thermodynamics states that The specific total enthalpy is defined as h e pv = + which implies Derivative relationships: Assume ( ) . Comparing the coefficients of this equation to that of Equation ( 13) we get Similarly, assuming ( ) And comparing the coefficient of this equation with that of Equation ( 14) we get

Reciprocal relations involving internal energy e and entropy s:
Consider the internal energy and entropy to be a function of temperature and specific volume, i.e, ( ) The coefficient of dT, in the first equation, is by definition the heat capacity at constant volume, v c .Substi- tute these two equations in (13) to get , Differentiating the first equation of (18) with respect to T and the second with respect to v gives us Substituting (19) in the first equation of (18) yields One useful form involving internal energy is obtained by substituting v c for the coefficient of dT in (20) for the coefficient of dv in the first equation of (17).
Reciprocal relations involving enthalpy h Assume ( ) The coefficient of dT is by definition the heat capacity at constant pressure, p c .In a similar procedure as in the internal energy and entropy case, above we get the following relationships. , By double differentiating we do get

Heat capacities
By equating the difference of ( 13) and ( 14) to the difference of ( 21) and (25) we get Dividing by dT and holding p constant gives ( )

Numerical Methods: Godunov Scheme with Roe Solver
In this section we will consider a numerical scheme to solve homogeneous Euler equation with initial condition by employing different EOS.The Euler equation in vector form: x a = and Let us define ( ) . Then we solve exactly the local Riemann problem for , , for Let us denote the solution by ( ) Since v is an exact solution on With the numerical flux ( ) This scheme is called Godunov scheme.Solving a Riemann problem exactly is not always an easy task.Then we may need to consider an approximate solution of the Riemann problem.

Riemann Problem for a Linear System
Suppose we have a linear system 0 λ λ λ < < are the eigenvalues and 1 2 3 , , r r r are the corresponding eigenvectors.Define , 1, 2, Then the solution of the Riemann problem is given by ( ) A variety of approximate Riemann solvers have been proposed that can be applied more easily than the exact Riemann solver.One of the most popular Riemann solvers currently in use is due to Roe.

Godunov scheme with Roe approximation.
The idea is to replace the non-linear Riemann problem solved at each interface by an approximate one.
( ) where l U and r U are the left and right values and ( ) U U U → Conservation form of the Roe scheme.The Roe scheme can be written in conservation form as ( ) ( )

∑
where i λ and i r are the eigenvalues and eigenvectors of ( ) , A u w and To solve our problem with the Roe scheme, we need to calculate the eigenvalues and their eigenvectors of the Jacobian matrix ( ) DF U which are needed to compute the Roe flux.But for complex EOS the determination of these eigenvectors may not be simple.One way of determining the eigenvectors of this Jacobian is by expressing the Euler equation in terms of primitive variables ( ) We choose the temperature T as one of primitive variables than the pressure p, because in most equation of state p is expressed in terms of T.
Let 0 V BV + = be the Euler equation in terms of the primitive variables V and ( ) 0 t x U F U + = be in conservative variables.The approximate linear system is ( ) , so that we can apply the Roe scheme easily.
Continuity equation: , the momentum equation in terms of the primitive variables is 0 Now, using Then the Euler equation in primitive variables is written as Or in vector form 0 Eigenvalues and eigenvectors of the coefficient matrix B of (31) are computed as follows.
where the local speed of sound c is given by The matrix of the corresponding eigenvectors is: To compute the eigenvectors of the Jacobian ( ) DF U we need to compute the matrix The matrix R of eigenvectors of ( ) DF U is given by: ( ) Since the total specific enthalpy h is given by we can write the eigenvectors in terms of h as ( )

Solving Euler Equation Using the Van der Waals (VDW) EOS
Here we solve one dimensional Euler equation with VDW EOS.Consider again the euler equation ( 28 Again using (21), the change of internal energy is given by: Integrating the above differential equation gives the internal energy v e c T aρ = − .The total energy  is given by: , and Hence the momentum equation is reduced to Energy Equation ( ) , and and the coefficient of x p is 0. Then (32) reduces to ( ) The Euler equation is written as where ( ) Eigenvalues and eigenvectors of the coefficient matrix B of (34) are computed as follows.
21 23 32 0 0 , and where the local speed of sound c is defined as The matrix of the corresponding eigenvectors is: To compute the eigenvectors of the Jacobian ( ) DF U we need to compute the matrix ( ) The matrix R of eigenvectors of ( )   is the accentric factor R is gas constant, c P critical pressure, c T critical temperature, and r T is the reduced temperature.

R u c u u c h uc r h uc
The internal energy is given by: where ( ) After integrating the differential equation of the internal energy, we get ( ) The total energy  is given by: ( ) ) momentum equation is written as Energy Equation: ( ) The Euler equation is written as ) Eigenvalues and eigenvectors of the coefficient matrix A of Equation (37) are given as follows.The matrix of the corresponding eigenvectors is: To compute the eigenvectors of the Jacobian ( ) DF U we need to compute the matrix where ( ) The matrix R of eigenvectors of ( ) DF U is given by: Since the specific enthalpy h is given by  , w is the accentric factor R is gas constant, c P critical pressure, c T critical temperature, and r T is reduced temperature.The internal energy is given by: The total energy  is given by: , and The momentum equation is written as Energy Equation ( )  The Euler equation is written as The matrix of the corresponding eigenvectors is: To compute the eigenvectors of the Jacobian ( ) Since the specific enthalpy h is given by ( ) ( ) The internal energy is given by: The total energy  is given by: ( ) The momentum equation is written as and the coefficient of     The conservative variables ( ) , , u ρ ρ  are determined by the scheme.The velocity is obtained from ρ and u ρ .But to determine the value of the temperature T we use an iteration method (especially for the cases of complex EOS).Then the pressure P is computed from the EOS

Numerical Results
In this section we present some numerical results.We consider a tube of length 1, filled by Methane gas, the initial discontinuity is located at 0 0.5 x = .In our simulation the following initial data is used. .In Figure 1, we have plotted the density, pressure, velocity, temperature, and the real gas compressibility factor computed by using each of EOS we discussed.
Figure 2 depicts results of (6), i.e, the Euler equation with the source term included, obtained by applying PR, and BWRS EOS.

Conclusion
The model that describes the flow of gas in a pipe is presented.Simplifications to the equations are made using appropriate assumptions.Several Equations of states that close the system of equations are examined and the results obtained for each equation of state are compared.

(
By transport theorem) relative strength of viscosity (the diffusion of momentum) to that of heat.It is entirely a property of the fluid not the flow.In our case the value of Pr is about 0.7, this enables us to regard the flow as inviscid flow.For gas flow typical values of Pr are between 0.7 and 1.

w
is the accentric factor which is a measure of the gas molecules deviation from the spherical symmetry, R is gas constant, c P critical pressure, c T critical temperature, and r BWRS can be adapted for mixtures by the rules:

=
One of the methods to solve a 1D nonlinear hyperbolic systems is the Godunov scheme 4.1.Godunov Scheme Suppose we have subdivided our domain [ ] , a b in to N subintervals with 1 of the local Riemann problems are used to define the global solution v as

.
The main task in the Roe scheme is the determination of the matrix of linearization A. Now let us consider our equation (28) together with an equation of state of the form this non-linear system with an approximate linear system as follows: enthalpy.These averages are called the Roe mean values.( )

4 . 3 .
Solving Euler Equation Using the Ideal Gas LawIn this section we solve one dimensional Euler equation with Ideal gas EOS.Consider the Euler equation (28) with the ideal gas law p RT ρ = .Using (21), the change of internal energy is given by d d Now let us express (28) in terms of the primitive variables

4 . 6 .
write the eigenvectors in terms of h as Solving Euler Equation Using the Peng-Robinson (PR) EOS Let us consider (28) with PR EOS.
Eigenvalues and eigenvectors of the coefficient matrix B of Equation (43) are computed as follows.

Figure 1 .
Figure 1.The results obtained by solving the homogeneous Euler equation by employing the ideal gas law and the other four equation of states.

Figure 2 .
Figure 2. The results obtained by solving the Euler equation (including the source term) by employing the PR and BWRS EOS.

.
The last equation is a system of simultaneous algebraic equations for the variables i α .

3 l
d

4.8. Application of the Roe solver
Now to apply the Roe scheme on (28), on each cell [ ]