High order finite volume schemes based on reconstruction of states for solving hyperbolic systems with nonconservative products. Applications to shallow-water systems
By Manuel Castro, José M. Gallardo, Carlos Parés
Abstract
This paper is concerned with the development of high order methods for the numerical approximation of one-dimensional nonconservative hyperbolic systems. In particular, we are interested in high order extensions of the generalized Roe methods introduced by I. Toumi in 1992, based on WENO reconstruction of states. We also investigate the well-balanced properties of the resulting schemes. Finally, we will focus on applications to shallow-water systems.
1. Introduction
The motivating question of this paper was the design of well-balanced high order numerical schemes for PDE systems that can be written under the form
where the unknown $w(x,t)$ takes values on an open convex subset $D$ of $\mathbb{R}^{N}$,$F$ is a regular function from $D$ to $\mathbb{R}^{N}$,$\mathcal{B}$ is a regular matrix-valued function from $D$ to $\mathcal{M}_{N \times N}(\mathbb{R})$,${S}$ is a function from $D$ to $\mathbb{R}^{N}$, and $\sigma (x)$ is a known function from $\mathbb{R}$ to $\mathbb{R}$.
System Equation1.1 includes as particular cases: systems of conservation laws ($\mathcal{B}=0$,$S=0$), systems of conservation laws with source term or balance laws ($\mathcal{B}=0$), and coupled systems of conservation laws.
More precisely, the discretization of the shallow-water systems that govern the flow of one layer or two superposed layers of immiscible homogeneous fluids was focused. The corresponding systems can be written respectively as a balance law or a coupled system of two conservation laws. Systems with similar characteristics also appear in other flow models, such as boiling flows and two-phase flows (see Reference13).
It is well known that standard methods that correctly solve systems of conservation laws can fail in solving Equation1.1, especially when approaching equilibria or near to equilibria solutions. In the context of shallow-water equations, Bermúdez and Vázquez-Cendón introduced in Reference2 the condition called conservation property or $\mathcal{C}$-property: a scheme is said to satisfy this condition if it correctly solves the steady state solutions corresponding to water at rest. This idea of constructing numerical schemes that preserve some equilibria, which are called in general well-balanced schemes, has been extended in different ways; see, e.g., Reference3, Reference6, Reference7, Reference8, Reference12, Reference14, Reference17, Reference18, Reference21, Reference22, Reference26, Reference28, Reference29, Reference30, Reference35.
Among the main techniques used in the derivation of well-balanced numerical schemes, one of them consists in first choosing a standard conservative scheme for the discretization of the flux terms and then discretizing the source and the coupling terms in order to obtain a consistent scheme which correctly solves a predetermined family of equilibria. This was the approach in Reference2 where the authors proved, in the context of shallow-water equations, that numerical schemes based on Roe methods for the discretization of the flux terms and upwinding the source term exacly solve equilibria corresponding to water at rest. In Reference12 it was shown that the technique of modified equations can be helpful in the deduction of well-balanced numerical schemes.
This procedure has been succesfully applied to obtain high order numerical schemes for some particular cases of Equation1.1 (see, for instance, Reference4, Reference38 and Reference39). The main disadvantage of this first technique is its lack of generality: the calculation of the correct discretization of the source and the coupling terms depends on both the specific problem and the conservative numerical scheme chosen.
Another technique to obtain well-balanced first order schemes for solving Equation1.1 consists in considering piecewise constant approximations of the solutions that are updated by means of Approximate Riemann Solvers at the intercells. In particular, Godunov’s methods, i.e., methods based on Exact Riemann Solvers, have been used in the context of shallow-water systems in Reference1, Reference9, Reference10, Reference15, Reference21, Reference22. This approach was also used in Reference5, where the flux and the coupling terms of a bilayer shallow-water system were treated together by using a generalized Roe linearization.
If this second procedure is followed, the main difficulty both from the mathematical and the numerical points of view comes from the presence of nonconservative products, which makes difficult even the definition of weak solutions: in general, the product $\mathcal{B}(w)w_x$ does not make sense as a distribution for discontinuous solutions. This is also the case for the product ${S}(w)\sigma _x$ when piecewise constant approximations of $\sigma$ are considered.
A helpful strategy in solving these difficulties consists in considering system Equation1.1 as a particular case of a one-dimensional quasilinear hyperbolic system:
system Equation1.1 can be easily rewritten under this form (see Reference17, Reference18, Reference21, Reference22).
In Reference11, Dal Maso, LeFloch, and Murat proposed an interpretation of non-conservative products as Borel measures, based on the choice of a family of paths in the phases space. After this theory it is possible to give a rigorous definition of weak solutions of Equation1.2. Together with the definition of weak solutions, a notion of entropy has to be chosen as the usual Lax’s concept or one related to an entropy pair. Once this choice has been done, the classical theory of simple waves of hyperbolic systems of conservation laws and the results concerning the solutions of Riemann problems can be extended to systems of the form Equation1.2.
The introduction of a family of paths does not only give a way to properly define the concept of weak solution for nonconservative systems, it also allows us to extend to this framework some basic concepts related to the numerical approximation of weak solutions of conservation laws. For instance, in Reference36 a general definition of Roe linearizations was introduced, also based on the use of a family of paths. In Reference28 a general definition of well-balanced schemes for solving Equation1.2 was introduced. It was shown that the well-balanced properties of these generalized Roe methods depend on the choice of the family of paths. Moreover, this general methodology was applied to some systems of the form Equation1.1 related to shallow-water flows, recovering some known well-balanced schemes, or resulting in new schemes.
The goal of this paper is to obtain the general expression of a well-balanced high order method for Equation1.2 based on the use of a first order Roe scheme and reconstruction of states. The interest of such a general expression is that, once obtained, particular schemes can be deduced for any system of the form Equation1.1, where the numerical treatment of source and coupling terms is automatically derived. To our knowledge, the present work is the first attempt to obtain well-balanced high order numerical schemes following this procedure.
The paper is organized as follows. In Section 2 we give some basic definitions and results about nonconservative systems, Roe linearizations and generalized Roe schemes, for which we will follow Reference28 closely. High order versions of the Roe schemes, based on reconstruction operators, are introduced in Section 3. Next, Section 4 is devoted to the analysis of the well-balanced properties of the high order schemes previously constructed. In Section 6, the WENO method is applied to build the reconstruction operators. Applications to a family of systems that generalize Equation1.1 are presented in Section 5, with particular interest in some shallow-water systems with one and two layers of fluid. Finally, Section 7 contains numerical results to test the performances of our high-order schemes. In particular, the high order well-balanced property is numerically verified.
2. Roe methods for nonconservative hyperbolic systems
where we suppose that the range of $W(x,t)$ is contained inside an open convex subset $\Omega$ of $\mathbb{R}^N$, and $W\in \Omega \mapsto \mathcal{A}(W)\in \mathcal{M}_N(\mathbb{R})$ is a smooth locally bounded map. The system Equation2.1 is assumed to be strictly hyperbolic: for each $W\in \Omega$ the matrix $\mathcal{A}(W)$ has $N$ real distinct eigenvalues $\lambda _1(W)<\cdots <\lambda _N(W)$. We also suppose that the $j$th characteristic field $R_j$ is either genuinely nonlinear:
For discontinuous solutions $W$, the nonconservative product $\mathcal{A}(W)W_x$ does not make sense as a distribution. However, the theory developed by Dal Maso, LeFloch and Murat (Reference11) allows us to give a rigorous definition of nonconservative products, associated to the choice of a family of paths in $\Omega$.
Definition 2.1
A family of paths in $\Omega \subset \mathbb{R}^N$ is a locally Lipschitz map
$\Phi (0;W_L,W_R)=W_L$ and $\Phi (1;W_L,W_R)=W_R$, for any $W_L,W_R\in \Omega$.
(2)
Given an arbitrary bounded set $\mathcal{B}\subset \Omega$, there exists a constant $k$ such that$$\begin{equation*} \bigg | \frac {\partial \Phi }{\partial s}(s;W_L,W_R)\bigg | \leq k|W_L-W_R|, \end{equation*}$$
for any $W_L,W_R\in \mathcal{B}$ and for almost every $s\in [0,1]$.
(3)
For every bounded set $\mathcal{B}\subset \Omega$, there exists a constant $K$ such that$$\begin{equation*} \bigg | \frac {\partial \Phi }{\partial s}(s;W_L^1,W_R^1)-\frac {\partial \Phi }{\partial s}(s;W_L^2,W_R^2)\bigg | \leq K(|W_L^1-W_L^2|+|W_R^1-W_R^2|), \end{equation*}$$
for each $W_L^1,W_R^1,W_L^2,W_R^2\in \mathcal{B}$ and for almost every $s\in [0,1]$.
Suppose that a family of paths $\Phi$ in $\Omega$ has been chosen. Then, for $W\in (L^\infty (\mathbb{R}\times \mathbb{R}^+)\cap BV(\mathbb{R}\times \mathbb{R}^+))^N$, the nonconservative product can be interpreted as a Borel measure denoted by $[\mathcal{A}(W)W_x]_\Phi$. When no confusion arises, we will drop the dependency on $\Phi$.
A weak solution of system Equation2.1 is defined as a function $W\in (L^\infty (\mathbb{R}\times \mathbb{R}^+)\cap BV(\mathbb{R}\times \mathbb{R}^+))^N$ that satisfies the equality
In particular, a piecewise $\mathcal{C}^1$ function $W$ is a weak solution of Equation2.1 if and only if the two following conditions are satisfied:
(i)
$W$ is a classical solution in the domains where it is $\mathcal{C}^1$.
(ii)
Along a discontinuity $W$ satisfies the jump condition$$\begin{equation} \int _0^1\big (\xi \mathcal{I}-\mathcal{A}(\Phi (s;W^-,W^+))\big )\frac {\partial \Phi }{\partial s}(s;W^-,W^+)\, ds=0, \tag{2.2}\cssId{RHg}{} \end{equation}$$
where $\mathcal{I}$ is the identity matrix, $\xi$ is the speed of propagation of the discontinuity, and $W^-$,$W^+$ are the left and right limits of the solution at the discontinuity.
Note that in the particular case of a system of conservation laws (that is, $\mathcal{A}(W)$ is the Jacobian matrix of some flux function $F(W)$) the jump condition Equation2.2 is independent of the family of paths, and it reduces to the usual Rankine-Hugoniot condition:
In the general case, the selection of the family of paths has to be based on the physical background of the problem under consideration. Nevertheless, it is natural from the mathematical point of view to require this family to satisfy some hypotheses concerning the relation of the paths with the integral curves of the characteristic fields. For instance, if $W_L$ and $W_R$ are linked by an integral curve of a linearly degenerate field, the natural choice of the path is a parametrization of that curve, as this choice assures that the contact discontinuity
where $\xi$ is the (constant) value of the corresponding eigenvalue through the integral curve, is a weak solution of the problem, as would be the case for a system of conservation laws.
Due to these requirements, the explicit calculation of the path linking two given states $W_L$ and $W_R$ can be difficult: in most cases, the explicit expression of the solution of the Riemann problem related to the states is needed (see Reference28).
Together with this definition of weak solutions, a notion of entropy has to be chosen, either as the usual Lax’s concept or one related to an entropy pair (see Reference16 for details). Once this choice has been done, the theory of simple waves of hyperbolic systems of conservation laws and the results concerning the solutions of Riemann problems can be naturally extended to systems of the form Equation2.1 (see Reference11).
Some of the usual numerical schemes designed for conservation laws can be adapted to the discretization of the more general system Equation2.1. This is the case of Roe schemes (see Reference31): in Reference36 a general definition of Roe linearization was introduced, based again on the use of a family of paths.
Definition 2.2
Given a family of paths $\Psi$, a function $\mathcal{A}_\Psi \colon \Omega \times \Omega \to \mathcal{M}_N(\mathbb{R})$ is called a Roe linearization of system Equation2.1 if it verifies the following properties:
(1)
For each $W_L,W_R\in \Omega$,$\mathcal{A}_\Psi (W_L,W_R)$ has $N$ distinct real eigenvalues.
(2)
$\mathcal{A}_\Psi (W,W)=\mathcal{A}(W)$, for every $W\in \Omega$.
(3)
For any $W_L,W_R\in \Omega$,$$\begin{equation} \mathcal{A}_\Psi (W_L,W_R)(W_R-W_L)=\int _0^1\mathcal{A}(\Psi (s;W_L,W_R))\frac {\partial \Psi }{\partial s}(s;W_L,W_R)\,ds. \tag{2.5}\cssId{Roe-prop-g}{} \end{equation}$$
Note again that if $\mathcal{A}(W)$ is the Jacobian matrix of a smooth flux function $F(W)$,Equation2.5 is independent of the family of paths and reduces to the usual Roe property:
Once a Roe linearization $\mathcal{A}_\Psi$ has been chosen, in order to construct numerical schemes for solving Equation2.1, computing cells $I_i=[x_{i-1/2},x_{i+1/2}]$ are considered; let us suppose for simplicity that the cells have constant size $\Delta x$ and that $x_{i+\frac {1}{2}}=i\Delta x$. Define $x_i=(i-1/2)\Delta x$, the center of the cell $I_i$. Let $\Delta t$ be the constant time step and define $t^n=n\Delta t$. Denote by $W_i^n$ the approximation of the cell averages of the exact solution provided by the numerical scheme, that is,
Then, the numerical scheme advances in time by solving linear Riemann problems at the intercells at time $t^n$ and taking the averages of their solutions on the cells at time $t^{n+1}$. Under usual CFL conditions, it can be written as follows (see Reference28):
$$\begin{equation} W_i^{n+1}=W_i^n-\frac {\Delta t}{\Delta x}\big (\mathcal{A}^+_{i-1/2}(W_i^n- W_{i-1}^n) +\mathcal{A}^-_{i+1/2}(W_{i+1}^n-W_i^n)\big ). \tag{2.7}\cssId{Roe}{} \end{equation}$$ Here, the intermediate matrices are defined by
and $\mathcal{K}_{i+1/2}$ is an $N\times N$ matrix whose columns are associated eigenvectors.
In the particular case of a system of conservation laws, Equation2.7 can be written under the usual form of a conservative numerical scheme. First, the numerical flux is defined by
The best choice of the family of paths $\Psi$ appearing in the definition of Roe linearization is the family $\Phi$ selected for the definition of weak solutions. In effect, Roe methods based on the family of paths $\Phi$ can correctly solve discontinuities in the following sense: let us suppose that $W_i^n$ and $W_{i+1}^n$ can be linked by an entropic discontinuity propagating at speed $\xi$; then, from Equation2.5 and Equation2.2 we deduce that
i.e., $\xi$ is an eigenvalue of the intermediate matrix and $W_{i+1}^n -W_i^n$ is an associated eigenvector. As a consequence, the solution of the linear Riemann problem corresponding to the intercell $x_{i+1/2}$,
Both solutions consist of a discontinuity linking the states and propagating at velocity $\xi$.
Nevertheless, the construction of these schemes with $\Psi =\Phi$ can be difficult or very costly in practice. In this case, a simpler family of paths $\Psi$ has to be chosen as the family of segments
In Reference28 it was remarked that, in this case, the convergence of the numerical scheme can fail when the weak solution to approach involves discontinuities connecting states $W^-$ and $W^+$ such that the paths of the families $\Phi$ and $\Psi$ linking them are different.
As in the case of a system of conservation laws, the scheme Equation2.7 has to be used with a CFL condition:
$$\begin{equation*} \max \left\{ |\lambda ^{i+1/2}_{l}|,\ 1 \leq l \leq N,\ i \in \mathbb{Z} \right\}\frac {\Delta t}{\Delta x}\leq \gamma , \end{equation*}$$
with $0<\gamma \leq 1$. An entropy fix technique, as the Harten-Hyman one (Reference24, Reference25), also has to be included.
3. High order schemes based on reconstruction of states
there are several methods to obtain higher order schemes based on the use of a reconstruction operator. In particular, methods based on the reconstruction of states are built using the following procedure: given a first order conservative scheme with numerical flux function $G(U,V)$, a reconstruction operator of order $p$ is considered, that is, an operator that associates to a given sequence $\{W_i\}$ two new sequences, $\{W_{i+1/2}^-\}$ and $\{W_{i+1/2}^+\}$, in such a way that, whenever
Once the first order method and the reconstruction operator have been chosen, the method of lines can be used to develop high order methods for Equation3.1: the idea is to discretize only in space, leaving the problem continuous in time. This procedure leads to a system of ordinary differential equations which is solved using a standard numerical method. In particular, we assume here that the first order scheme is a Roe method.
Let $\overline{W}_j(t)$ denote the cell average of a regular solution $W$ of Equation3.1 over the cell $I_i$ at time $t$:
where $W_i(t)$ is the approximation to $\overline{W}_i(t)$ provided by the scheme, and $W_{i+1/2}^\pm (t)$ is the reconstruction associated to the sequence $\{W_j(t) \}$.
Let us now generalize this semi-discrete method for a nonconservative system Equation2.1. Observe that, in Section 2, the key point to generalize both the Rankine-Hugoniot condition Equation2.3 and the Roe property Equation2.6 to system Equation2.1 was to replace a difference of fluxes by an integral along a path. Let us apply the same technique here. First of all, as the first order scheme is a Roe method, using the equalities Equation2.8 and Equation2.9 (replacing $W_i^n$ and $W_{i+1}^n$ by $W_{i+1/2}^-(t)$ and $W_{i+1/2}^+(t)$, respectively) we can rewrite Equation3.2 as follows:
where the intermediate matrices are defined by means of a Roe linearization based on a family of paths $\Psi$ and $P^t_i$ is again a regular function satisfying Equation3.4.
Remark 3.1
It is important to note that for conservative problems, the numerical scheme Equation3.6 is equivalent to the conservative numerical scheme Equation3.3 if, and only if, the integral term is computed exactly. However, the formulation Equation3.6 is useless when working with conservative problems, as we would get involved with a more complex expression of the numerical scheme. The numerical scheme Equation3.6 is useful only for problems with nonconservative products, as it allows us to deduce numerical schemes for particular problems, using numerical quadratures if necessary.
There is an important difference between the conservative and the nonconservative case: in the conservative case the numerical scheme is independent of the functions $P^t_i$ chosen at the cells, but this is not the case for nonconservative problems. As a consequence, while the numerical scheme Equation3.2 has the same order of the reconstruction operator, in the case of the scheme Equation3.6 it seems clear that, in order to have a high order scheme, together with a high order reconstruction operator, the functions $P^t_i$ and their derivatives have to be high order approximations of $W(\cdot , t)$ and its partial derivative $W(\cdot , t)_x$.
In practice, the definition of the reconstruction operator gives the natural choice of the function $P^t_i$, as the usual procedure to define a reconstruction operator is the following: given a sequence $\{W_i\}$ of values at the cells, first an approximation function is constructed at every cell $I_i$, based on the values of $W_i$ at some of the neighbor cells (the stencil):
for some natural numbers $l,r$. These approximations functions are calculated by means of an interpolation or approximation procedure. Once these functions have been constructed, $W_{i+1/2}^-$ (resp. $W_{i+1/2}^+$) is calculated by taking the limit of $P_i$ (resp. $P_{i+1}$) to the left (resp. to the right) of $x_{i+1/2}$. If the reconstruction operator is built following this procedure (as we will assume in the sequel), the natural choice of $P^t_i$ is
Let us now investigate the order of the numerical scheme Equation3.6. Note first that, for regular solutions $W$ of Equation2.1, the cell averages at the cells $\{\overline{W}_j(t)\}$ satisfy
Thus, Equation3.6 is expected to be a good approximation of Equation3.7. This fact is stated in the following result:
Theorem 3.2
Let us assume that $\mathcal{A}$ is of class $\mathcal{C}^2$ with bounded derivatives and $\mathcal{A}_\Psi$ is bounded. Let us also suppose that the $p$-order reconstruction operator is such that, given a sequence defined by
for every smooth enough solution $W$,$W_{i+1/2}^\pm (t)$ being the associated reconstructions and $P^t_i$ the approximation functions corresponding to the sequence
The equality Equation3.8 is easily deduced from the above equalities.
■
Remark 3.3
For the usual reconstruction operators one has $r \leq q \leq p$, and thus the order of Equation3.6 is $r + 1$ for nonconservative systems and $p$ for conservation laws. Therefore a loss of accuracy can be observed when a technique of reconstruction giving order $p$ for systems of conservation laws is applied to a nonconservative problem.
4. Well-balanced property
In this paragraph we investigate the well-balanced properties of schemes of the form Equation3.6. Well-balancing is related with the numerical approximation of equilibria, i.e., steady state solutions. System Equation2.1 can only have nontrivial steady state solutions if it has linearly degenerate fields: if $W(x)$ is a regular steady state solution it satisfies
and then 0 is an eigenvalue of $\mathcal{A}(W(x))$ for all $x$ and $W'(x)$ is an eigenvector. Therefore, the solution can be interpreted as a parametrization of an integral curve of a linearly degenerate characteristic field whose corresponding eigenvalue takes the value 0 through the curve. In order to define the concept of well-balancing, let us introduce the set $\Gamma$ of all the integral curves $\gamma$ of a linearly degenerate field of $\mathcal{A}(W)$ such that the corresponding eigenvalue vanishes on $\Gamma$. According to Reference28, given a curve $\gamma \in \Gamma$, a numerical scheme is said to be exactly well-balanced (respectively well-balanced with order$k$) for $\gamma$ if it solves exactly (respectively up to order $O(\Delta x^k)$) regular stationary solutions $W$ satisfying $W(x) \in \gamma$ for every $x$. The numerical scheme is said to be exactly well-balanced or well-balanced with order $k$ if these properties are satisfied for any curve $\gamma$ of $\Gamma$ (see Reference28 for details).
In the cited article, it has been shown that a Roe scheme Equation2.7 based on a family of paths $\Psi$ is exactly well balanced for a curve $\gamma \in \Gamma$ if, given two states $W_L$ and $W_R$ in $\gamma$, the path $\Psi (s;W_L,W_R)$ is a parametrization of the arc of $\gamma$ linking the states. In particular, if the family of paths $\Psi$ coincides with the one used in the definition of weak solutions $\Phi$, the numerical scheme is exactly well balanced. On the other hand, the numerical scheme is well balanced with order $p$ if $\Psi (s;W_L,W_R)$ approximates with order $p$ a regular parametrization of the arc of $\gamma$ linking the states. In particular, a Roe scheme based on the family of segments Equation2.10 is well balanced with order 2. Moreover, it is exactly well balanced for curves of $\Gamma$ that are straight lines.
The definition of a well-balanced scheme introduced in Reference28 can be easily extended for semi-discrete methods.
Definition 4.1
Let us consider a semi-discrete method for solving Equation2.1:
where $\mathbf{W}(t) = \{ W_i(t) \}$ represent the vector of approximations to the cell averages of the exact solution, and $\mathbf{W}_0 = \{W^0_i \}$ the vector of initial data. Let $\gamma$ be a curve of $\Gamma$. The numerical method Equation4.1 is said to be exactly well balanced for $\gamma$ if, given a regular stationary solution $W$ such that
Finally, the semi-discrete method Equation4.1 is said to be exactly well balanced or well balanced with order$p$ if these properties are satisfied for every curve $\gamma$ of the set $\Gamma$.
For the particular case of the numerical schemes based on reconstruction of states Equation3.6 we have
where $W_{i+1/2}^\pm$ are the reconstructions associated to the sequence $\mathbf{W}$ and $P_i$ the corresponding approximation functions. Hereafter, we give two results concerning the well-balanced property of this scheme, but first we introduce a new definition.
Definition 4.2
A reconstruction operator based on smooth approximation functions is said to be exactly well balanced for a curve $\gamma \in \Gamma$ if, given a sequence $\{W_i\}$ in $\gamma$, the approximation functions satisfy
Let $\gamma$ belong to $\Gamma$. Let us suppose that both the generalized Roe method and the reconstruction operator chosen are exactly well balanced for $\gamma$. Then the numerical scheme Equation3.6 is also exactly well balanced for $\gamma$.
Proof.
Let $W$ be a regular stationary solution satisfying
On the other hand, using Equation4.2, $P_i$ can be understood as a parametrization of an arc of $\gamma$, which is an integral curve of a linearly degenerate field whose corresponding eigenvalue is zero. Therefore,
Here, $\sigma (x)$ is a known function from $\mathbb{R}$ to $\mathbb{R}$,$F$ is a regular function from $\Omega \times \mathbb{R}$ to $\mathbb{R}^{N}$,$\Omega$ is an open convex subset of $\mathbb{R}^{N}$,$\mathcal{B}$ is a regular matrix-valued function from $\Omega \times \mathbb{R}$ to $\mathcal{M}_N(\mathbb{R})$, and $\widetilde{S}$ is a function from $\Omega$ to $\mathbb{R}^{N}$. We can assume without loss of generality that $\widetilde{S}$ has the form
System Equation5.1 includes as particular cases systems of conservation laws ($\mathcal{B}=0$,$S=0$) whose flux function may depend on $x$ via the function $\sigma$, systems of conservation laws with source term or balance laws ($\mathcal{B}=0$), or coupled systems of conservation laws as defined in Reference5. In this latter case, $\mathcal{J}$ is block-diagonal and the blocks of $\mathcal{B}$ corresponding to the nonzero diagonal blocks of $\mathcal{J}$ are zero.
Following the idea developed in Reference17, Reference18 for conservation laws with source terms, if we add to Equation5.1 the trivial equation
and associated eigenvectors $R_j(\widetilde{W})$,$j=1,\dotsc ,N$. If these eigenvalues do not vanish, Equation5.2 is a strictly hyperbolic system: $\mathcal{\widetilde{A}}(\widetilde{W})$ has $N+1$ distinct real eigenvalues
Clearly, the $(N+1)$-th field is linearly degenerate and, for the sake of simplicity, we will suppose that it is the only one. The integral curves of the linearly degenerate field are given by those of the ODE system
Note that, in this case, the set $\Gamma$ defined in the previous section is simply the set of all the integral curves of the linearly degenerate field, as the corresponding eigenvalues always take the value 0. Let us illustrate in this case the relation between these integral curves and the stationary solutions. Let $\gamma$ be an integral curve of the linearly degenerate field and let us suppose that it can be described implicitly by a system of $N$ equations:
As $\sigma$ is a known function, for every $x$,Equation5.3 is a system of $N$ equations with $N$ unknowns $w_1, \dotsc , w_N$. The stationary solutions associated to the curve $\gamma$ are obtained by searching solutions $\{w_1(x), \dotsc , w_N(x)\}$ of system Equation5.3 which depend smoothly on $x$.
For the definition of weak solutions of system Equation5.2 and the choice of the family of paths, we refer the interested reader to Reference28 and the references therein. Let us only mention that the complete definition of the path linking two states is not easy, as it requires the explicit knowledge of the solution of the corresponding Riemann problem. Therefore, the construction of Roe schemes based on the family of paths used in the definition of weak solutions is, in general, a difficult task.
Thus we consider the general case in which the family of paths $\widetilde{\Psi }$ used for the construction of Roe matrices is different to that used in the definition of weak solutions. In particular, in the applications the family of segments Equation2.10 has been considered.
Let us suppose that, for any fixed value of $\sigma$, Roe matrices can be calculated for the system of conservation laws corresponding to $\mathcal{B}=0$ and $S=0$, i.e., we assume that, given ${W}^n_i$,${W}^n_{i+1}$ and $\sigma$, we can calculate a matrix $\mathcal{J}^{\sigma }_{i+1/2}$ such that
Let us also suppose that it is possible to calculate a value $\sigma _{i+1/2}$ of $\sigma$, a $N \times N$ matrix $\mathcal{B}_{i+1/2}$, and a vector $S_{i+1/2}$, such that the following identities hold:
where $\mathcal{K}_{i+1/2}$ is the $N\times N$ matrix whose columns are the eigenvectors $R_{i+1/2,1}$, …, $R_{i+1/2,N}$ and $\mathop {\operatorname {sgn}}\nolimits (\mathcal{L})_{i+1/2}$ is the diagonal matrix whose coefficients are the signs of the eigenvalues $\lambda _{i+1/2,1}$,…,$\lambda _{i+1/2,N}$. Besides,
In this context, the meaning of the well-balanced property of the reconstruction operator can be understood as follows: let us suppose, as in Remark 5.1, that an integral curve $\gamma$ of the linearly degenerate field can be described by a system of equations Equation5.3. Let us suppose that $\widetilde{W}(x)=(w_1(x), \dotsc , w_N(x), \sigma (x))$ is a stationary solution such that $\widetilde{W}(x) \in \gamma$ for all $x$, i.e.,
If we now apply a well-balanced reconstruction operator to the sequence $\{ \widetilde{W}(x_i) \}$, then the approximation functions $\widetilde{P}_i$ have to satisfy
In this section we consider numerical schemes of the form Equation3.6, in which the approximation functions used in the reconstruction operator are built by means of a WENO interpolation procedure using stencils with $r$ points; we denote this method simply as $r$-WENO, and the resulting scheme as $r$-WENO-Roe. For the details about WENO interpolation, see Reference23, Reference27, Reference33, Reference34. The reconstructions proposed in the $r$-WENO method are as follows:
These weights are calculated so that, on the one hand, the reconstruction operator is of order $2r - 1$ and, on the other hand, the weight $\omega _k$ is near to zero when the data on the stencil $S_k^r$ indicate the presence of a discontinuity.
In order to construct the approximation function at the cells, let us first define
Due to this fact, if a WENO-Roe scheme Equation3.6 is used to design a high order numerical method for a problem of the form Equation5.1, when the numerical scheme is written under the form Equation5.5, an extra term has to be added at the right-hand side:
Due to the definition of the reconstruction operator, the functions $P_i$ given by Equation6.1 or Equation6.2 provide only approximations of order $r$ at the interior points of the cells, while their derivatives give approximations of order $r-1$. Therefore, applying Theorem 3.2, the method Equation3.6 has only order $r$, while it has order $2r-1$ when it is applied to systems of conservation laws.
Remark 6.1
If, instead of a WENO method the $r$-ENO reconstruction operator is chosen, the expected order of the numerical scheme is $r$, since in this case the approximation functions coincide with interpolation polynomials constructed on the basis of stencils with $r$ points. Nevertheless, as commented in Reference34, the use of WENO approximations has several advantadges: it gives smoother operators, it is less sensible to round-off errors, and it avoids the use of conditionals in its practical implementation, being optimal for the vectorization of the algorithms.
It is however possible, performing some slight modifications on the WENO interpolation procedure, to obtain a method of order $2r-1$. The idea is as follows: instead of choosing the usual WENO reconstructions we consider the functions
where the weights now depend on $x$ and are calculated following the usual procedure in WENO reconstruction, so that the order of approximation is $2r-1$ in the cell. Unfortunately, the derivatives of these approximation functions are not easy to obtain. Instead, we substitute these derivatives by new WENO approximation functions
where, again, the weights are calculated, for every $x$, following the usual procedure in WENO reconstruction. Therefore, we again obtain order $2r-2$ in the cell.
Once these functions have been defined, we introduce the new approximation functions at the cells given either by
In practice, this integral is approached by means of a Gaussian quadrature of order at least $2r-1$. As a consequence, the weights $\omega ^\pm (x)$ and $\gamma ^\pm (x)$ have to be calculated only at the quadrature points.
Following the same steps as in the case of the $r$-WENO-Roe method, it can be easily shown that the resulting scheme (that will be denoted as modified$r$-WENO-Roe) is well balanced with order $2r-1$.
The computational cost of this modified numerical scheme is higher than those corresponding to standard WENO reconstructions, as two set of weights have to be calculated at every quadrature point. Moreover, the positivity of the weights is only ensured at the intercells, due to the choice of stencils. Therefore, in some cases negative weights may appear at interior quadrature points giving rise to oscillations and instabilities. For handling these negative weights, if necessary, the splitting technique of Shi, Hu and Shu (Reference32) can be applied. However, in some cases (see, e.g., Section 7.7) this technique does not completely remove the oscillations, and the scheme eventually crashes. The causes of this problem are currently under investigation.
We finish this section with a remark about time-stepping. As is usual in WENO interpolation based schemes, in order to obtain a full high resolution scheme it is necessary to use a high order method to advance in time. In the schemes considered here we have taken optimal high order TVD Runge-Kutta schemes (Reference19, Reference33).
6.1. Shallow-water equations with depth variations
The equations governing the flow of a shallow-water layer of fluid through a straight channel with constant rectangular cross-section can be written as
The variable $x$ makes reference to the axis of the channel and $t$ is time, $q(x,t)$ and $h(x,t)$ represent the mass-flow and the thickness, respectively, $g$ is gravity, and $H(x)$ is the depth function measured from a fixed level of reference. The fluid is supposed to be homogeneous and inviscid.
The system Equation6.3 can be rewritten under the form Equation5.1 with $N=2$,
Therefore, solutions corresponding to still water define straight lines in the $h$-$q$-$H$ space. As a consequence, Roe methods based on the family of segments are exactly well balanced for still-water solutions and well balanced with order 2 for general stationary solutions (see Reference2, Reference28).
The reconstruction operator proposed here to get higher order schemes is based on WENO reconstruction related to the variables $q$,$H$ and $\eta =h-H$ (this variable represents the water surface elevation). That is, given a sequence $(q_i,h_i,H_i)$ we consider the new sequence $(q_i, \eta _i, H_i)$ with $\eta _i=h_i-H_i$ and apply the $r$-WENO reconstruction operator to obtain polynomials
This reconstruction is exactly well balanced for stationary solutions corresponding to water at rest. In effect, if the sequence $(q_i, h_i, H_i)$ lie on the curve defined by Equation6.5, then $q_i=0$ and $\eta _i=C$. As a consequence,
and thus the reconstruction operator is well balanced (see Remark 5.2).
Applying Theorems 4.3 and 4.4, we deduce that the corresponding WENO-Roe schemes satisfy the $\mathcal{C}$-property, i.e., they are exactly well balanced for still-water solutions, and well balanced with order $r$ for general stationary solutions. To obtain a well-balanced numerical scheme with order $2r-1$, we have to add to the numerical scheme the modifications proposed in Section 6.
6.2. The two-layer shallow-water system
We now consider the equations of a one-dimensional flow of two superposed inmiscible layers of shallow-water fluids studied in Reference5: