An approximate well‐balanced upgrade of Godunov‐type schemes for the isothermal Euler equations and the drift flux model with laminar friction and gravitation

In this article, approximate well‐balanced (WB) finite‐volume schemes are developed for the isothermal Euler equations and the drift flux model (DFM), widely used for the simulation of single‐ and two‐phase flows. The proposed schemes, which are extensions of classical schemes, effectively enforce the WB property to obtain a higher accuracy compared with classical schemes for both the isothermal Euler equations and the DFM in case of nonzero flow in the presences of both laminar friction and gravitation. The approximate WB property also holds for the case of zero flow for the isothermal Euler equations. This is achieved by defining a relevant average of the source terms which exploits the steady‐state solution of the system of equations. The new extended schemes reduce to the original classical scheme in the absence of source terms in the system of equations. The superiority of the proposed WB schemes to classical schemes, in terms of accuracy and computational effort, is illustrated by means of numerical test cases with smooth steady‐state solutions. Furthermore, the new schemes are shown numerically to be approximately first‐order accurate.

decision-making process; for instance, pipelines are usually designed to operate in the steady condition. 3 Moreover, mathematical models are used for leak detection in pipelines by comparing the measurements with the numerical steady-state solution. 4 This highly relies on the accuracy of the numerical solution. In addition, transient multiphase flows commonly occur in pipelines when changes in operational conditions, such as inlet and outlet flow rates, and set-point pressures, are induced. These changes are usually exerted to reach a new steady condition in the system. All these points dictate that a reliable simulator should predict not only transients accurately but also the steady-state solution.
Many numerical methods have been used for solving the equations governing the physics of a phenomenon. When there is no effect from external sources, these methods are often highly accurate in predicting steady-state behaviors of the system. However, many realistic industrial systems, such as, for example, managed pressure drilling systems, are inevitably affected by external sources such as friction and gravitation. 5 It has been observed that classical finite-volume schemes do not preserve the analytical (or the trustworthy numerical) smooth (i.e., continuous and differentiable) steady-state solution of systems in the presence of such source terms. [6][7][8] To resolve this issue, much effort has been put into deriving schemes capable of preserving the analytical steady-state solution. Such schemes are called well-balanced (WB) schemes. 9,10 The isothermal Euler equations (henceforth called Euler equations) 11 and the drift flux model (DFM) 12 have been widely used for modeling single and two-phase flows in pipelines. The accuracy of Euler equations and the DFM for single-phase and two-phase flows in drilling scenarios has been verified by comparing it to real-life field data 5,13 . Thus, these models are trustworthy for simulation of single-phase and two-phase flows in pipelines and in drilling. Developing WB schemes for Euler equations and the DFM is necessary to approximate the correct steady behavior of the flow in pipelines. To capture the analytical steady-state solution of Euler equations, a few studies have been carried out for special cases. [14][15][16][17] The developed schemes in these studies are all WB only with respect to gravity, and not friction. Moreover, the proposed solution in References 14-16 relies on the analytical steady-state solution, which is, in general, not available or is computationally expensive to be computed at each time step. Furthermore, schemes in References 14-16 are only WB in case of zero flow (stationary steady-state solutions), which restricts the applicability of these schemes. Stationary steady-state solution of Euler equations is, however, much easier to capture compared with moving (nonzero flow) steady-state solutions as the mass conservation law is automatically satisfied in the stationary situation. Moreover, the analytical stationary steady-state solution replicates a trivial hydrostatic solution, which might not be the most important scenario in drilling for oil and gas and also the transportation of liquid and gas through pipelines. The solution in Reference 17, although applicable for moving steady-states, is based on neglecting the diffusive part of Rusanov scheme 18 when the system is in its steady state (Rusanov scheme then converts to a centered scheme when their algorithm detects that the system has reached steady condition). Nonetheless, the diffusive part of this scheme is essential for the stability of the solution during transients. Moreover, the mechanism that detects the solution is now steady and the diffusive part should be neglected is not well explained. Therefore, in this study, a different method is proposed to solve Euler equations in a WB manner in a general scenario (zero and nonzero flow) in the presence of both laminar friction and gravitation with an accuracy much higher than classical schemes; however, our method does not capture steady-state solution exactly and therefore we recover the actual steady-state solution approximately. This leads to an approximate WB scheme. Moreover, we also propose a scheme for the case of two-phase flows governed by the DFM. To the best of authors' knowledge, no study has been performed on capturing the smooth steady-state solution of multiphase flows, especially the DFM.
Contrary to Euler equations and the DFM, the shallow water equations have attracted many researchers for developing a WB scheme. 9,10,16,19,20 These efforts were pioneered by developing a WB scheme for a lake at rest, 16,20 further extended to nonzero velocity with topography 10 and friction source terms. 9 In case of nonzero velocity, the schemes become approximately WB as many assumptions have to be made. Moreover, the convergence rate of the scheme also reduces due to these simplifying assumptions. 9,10 Fortunately, Euler equations, 17 the DFM, 21 and the shallow water equations 10 share many common features; for example, all three systems are hyperbolic (for the DFM, it is hyperbolic over a wide region of the state variables 22 ) and these systems are typically discretized by finite-volume techniques. Therefore, we propose to generalize the ideas proposed for the shallow water equations and extend these techniques to Euler equations and the DFM. Compared with the studies on shallow water equations, apart from the application for other classes of partial differential equations (PDEs), this study differs in another crucial aspect. In References 9,10, the intermediate values in the approximate Riemann solver are found to force the system to be approximately WB. Here, we modify an existing scheme to upgrade it approximately WB. All in all, this research aims to develop a scheme-and model-dependent framework for increasing the accuracy of numerically obtained steady-state solutions of Euler equations and the DFM. Nonetheless, this approach can be extended straightforwardly to other numerical schemes and also to other systems of hyperbolic PDEs. This article is organized as follows. In Section 2, Euler equations and the DFM together with their steady-state solutions are introduced. In Section 3, a methodology to upgrade a finite-volume scheme to become WB is discussed. Moreover, the application of this methodology to an advection equation, Euler equations and DFM is elaborated (the exact WB property of the scheme is proved for the advection equation). In Sections 4 and 5, the proposed schemes are supported by illustrative numerical test cases for Euler equations and the DFM, respectively. In Section 6, the error convergence of the schemes is studied. Finally, Section 7 concludes the article.

SINGLE-AND TWO-PHASE FLOW MODELS
In this section, Euler equations and the DFM, together with the corresponding steady-state solutions, are introduced.

Isothermal Euler equations
Single-phase flow inside a pipe can be modeled by the isothermal Euler equations. 23 This system of equations is as follows: where w, q(w, x) and f (w) are, respectively, the vector of conservative variables, the source terms and the mathematical flux function defined as with (t, x), u(t, x), and p(t, x) denoting density, velocity, and pressure of the fluid, respectively. The temporal and spatial variables are denoted by t and x while T and L are the final time of the simulation and the length of the spatial domain (i.e., the length of the pipe). Moreover, F(w, x) is the laminar friction and G(w, x) is the gravitational source term. In this article, we consider where , d, g, and are, respectively, the viscosity of the fluid, the hydraulic diameter of the pipe, the gravitational acceleration, and the pipe inclination with respect to the horizontal plane. System (1) is completed by the equation of state (EOS) of the fluid as follows: where 0 and p 0 are the reference density and pressure for defining the EOS of the fluid and c is the constant speed of sound in the medium occupied by the fluid. As the dynamics for single-phase flow are now uniquely described, we can proceed to find a numerical steady-state solution of system (1).

Steady-state solution of the isothermal Euler equations
The equations describing the steady-state solution are obtained by using w/ t = 0 in (1) and substituting p from (4) into (1), which read as dm where m = u is the momentum. For subsonic flows (which is the common case for drilling applications and transport of gas and liquid), the two boundary conditions to be specified to solve system (5) lead to a two-point boundary value problem (BVP), which is hard to be solved analytically, especially because both friction and gravity are present in q. However, this BVP can be solved numerically by the bvp4c solver of MATLAB. This solver approximates the solution to (5) in an iterative way while considering the boundary values at both ends of the computational domain. 24 This numerical solution later serves as a reference solution to evaluate the accuracy of the proposed scheme in predicting the steady-state solution. However, due to the nature of the BVP, the steady-state solution is computationally expensive and thus the approach from References 14-16 is not applicable here to develop a WB solver for Euler equations. Moreover, since friction and gravity are present simultaneously and we aim to maintain all characteristics of the scheme (not neglecting the diffusive part), the methodology from Reference 17 cannot be employed either. Thus, we try to find a more general approach to develop a WB numerical solver in Section 3. Next, we analyze the steady-state solution in case of only laminar friction and only gravitational force, which later supports our arguments in Section 3.2.2. (1) with (4) for the case of only friction source term (3a) has (a unique or two) steady-state solutions for a reference point

Lemma 1. The isothermal Euler Equations
Proof. In the case of only laminar friction as the source term, (5) reads as Since the momentum is constant along the spatial domain in the steady state, we consider m = m 0 and rewrite (7b) as follows: Consider x 0 ∈ R as a reference point with (x 0 ) =̂0 and integrating (8) over (x 0 , x) and denoting (x) = , we have: The minimum of occurs at Function is strictly increasing on ∈ ( c , ∞) and decreasing on ∈ (0, c ). Therefore, function admits a unique minimum for = c . As → ∞ when → 0 and → ∞, to have any steady solution, ( c ) ≤ 0 must hold. This leads to the condition (6).
Proof. Similar to Lemma 1, we have Consider x 0 ∈ R as a reference point with (x 0 ) =̂0 and integrating above over (x 0 , x) and showing (x) = , we have: The minimum of occurs at Function is strictly increasing on ∈ ( c , ∞) and decreasing on ∈ (0, c ). Therefore, function admits a unique minimum for = c . As → ∞ when → 0 and → ∞, to have any steady solution, ( c ) ≤ 0 must hold. This gives the condition (11). ▪

The DFM
The DFM is described by the following system of PDEs, 25 which describes a two-phase flow inside a pipe: where the subscripts l and g denote the liquid and gaseous phase, respectively, and l (t, x) and g (t, x) denote the volume fraction of each phase. Frictional and gravitational terms are given by with mix = l l + g g the mixture viscosity, u mix = l u l + g u g the mixture velocity ( l u l and g u g are the superficial velocities of each phase), and mix = l l + g g the mixture density of the gas and liquid. The DFM is completed by closure relations, as listed below: 21,26 where K and S are two constant parameters depending on the flow regime. 27,28 Moreover, c l and c g are the constant speed of sound in the liquid and gas medium, respectively. Now, the steady-state solution of system (15) can be computed numerically.

Steady-state solution of the DFM
The steady-state solution of the DFM can be obtained by solving the following system of equations, which is obtained by setting w/ t = 0 in (15): where m i = i i u i , i ∈ {l, g}, is the momentum of phase i. By embedding the closure relations (18) in (19), this system, which also leads to a two-point BVP due to the boundary conditions for subsonic flow, can be solved numerically. This numerical steady-state solution will be used as the reference solution to assess the WB property of the numerical solver to be proposed. Similar to the discussion for Euler equations in Section 2.1.1, the approaches proposed in References 14 and 17 cannot be employed for solving system (15) and new approaches should be developed. In the following section, the novel WB schemes are introduced.
Remark 1. Notably, other friction functions rather than laminar friction can also be studied; however, more complicated frictional source terms complicates the analysis presented in this article, if not impossible.
Remark 2. The analysis in this article holds true if the parameters present in the source terms such as , hydraulic diameter and viscosity vary smoothly along the spatial domain. In this case, the steady-state solution of Euler equations and the DFM will also be smooth, which will be used in the proofs presented in Sections 3.2.2 and 3.2.3. For the sake of simplicity, these parameters are assumed to be constant in this study.

WELL-BALANCED FINITE-VOLUME SCHEME
Classical cell-centered finite-volume schemes are reliable for solving systems in the absence of source terms. 11,18,21 When source terms appear in the governing equations, these numerical methods are no longer WB; that is, the steady-state solutions predicted by the finite-volume solvers differ from the analytical steady-state solution. As a consequence, the solution of the finite-volume method deviates from the analytical steady-state solution even when initialized on the steady-state solution. This deviation can be made smaller by increasing the number of grid cells to a computationally unfeasible large number (much more than the required spatial resolution). Here, we propose a method that can achieve a significantly higher accuracy in predicting the steady-state solution with relatively low number of grid cells. As the method we are suggesting is scheme-dependent, in this section, first Rusanov scheme 18 is introduced as a reference classical scheme, and then a modification of the scheme is proposed, which is able to compute an accurate approximation of the analytical steady-state solution. Reasons for choosing Rusanov scheme are its simple formulation compared with other numerical schemes and the independence of its diffusivity properties on the Courant-Friedrichs-Lewy (CFL) number. 18 These features yield less complicated nonlinearities. However, the methodology introduced in this article can be applied straightforwardly to other Godunov-type schemes as well.

Rusanov scheme
Let Δt and Δx refer to the temporal and spatial discretization intervals over time and space, respectively. The spatial discretization consists of cells spatially located between two interfaces (x i − 1/2 , x i + 1/2 ) with the length of Δx centered at x i = x i−1∕2 + Δx∕2. Time discretization is performed using a forward Euler integration method. Finally, first-order Gudonov-type schemes are used to numerically solve systems (1) and (15) by where W n i and q(W n i , x i ) are approximate averages of the conservative variables and the source terms within the ith cell at the time instant t n ∶= nΔt, respectively. Variables W n i and W n i+1 are approximations of the conservative variables at the left and right sides of the interface x i + 1/2 at the time instant t n , respectively. Moreover,  (⋅, ⋅) is the scheme-specific numerical flux function. Various numerical flux functions have been introduced in the literature. 11,18,21,[29][30][31] As a case study, the flux function for Rusanov scheme is defined as follows where n ∶= (W n i+1 , W n i ) is half of the absolute value of the largest eigenvalue of the Jacobian matrix of the system of Equations (1) or (15) (the Jacobian matrix is (f (w))∕ w in these equations). This state-dependent eigenvalue is is calculated locally at the left and right side of each interfaces. For instance, for system (1), and for system (15), with n i denoting the speed of sound in the mixture of liquid and gas. 21 In the case of no-slip, that is, K = 1, S = 0, the speed of sound in the mixture, n i , can be analytically written as follows, known as the Wood or Wallis speed of sound: 32 .
When slip occurs between the two phases, computing the analytical sound velocity, due to the effect of slippage between the two phases and its effect on the wave propagation speed, is mathematically involved, if not impossible. Thus, simplified surrogates have been suggested in Reference 21 for cases with g g ≪ l l and 0 < g < 1, such as .
These surrogates are not exact and may lead to inaccurate solutions. For this reason, only the case of no-slip is considered in this article. The reader is referred to References 21,29 for a detailed analysis of the speed of sound in the mixture of the gas and liquid.

Modified Rusanov scheme
In the presence of source terms, the numerical solutions obtained by classical finite-volume schemes may drift from the actual solution depending on the contribution of the source terms to the solution. Resolving the issue of generating nonphysical steady-state solutions thus requires further adjustments of the scheme by considering the effects of the source terms, leading to the definition of a WB scheme. By definition, a WB scheme preserves the actual steady-state equation. 8 In the following, we introduce a framework for an approximately WB scheme and provide the motivation for choosing such a framework. This framework differs from the one in (20) in one crucial aspect: the effects of the source term q(W n i , x i ) are incorporated in the input arguments of the numerical flux function  (⋅, ⋅). The proposed modification is inspired by the work in Reference 9. The applicability of that work is, however, limited to basic shallow water equations. This article extends the introduced framework in that article to more advanced and generic models, such as the isothermal Euler equations and DFM.
The proposed structure for the WB scheme: The proposed WB solver of the PDEs (1) and (15) has the following structure where  (⋅, ⋅) can be any numerical flux function and specifically for this study it is defined in (21) as in Rusanov scheme. Moreover, W * ,n i±1 , henceforth called the intermediate variable, satisfies a consistency condition defined in References 9,35. In order to define W * ,n i±1 uniquely, complementary equations alongside the consistency condition, which are source-and model-dependent, are defined for an advection equation, for Euler Equations (1) and for the DFM (15) in Sections 3.2.1,3.2.2, and 3.2.3, respectively. Specifically, the computational steps for the numerical simulation of Euler equations and the DFM are summarized in Algorithms 1 and 2. In the following, we provide a motivation for the proposed structure in (26).
Motivation. Comparing the newly proposed scheme in (26) with the classical one in (20) reveals that the effect of the source terms in (26) are hidden in the intermediate variables while in (20) the source terms are simply integrated over a grid cell. The latter leads to an erroneous steady-state prediction, which originates from the underlying nature of classical schemes that approximate the source terms with a simple average. The problem in predicting a wrong steady-state solution is not resolved by using a different time integrator. To resolve this issue, we embed the effect of the source terms in the intermediate variables in (26) while the intermediate variables satisfy some algebraic constraints derived from steady-state equations. This aspect is the main difference between this study and other studies in literature, 9,10 where the intermediate variables in the approximate Riemann solver have been defined such that the scheme becomes approximately WB. Instead, we modify an existing scheme to make it approximately WB.
Consistency. When manipulating the scheme to contain the effect of source terms, some properties of the scheme should remain intact. Most importantly, the modified scheme should still simulate the transients accurately. In addition, when reaching the steady-state solution, the analytical steady-state solution should be approximated as accurately as possible. Besides, the scheme should always be consistent, that is, , which is the case for the original Rusanov scheme in (21). This is also the case for the proposed WB scheme as we are not changing the definition of the flux function itself, but only modify its input arguments. Another consistency condition, defined in References 9,35, states that the average of the conservative variables obtained by the scheme should be equal to the average over the same cell of the exact solution of the Riemann problem over a length of Δx. If we focus on the spatial interval of (x i , x i + 1 ), this consistency condition imposes the following equality: where W  (W n i , W n i+1 ) gives the conservative variables at the time instant t n + 1 obtained from the exact solution of the Riemann problem at the interface x i + 1/2 , which is dependent on the solutions at the neighboring cell of this interface at time instant t n . Moreover, W n+1 i is obtained from (20) together with (21). To compute the right-hand side of (27), the exact solution of the Riemann problem should be defined over the spatial domain of (x i , x i + 1 ) and in the temporal domain (t n , t n + Δt). It should be noted that the exact solution varies continuously over the spatial and temporal coordinate. For obtaining the above integral for the exact solution of the Riemann problem, one can use the following equivalent equation: 35 1 Δx where f (⋅) is the mathematical flux function as in (1) or (15). The integral of the numerical solution on the left-hand side of (27) depends on the order of the accuracy of the scheme, that is, how the solution changes within a grid cell (for first-order accurate schemes, the solution is constant within a grid cell). We proceed with the computation below by considering first-order accurate schemes. This changes the left side of (27) to: 1 Δx By embedding (20) into (29), we obtain: Now, we propose to accommodate the effects of the source terms q into those input arguments of the numerical flux functions that are not in the neighborhood of the interface x i + 1/2 (in this case, W i + 2 and W i − 1 and change their subscripts to enable the solution locally at each interface). This leads to the definition of the intermediate variables, W * ,n i+1 and W * ,n i , that encompass the effect of the source terms. Embedding these into (30) leads to Finally, the consistency equations at each interface are obtained by equating (31) and (28): Equation (32) is the consistency equation that we consider for the proposed scheme (26).  (31), contain the effect of the source terms, q(W n i , x i ) and q(W n i+1 , x i+1 ), and the effect conservative variables, W n i+2 and W n i−1 . These variables are computed based on the governing equations of the physical phenomenon and also the source terms, which will be explained in Sections 3.2.1-3.2.3 for each case study.
Remark 8. In this study, we focus only on first-order schemes. To upgrade the scheme to second-order, MUSCL approaches (see, e.g., Reference 36) can be followed, which will be the topic of future studies.
The main question is how to approximate the integral in (32) which contains the exact solution to the Riemann problem. This can be resolved by exploiting the governing PDEs that encompass the exact solution, which will be explained in Sections 3.2.2 and 3.2.3. Moreover, (32) contains many unknowns and the equation is not uniquely solvable in isolation. All the complementary equations, necessary to obtain the intermediate variables such as W * ,n i uniquely, will be introduced later for each system of equations in Sections 3.2.2 and 3.2.3.
The unknown terms in (32) are the intermediate variables and the average of the source term (the integral term in (32)). The average of source terms should be determined such that the numerical steady-state solution approximately recovers the analytical one. Since we do not have the exact solution to the Riemann problem, we can exploit the original PDEs to find the average of the source terms. To this end, this average will be approximated by exploiting the algebraic relations originated from the system of Equations (1) and (15), which will be clarified in Sections 3.2.2 and 3.2.3. Henceforth, the integral 1∕Δt 1∕Δx It should be noted that this treatment of the scheme is dependent on the nature of the source terms, which is explained later in the aforementioned sections. The intermediate variables should be obtained by considering two essential properties that are required to be satisfied by the WB solver: consistency with the actual system and the WB property.
The consistency with the actual system of equations is ensured by satisfying algebraic relations originating from the steady-state model. These algebraic relations are obtained in the following section for the isothermal Euler equations and the DFM. For the WB property, from (26), we deduce that the solution is stationary, W n+1 i = W n i , if the solution is assumed constant within each grid cell (the underlying scheme is first-order accurate) and W * ,n i+1 = W n i and W * ,n i = W n i+1 . Therefore, we seek W * ,n i+1 = W n i and W * ,n i−1 = W n i as soon as W n i and W n i+1 define a steady state. Then, the right-hand side of (32) becomes zero and the left-hand side resembles the steady-state model at the discrete level. Here, the pair (W i , W i + 1 ) is said to define a steady state if Equations (5) and (19)  The idea proposed in Reference 14 consists of modifying the effect of the source term by knowing the difference between the numerical steady-state solution and the analytical one. This means that finding the averaged contribution of the source terms requires finding the analytical steady-state solution, which is, however, challenging and expensive due to the BVP structure of the steady-state problem. Instead, one can use algebraic relations that are valid during the steady state without prior knowledge of the steady-state solution itself. Now, the methodology introduced in References 9,10 is employed and modified for any general scheme and applied to Euler equations and the DFM. But first, in the following, we prove that the proposed scheme leads to an exact WB solution for an advection equation.
Remark 9. In this study, we investigate laminar friction characterizations. In general, the approach in this article is applicable to turbulent friction functions as well. However, with turbulent friction, the analysis is highly demanding and the numerical solution is generally hard to obtain, even in the case of a classical numerical solver. As a result, the WB solution would be even more complex.

An advection equation with a source term
As Euler equations and DFM are coupled and the corresponding source terms are sometimes nonlinear with respect to the conservative variables, the analytical assessment of the performance of the proposed scheme on these equations is cumbersome if not impossible. Therefore, we provide the assessment for a simple, though relevant, test case, a scalar PDE governing an advection phenomenon as below: where we set q(w) = w for simplicity of the assessment. (33). The numerical solution of this system for q(w) = w obtained by the solver (26) with numerical flux function in (21) and the intermediate variables satisfying (32) and Δt∕Δx < 1 leads to zero error in approximating the steady-state solution while the solver (20) together with (21) yields an erroneous steady-state solution.

Theorem 1. Consider the advection equation in
Proof. Applying the Rusanov scheme (21) yields ( n i = 1∕2) The consistency condition (32) adapted for this advection equation at any interface is obtained as follows Embedding (34) into the above equation gives, Now, we have to exploit (33) in the steady-state condition to compute Q, as below. In steady state, it holds that Moreover, as we know the source term, for the steady-state solution it holds that Therefore, we can find the expression for QΔx = W i (e Δx − 1). Using the Taylor expansion of e Δx , it can be verified that QΔx is consistent with W i Δx. Then, at each interface, (36) yields Due to the specific form of the advection equation and the Rusanov scheme, the calculation of W * ,n i+1 is not required. The proposed WB solver (26) is repeated here: When reaching the exact steady-state profile, the relation W n+1 i = W n i should hold. We define, = W n+1 i − W n i and we compute for the advection equation with the classical way of dealing with source term as in (20) and the proposed WB way as in (26).
For the WB solver (26) with the Rusanov scheme (21), it holds that By incorporating (39), we obtain Due to the stability of the scheme and the positive numerical diffusion coefficient (because of the CFL condition Δt∕Δx < 1), 37 as n → ∞, → 0. Therefore, W n i = W n i−1 , which corresponds to the actual steady-state solution (38). This shows that the new proposed scheme will lead to zero error at the actual steady-state solution and if the solver starts from the analytical steady-state solution, it remains there. In addition, it can be inferred that the only solution of the WB scheme that yields W n+1 i = W n i is the analytical steady-state solution. Next, we prove that the classical scheme does not preserve the solution. This is proved by contradiction.
Recalling the classical solver (20) with Rusanov scheme (21) as below: we obtain Then, if the solver starts from the analytical steady-state solution, that is, W i = W 0 e iΔx , we obtain, Clearly, the right-hand side of (45) is nonzero, meaning that W i = W 0 e iΔx is not the steady-state solution of (43). This contradiction completes the proof. The error in steady state approximation tends to zero only by making the spatial and temporal grid size (Δx and Δt) smaller. ▪ Supported by Theorem 1, we project that the proposed scheme also leads to better results for coupled equations such as Euler equations and DFM.

Modified scheme for the isothermal Euler equations
To compute the average source terms Q that satisfy steady-state Equation (5), we exploit the discrete version of the steady-state equations. To this end, we integrate (5) at each interface over ( where [⋅] = (⋅) R − (⋅) L denotes the difference of variables between the right (subscript R) and left (subscript L) side of the interface x i + 1/2 . Considering (46b), there are two unknowns, m 0 and Q. In this article, Q = F + G where F and G are the average of frictional and gravitational source terms. One more equation is thus required to solve this equation. The steady-state solution associated to the full source term Q does not admit an algebraic expression in the presence of both friction and gravity and therefore the source terms should be decomposed into individual source terms. So, instead of finding Q such that (46b) is satisfied, we find F and G satisfying other equations with similar structure to (46b), and then set Q = F + G. It should be noted the Q found in this way might not satisfy (46b) exactly an therefore it leads to some errors in the steady-state solution. We first explain this step for the friction-related terms and then for the gravity contribution to the source term. To find consistent (relevant) source terms F and G even in transient case, we exploit the discrete steady-states equations in case of only friction and only gravity, respectively. These source terms entail defining new parameters in the source terms, which converge to the corresponding steady values when the system is reaching the steady condition.

Frictional source terms
Considering only laminar friction (3a) in (5b), to define an average friction source term that is consistent with the steady-state equations, we set: where m can be interpreted as an average of the momentum over the left and right side of the interface and it should be defined such that it converges to m 0 as the solution reaches the steady state. To this end, the following description for m is assumed, This average indeed ensures that when the system reaches its steady state, m L = m R , then m = m L = m R = m 0 . By integrating (47) over the interval of (x i , x i + 1 ), we obtain: When reaching steady-state, m = m 0 and (49) becomes a discrete version of (8) and the same conditions for the availability of the solution holds as in Lemma 1 (with x 0 = x L , x = x R ,̂0 = L , and = R ). Finally solving (49) gives m 0 , and by substituting this m 0 into (46b) in case of only laminar friction as the source term, we obtain Now, by substituting m 0 from (49) into (50), we obtain This source term should be equivalent to (consistent with) the actual friction source term, which is stated below.

Proposition 1. Under the assumption of smooth (steady and transient) solution, F obtained in (51) is consistent with the actual friction defined in (3a).
Proof. For smooth solutions, there exists i and j with Substitution of (52) in (51) yields: For smooths solutions, by increasing the number of the cells, the values of R , L , i , j all converge to a single value, let us say . Then, (53) changes to: In smooth solutions, it holds that m∕ = u. This proves that under the assumption of smooth solutions, F, obtained from (51), indeed approximates (3a). ▪

Gravitational source terms
For the case of only gravitational source term (3b), (5b) changes to: By integrating (55) over the interval of (x i , x i + 1 ), we obtain Notably, (56) is a discrete version of (12) and the same conditions for the availability of the solution holds as in Lemma 2 (with x 0 = x L , x = x R ,̂0 = L , and = R ). After obtaining m 0 from the above equation, G can also be computed from (46b) as follows: Finally, we obtain the following expression for G: Now, we study the equivalence of this source term with the actual gravitational source term. (58) is equivalent to (3b).

Proposition 2. Under the assumption of smooth (steady and transient) solution, G in
) .
Increasing the number of the cells for smooth solutions yields L , R → and therefore (59) changes to: This proves that under the assumption of smooth solution, the term G indeed approximates (3b Remark 11. Note that the expressions of the averaged source terms (51) and (58) have been obtained by considering W L and W R satisfying steady-state models. Since these expressions only depend on the left and right states and it has been proved in Propositions 1 and 2 that the averaged source terms are equivalent to the actual source terms, it is relevant to extend the usage of the averaged source terms to the case where these states do not define a steady state, and actually use the WB expressions (51) and (58) for all W L and W R .
The average source term is then specified as Q = F + G. After obtaining the average source term, the modified scheme can be completed by calculating the intermediate values m * L , m * R , * L , * R needed in the modification of the scheme. To this end, by rewriting the consistency conditions (32) for the modified Rusanov scheme (26) and for (1), we have: As we are focusing on subsonic scenarios, L and R are dominantly governed by the sound velocity apparent from (22) as u ≪ c, not by the state values, we can consider ∶= L = R . Although this assumption leads to some errors in the end, it will significantly simplify the calculations and the nonlinear equations can be solved more easily by Newton-based methods. Due to the specific form of the steady-state solution (constant momentum over the spatial domain) and after, 35 we set m * ∶= m * R = m * L (this choice is also used in case of transients and helps to satisfy the steady equations more easily). Therefore, the system of algebraic relations in (61) simplifies to * Still one more equation is needed to be able to compute the variables uniquely. This last equation should be defined such that when we are on the steady-state profile, the intermediate variables satisfy the steady-state equation (and also should be usable during transients). To do so, we suggest an equation which can be used both at the steady-states and transients. 9 This is carried out by adapting the steady Equation (46b) as follows: The intermediate values should also satisfy the last relation in the above equation, 9,10 meaning that: Now, Equations (62)

Modified scheme for the DFM
Integrating system (19) at each interface over the interval (x i , x i + 1 ), we have: Three unknowns, m l0 , m g0 , and Q, are present in (65c). Two more source-specific equations are thus needed to be coupled with (65c) to uniquely find these unknowns. Since the additional equations are source-dependent, friction and gravity are treated separately. Similar arguments as put forward for the isothermal Euler equations also hold here. So, instead of finding Q such that (65c) is satisfied, we find F and G satisfying other conditions, and then set Q = F + G. It should be noted the Q found in this way might not satisfy (65c) exactly. We will first do the decomposition for the friction-related terms and then for the gravity contribution to the source term. Due to the complicated nature of the DFM, the following assumption is made only to attain the average source terms.
Assumption 1 (39). Only to obtain the average source terms, it is assumed that at interfaces of the discrete DFM as in (19), the volume and mass composition of the mixture do not change.
This assumption is approximately valid when smooth solutions are considered; otherwise, this assumption is less accurate. Assumption 1 results in constant volumetric fraction and mass fraction of each phase at the interface, respectively. This implies that at each interface, the following conditions hold: This assumption may restrict the applicability of the method to systems with smooth solutions; however, the analysis of the applicability is beyond the scope of this study. Again we note that Assumption 1 is only used for obtaining the average source terms and not applied into the DFM (15) itself.

Frictional source terms
Considering only friction, from (19c), we have: A similar approach to (47) cannot be followed here to find similar algebraic relations. Therefore, simplifying Assumption 1 (Equation 66) is considered. Following Reference 9, rearranging (67) and having the constant variables of (66) in mind, we have: where (⋅) represents an average of the variable over the left and right side of the interface, which is defined similar to (48). After some steps of straightforward computations and integration over each interface, we obtain: where From (70), the term B can be computed. By substituting the value of B in (69), the value of A can be computed. Moreover, for computing F in (65c) in the case of only friction using: the values of m l0 and m g0 are required. By knowing the value of A, taking into account the structure of A specified in (70) and assumptions (66), F is obtained as follows: Finally, we obtain

Proposition 3. Under the assumption of smooth (steady and transient) solutions, F obtained in (72) is approximately equivalent to the actual friction (17a).
Proof. For smooth solution, there exists Substituting (73) and EOS (18c) to (72) and carrying out some straightforward simplification leads to: By increasing the number of grid cells, ( g ) R , ( g ) L , ( g ) i → g . Using this and also substituting the expression of B from (70) leads to: ▪ After obtaining the friction contribution, a similar treatment is applied for finding the contribution from the gravitational term.

Gravitational source terms
If gravity is the only source term, from (19c), we have: Similar to the friction source term and using Assumption 1 (Equation 66), we obtain: Finally, we obtain:

Proposition 4. Under the assumption of smooth (steady and transient) solution, G in (80) is approximately equivalent to the actual gravitation source term (17b).
Proof. For smooth solution, there exists Substituting (81) and EOS (18c) to (80) and carrying out some straightforward simplification leads to: ) .
By increasing the number of grid cells, ( g ) R , ( g ) L , ( g ) i → g . By using this and also substituting the expression of B from (78) leads to: This proves that G is an approximation of (17b). ▪ Remark 13. Assumption 1 is adopted only to find expressions for the average friction and gravity source terms. These are not applied to the DFM (15) and, from now on, these assumptions are only hidden in F and G and do not affect the structure of the governing equations.
The average source term is now specified as Q = F + G. Following the same steps introduced in the previous section for Euler equations, we can solve the following equations, consisting of three consistency equations, one algebraic steady-state condition and eight closure laws. Following the same line of reasoning for Euler equations, after imposing The short-hand notation, for instance, ( l l ) * R = ( l ) * R ( l ) * R is used to compact the equations. These 12 equations can be solved simultaneously to compute the intermediate primitive variables, The entire procedure is summarized in Algorithm 2. By running Algorithm 2 for both (W L , W R ) at x i + 1/2 and at x i − 1/2 , the variables that should be substituted into (26) are computed. ,̄l,̄g,m l , andm g 3: ComputeF andḠ via (72) and (80), respectively, 4: SetQ =F +Ḡ, 5: Solve the system (1)- (12) simultaneously to obtain the intermediate primitive variables, 6: Set Now, all the required components for implementing the modified scheme (26) are available and numerical simulations can be obtained.

NUMERICAL RESULTS FOR SINGLE-PHASE FLOW
In this section, numerical results of a single-phase flow inside a pipe are shown. First, preservation of the steady-state solution is considered and then, a transient simulation from an initial steady-state to another steady-state is carried out. The values of the parameters involved in system (1) are listed in Table 1 (q p and p R are the volumetric flow rate of the pump at the left boundary and pressure at the right boundary, respectively). Figure 1 shows the computational domain, where n i is given by (22) and N is the number of grid cells and N t is the number of time-steps. For all simulations in this section, we estimate max(|2 n i |) = c + 10 and CFL = 0.99. Then, according to the chosen Δx, Δt is specified by (85).

Preservation of the steady-state solution
In this section, the performance of the WB scheme in terms of the accuracy of the steady-state profile is evaluated in case of zero and nonzero flow. The initial condition in the test cases of this section coincides with the steady-state solution.
Results are also compared against the classical Rusanov scheme. Moreover, we set T = 10 s. After this time instant, the solution varies negligibly over time (the spatial 2-norm of the solution varies less than 0.1% relatively over time), indicating that the solution has reached its steady state.

Zero flow
Results of simulating a steady system with zero flow are shown in Figures 2 and 3. The differences between the steady-state velocity and the numerical velocity of both the WB and the classical schemes at the right boundary are depicted over time in Figure 4. This example bears relevance with the connection scenario commonly performed in the drilling context, which is shutting down the pump, waiting for the drilling liquid to become stagnant, and adding a new stand of pipe to the current configuration before resuming drilling ahead.
In the left plot of Figure 2, the classical scheme gives a linear change of mass flow over the spatial domain in the steady state. This is in contrast with the physical steady-state solution, that is, the mass flow rate should be constant over the spatial length. The slope of the line in that plot is significantly smaller for the WB scheme, close to zero as the physical governing equations show. However, the pressure in the WB scheme deviates slightly from the steady-state solution, in an extent comparable to the classical scheme. As shown in Figure 4, the error for approximating the velocity at the right boundary in the WB scheme decreases as time evolves. Nonetheless, this error remains unchanged for the classical scheme.
Most notably, the classical scheme implies that zero mass enters from the left boundary and approximately 0.004 kg per second exits the pipe from the right boundary (note that in Figure 2, momentum, u, at the right boundary is around 0.5 kg/m 2 s andṁ = uA = 0.5 × 4 d 2 = 0.004 kg/s), meaning that the mass inside the pipe decreases 0.004 kg per second, contrary to the physical governing condition of the test case.

Nonzero flow
For simulation scenario with nonzero flow, we set q p = 2000 L/min. Simulating such a system leads to the results shown in Figures 5 and 6, where the states and the error in approximating the correct steady-state solution are depicted, respectively. Evolution of the error in approximating the velocity at the right boundary for both schemes is shown in Figure 7 over time. This flow scenario is usually present in the pipeline networks while the flow is pumped from one location to another under a constant volumetric flow rate. This can also be observed in drilling with single-phase flow when the rate of penetration is too low.
Analyzing the results of Figures 5 and 6 reveals that, in case of nonzero flow, the WB scheme always outperforms the classical scheme, such that, its error is considerably lower than that of the classical one. The most important feature of the WB scheme is the preservation of the mass flow rate (see the plots related to u in Figures 5 and 6). Figure 7 shows that the WB scheme always remains closer to the analytical steady-state solution, which further confirms the accuracy of this scheme.

Transient simulation scenario
By observing the results of the previous section, it can be inferred that the WB scheme preserves the steady-state solution with a higher accuracy compared with the classical scheme. In this section, we numerically verify that the WB scheme also approximates the correct steady-state solution with a better accuracy during a transient going from one steady-state to another one by changing the inputs of the system. To this end, in the previous test case at t = 10 s, the right boundary pressure changes from p R = 1 bar to p R = 10 bar and the pump flow rate changes from q p = 2000 L/min to q p = 4000 L/min and the dynamics is simulated until T = 100 s. After this time instant, the state variables of the system do not vary with time, that is, the system has reached its steady state. This example shows a set-point change in the pipeline networks or in a drilling operation. Figures 8 and 9 show the initial condition and the final steady-state solution and the numerical steady-state solution of both schemes and the associated error, respectively. Results show that the WB scheme converges to the steady states with an error significantly smaller than the classical scheme. The approximation error of the velocity at the right boundary is depicted in Figure 10.
The WB scheme captures the new steady-state solution with a lower error. The small deviation of the WB solver from the steady-state solution is due to the simplifying assumption made in Assumption 1 and also the simplification mentioned before (62). The difference of the accuracy of both solvers is more striking if we highlight the source contributions by increasing the nominal density of the liquid ( 0 ) or viscosity of the fluid. The same procedure is carried out for the DFM in the following section.

NUMERICAL RESULTS FOR TWO-PHASE FLOW
In this section, numerical results of a two-phase flow inside a pipe, as shown in Figure 11, are shown. Similar to Section 4, preservation of the steady-state solution and then a transient simulation are provided. The values of the parameters involved in system (15) are summarized in Table 2. It should be mentioned thatṁ l,g = l,g l,g u l,g A and p R are the mass flow rate of the liquid and the gas at the left boundary and pressure at the right boundary, respectively, and A is the cross-sectional area of the pipe. Figure 11 shows the computational domain which is a vertical pipe with a constant cross-sectional area. For all simulations in this section, we set L = 1000 m and Δx = 10 m. For the temporal discretization, we use the CFL definition (85) by estimating max(|2 n i |) = c l and imposing CFL = 0.99. As it has been illustrated in Reference 21, the speed of sound in the mixture appeared in (23) affects the solution significantly. As we want to focus on the performance of the WB scheme compared with the classical scheme, we do not

Preservation of the steady state
For the case of zero flow in a two-phase system in a vertical pipe, due to the density difference of the two phases, the gaseous phase migrates up the pipe and the liquid goes down along the pipe in a nonhorizontal pipe; then, we cannot have a mixture of the gas and liquid for such a pipe. Moreover, the slip law does not permit the separation of the gas and the liquid since the velocity of the gas and liquid have to be equal. Thus, the system does not have any solution for stationary case of a nonhorizontal pipe containing two-phase flow governed by the DFM without any slip. For the horizontal pipe, both the gravitational and the friction term are zero and no source term affects the solution. Hence, in this section, unlike the single-phase flow, only the case of nonzero flow is studied for assessing the ability of the scheme to retain the equilibrium profile of the system. Again, the initial condition for the test cases in this section is the steady-state solution of the system.

Nonzero flow for a vertical pipe
Parameter values are mentioned in Table 2; besides that, T is set to 100 s. Comparison of the steady-state solution and the numerical solutions obtained from the WB and the classical solvers are depicted in Figure 12 and the difference between these solutions and the steady-state solution is shown in Figure 13. In addition, the time evolution of the error in approximating the velocity at the right boundary is shown in Figure 14. This scenario can occur in the two-phase pipelines as well as the flow inside a well-bore during an underbalanced drilling operation. 40 It can be observed that the WB scheme performs better than the classical scheme in terms of remaining on the steady-state profile when starting from the steady initial condition. The effectiveness of the WB scheme can be better observed in the plots related to mass flow rate of the liquid ( l l u l in Figures 12 and 13). Apparent from these figures, the classical scheme generates a nonphysical steady-state solution. Pressure and velocity predicted by the classical scheme deviate from the actual values significantly. Again, in the steady state, the classical scheme predicts that 0.3 kg of the liquid enters the pipe per second and around 0.37 kg of the liquid exits the pipe, which is in contrast with the physical governing equations. The results of the WB scheme and the actual steady-state solution correspond to each other with a high accuracy. Now, the scheme can be tested under transient simulations by going from one steady initial condition to another steady-state solution.

Transient simulation scenario
Analogous to the case of single-phase flow, boundary conditions are changed to excite the transients of the systems. At t = 1 s, the values of the mass flow rates are doubled instantly and the right boundary pressure changes to p R = 2 bar slowly over 10 s. This set-point change can also happen in the two-phase pipelines and also as a control action to harness the gas migration in a drilling well.
The results acquired at time instant T = 500 s are shown in Figures 15,16 and 17. These figures reveal that the WB scheme approximates the steady-state solution of the system with a higher accuracy compared with the classical scheme. As the system approaches to the steady condition, the errors in all state variables also tend to zero over time.
Results presented in Sections 4 and 5 confirm that the WB schemes yield more accurate results in both preserving and approximating the steady-state solution compared with the classical scheme. Moreover, the modified schemes have also been tested in the case of no source terms (the flow with zero viscosity in a horizontal pipe), which has shown that it produces the same results as the classical scheme (these results are not included in the article). Now, we can analyze the order of accuracy for the proposed scheme.

ORDER OF ACCURACY FOR THE MODIFIED SCHEME
In order to study the error convergence of the proposed scheme, a norm of the difference between the steady-state solution and the numerical steady-state solution generated by the proposed schemes is computed. The error-norm is velocity-based and it is defined as below: where N is the number of the grid cells, T, as mentioned earlier, is the last time instant of the simulation and u ss is the steady-state value of the velocity. In this study, we set r = {2, ∞}, meaning that we calculate the spatial 2−norm and ∞−norm of the difference between the trustworthy and the numerical steady-state solution divided by the spatial 2−norm or ∞−norm of the trustworthy steady-state solution. In this section, we denote e 2 for (86) calculated by r = 2 and e ∞ for  (86) calculated by r = ∞. Simulations are performed for the single-and two-phase flow (u: = u g = u l ) and the number of grid cells is varied to analyze the dependency of (86) on the number of grid cells.

Single-phase flow
The parameter settings in Section 4.1.2 are used for this section. Table 3 illustrates the dependency of the steady-state error measure in (86) on the number of grid cells for both the classical and WB (Table 3) schemes. In addition, the CPU time for carrying out these simulations are also reported in this table. The first-order convergence of the classical Rusanov scheme is clear, 18 as by doubling the number of the grid cells, the error is divided by two. Moreover, the error associated with the WB scheme with 100 grid cells is much less than the error of the classical scheme with 1600 grid cells. By extrapolation, the classical scheme with 12,800 grid cells yields the same accuracy as the WB scheme with 100 grid cells while the WB scheme with 100 grid cells is less expensive than using classical scheme even with 400 grid cells. This clearly shows the superiority of the WB scheme over the classical scheme in correct steady-state calculation. The proposed scheme is close to first-order accurate. To be precise, its accuracy is of order of 0.82 on average (see Figure 18 ). The main reason for this lower convergence rate can be attributed to the assumption ∶= L = R ; it seems that by increasing the number of grid cells, the error related to this assumption does not reduce linearly and therefore decreases the overall convergence rate of the proposed scheme. Another reason can be the point-wise calculation of the reference solution, rather than the average solution over a grid-cell. Approximating the averaged source terms as the summation of the averaged frictional and gravitational source terms and also the fact mentioned in Remark 7 contribute to this behavior. The data shown in Table 3 together with the order of accuracy can be seen in Figure 18. It can be observed that by increasing the number of grid cells, the incremental order of accuracy (the order of accuracy in each step of doubling the number of grid cells) does not follow a specific trend. By going from 100 to 200 grid cells, the error drops even more than a first-order accurate scheme. Then the order of accuracy decreases and then increases.

Two-phase flow
The parameter settings in Section 5.1.1 are used for this section. The effect of increasing the number of the grid cells on capturing the steady-state solution of the DFM is reported in Table 4 together with the CPU time allocated for the simulations.
The first-order accuracy of the classical Rusanov scheme can easily be interpreted from Table 4. Moreover, the error associated with the WB scheme with 100 grid cells is much less than the error of the classical scheme with 1600 grid cells. Similarly, by extrapolation, the classical scheme with 12,800 grid cells generates the same accuracy as the WB scheme with 100 grid cells while the WB scheme with 100 grid cells is less expensive than using classical scheme even with 1600 grid cells. Due to the nonlinearity of the equations in (84), the CPU time for the WB scheme is higher than the classical scheme. This, however, can be alleviated remarkably by embedding the linear equations into the nonlinear ones among (84), so a smaller set of equations has to be solved at each time step. In general, the comparison of the error (86) shows the superiority of the WB scheme over the classical scheme. However, this superiority comes at the expense of reduction in the error convergence rate; the proposed scheme is of order of 0.65 on average (see Figure 19). The main culprits for this reduction in the convergence rate are the assumption of ∶= L = R and the assumption (66). Approximating the averaged source terms as the summation of the averaged frictional and gravitational source terms and also the fact mentioned in Remark 7 contribute to this behavior. The data shown in Table 4 together with the order of accuracy can be seen in Figure 19. As can be observed in this figure, by increasing the number of grid-cells, the incremental order of accuracy increases and becomes closer to 1. This can be attributed to the fact that by increasing the number of grid cells, the conservative and primitive variables vary even more smoothly from one grid cell to another and the set of assumptions (66) becomes more realistic. In the high number of grid cells, the error generated by the assumption ∶= L = R dominates the error. Similarly, point-wise calculation of the reference steady-state solution can contribute to this behavior.

CONCLUSION
In this article, a novel extension of the Rusanov scheme has been proposed to improve the preservation of the steady-state solutions of Euler equations and the DFM. These schemes reduce to the original scheme when there is no source term. The proposed schemes capture the steady-state solution with significantly higher accuracy compared with the classical scheme in the presence of source terms. This is proved for an advection equation with a simple source term. Various test cases of zero and nonzero flow have been carried out and the improved performance of the WB schemes has been shown numerically, both for single-phase and two-phase scenarios. The modification is model-and scheme-dependent, such that a similar approach can be followed for other systems of PDEs solved by different schemes. The WB schemes lead to physical results while this is not the case for the classical finite-volume schemes in the presence of source terms. This can be interpreted from the presented results, especially those for the mass flow rate. Both Euler equations and DFM imply that the mass flow rate along the spatial domain remains constant in the steady condition. This property is preserved with significantly higher accuracy in the proposed WB schemes.