Key words: High-order reconstructions (interpolations); A-WENO schemes; well-balanced schemes; equilibrium variables; local characteristic decomposition.

AMS subject classification: 76M20, 65M06, 35L65, 35L67.

1 Introduction↩︎

This paper is focused on the development of a new local characteristic decomposition (LCD) of equilibrium variables for hyperbolic systems of balance laws, which, in the two-dimensional (2-D) case, read as \[\boldsymbol{U}_t+\boldsymbol{F}(\boldsymbol{U})_x+\boldsymbol{G}(\boldsymbol{U})_y=B^x(\boldsymbol{U})\boldsymbol{U}_x+B^y(\boldsymbol{U})\boldsymbol{U}_y+\boldsymbol{S}(\boldsymbol{U}). \label{1461}\tag{1}\] Here, \(x\) and \(y\) are the spatial variables, \(t\) is time, \(\boldsymbol{U}\in\mathbb{R}^d\) is a vector of unknowns, \(\boldsymbol{F}\) and \(\boldsymbol{G}\) are the fluxes, \(\boldsymbol{S}\) is the source term, and \(B^x(\boldsymbol{U})\boldsymbol{U}_x\) and \(B^y(\boldsymbol{U})\boldsymbol{U}_y\) are nonconservative product terms.

Development of high-order numerical methods for (1 ) is a challenging task for the following two main reasons. First, solutions of (1 ) may develop discontinuities (even for infinitely smooth initial data) and thus high-order schemes should rely on non-oscillatory high-order reconstructions (interpolations) of the solutions out of the computed discrete quantities. To make these reconstructions (interpolations) non-oscillatory, one typically needs to use nonlinear limiting techniques often applied to the local characteristic variables via the LCD; see, e.g., [1]–[8]. The LCD is applied by computing local linearizations of the matrices \(\frac{\partial{\boldsymbol{F}}}{\partial{\boldsymbol{U}}}-B^x(\boldsymbol{U})\) and \(\frac{\partial{\boldsymbol{G}}}{\partial{\boldsymbol{U}}}-B^y(\boldsymbol{U})\), locally switching to the corresponding characteristic variables, performing the reconstruction (interpolation) to these local variables, and then switching back to the original variables \(\boldsymbol{U}\).

The use of characteristic-wise reconstructions has played an important role in the development of high-order WENO-type schemes for hyperbolic systems. In particular, LCD is often essential for reducing spurious oscillations near shocks, contact discontinuities, and sharp transition regions; see, e.g., [2], [3], [5]–[7] and references therein. Recently, other transformed-variable approaches, such as WENO reconstructions based on Riemann invariants, have also been investigated as alternatives to the standard LCD in certain systems; see, e.g., [9]. These works demonstrate that a suitable choice of variables for the nonlinear reconstruction is crucial for the robustness of high-order schemes.

Second, many (physically relevant) solutions of (1 ) are, in fact, small perturbations of certain steady states, which are supposed to be exactly preserved by good high-order schemes—this is a so-called well-balanced (WB) property. A large body of work has been devoted to the design of WB schemes for balance laws, including a variety of shallow water models, the Euler equations with gravitation, and other systems with sources and nonconservative product terms. For the Euler equations with gravitation, for example, high-order WB finite-volume (FV) and finite-difference (FD) WENO-type schemes have been developed; see, e.g., [10]–[15] and references therein. These methods are often based on local hydrostatic reconstructions, equilibrium-perturbation reconstructions, or carefully balanced discretizations of the flux and source terms. In particular, characteristic-wise WB WENO schemes have been developed for the Euler equations with gravitation in [15].

One of the ways to ensure that the high-order reconstruction (interpolation) preserves these steady states is to reconstruct (interpolate) the equilibrium variables instead of the conservative ones. We note that in many cases, the equilibrium variables can be obtained by rewriting the system (1 ) in the following form (see, e.g., [16]): \[\boldsymbol{U}_t+M^x(\boldsymbol{U})\boldsymbol{E}^x(\boldsymbol{U})_x+M^y(\boldsymbol{U})\boldsymbol{E}^y(\boldsymbol{U})_y=\boldsymbol{0}, \label{1464}\tag{2}\] where \[\begin{align} &M^x(\boldsymbol{U})\boldsymbol{E}^x(\boldsymbol{U})_x=\boldsymbol{F}(\boldsymbol{U})_x-B^x(\boldsymbol{U})\boldsymbol{U}_x-\boldsymbol{S}^x(\boldsymbol{U}),\\[1.ex] &M^y(\boldsymbol{U})\boldsymbol{E}^y(\boldsymbol{U})_y=\boldsymbol{G}(\boldsymbol{U})_y-B^y(\boldsymbol{U})\boldsymbol{U}_y-\boldsymbol{S}^y(\boldsymbol{U}), \end{align} \label{1464a}\tag{3}\] and \(\boldsymbol{S}^x(\boldsymbol{U})+\boldsymbol{S}^y(\boldsymbol{U})=\boldsymbol{S}(\boldsymbol{U})\). In (3 ), \(M^x,M^y\in\mathbb{R}^{d\times d}\) and \(\boldsymbol{E}^x(\boldsymbol{U}(x,y))\), \(\boldsymbol{E}^y(\boldsymbol{U}(x,y))\) are equilibrium variables, since \[\boldsymbol{E}^x(\boldsymbol{U}(x,y))=\boldsymbol{E}^x(y)\quadand\quad\boldsymbol{E}^y(\boldsymbol{U}(x,y))=\boldsymbol{E}^y(x) \label{1463a}\tag{4}\] at steady states satisfying \(\boldsymbol{E}^x(\boldsymbol{U})_x=\boldsymbol{E}^y(\boldsymbol{U})_y\equiv0\). Therefore, one may prefer to reconstruct \(\boldsymbol{E}^x\) in the \(x\)-direction and \(\boldsymbol{E}^y\) in the \(y\)-direction and then to recalculate the corresponding values of \(\boldsymbol{U}\) to ensure that (4 ) is satisfied at the discrete level.

The reconstruction of equilibrium variables and the use of LCD can be, in principle, viewed as two separate components of the numerical solution algorithm. However, when the reconstruction is performed for the equilibrium variables rather than for the conservative variables, the LCD based on the Jacobians of the conservative-variable formulation is no longer a natural characteristic decomposition for the quantities being interpolated. To combine the WB reconstruction of equilibrium variables with the oscillation-suppressing effect of the LCD, one needs to construct the characteristic decomposition associated with the equilibrium-variable formulation itself.

In this paper, we introduce a new LCD of equilibrium variables. To this end, we first rewrite the system (2 ) in the following two equivalent (for smooth solution) formulations: \[\begin{align} &\boldsymbol{E}^x(\boldsymbol{U})_t+C^x(\boldsymbol{U})\boldsymbol{E}^x(\boldsymbol{U})_x+D^x(\boldsymbol{U})\boldsymbol{U}_y=\tilde{\boldsymbol{I}}^x(\boldsymbol{U}),\\[1.ex] &\boldsymbol{E}^y(\boldsymbol{U})_t+D^y(\boldsymbol{U})\boldsymbol{U}_x+C^y(\boldsymbol{U})\boldsymbol{E}^y(\boldsymbol{U})_y=\tilde{\boldsymbol{I}}^y(\boldsymbol{U}), \end{align} \label{1465}\tag{5}\] where the matrices \(C^x\) and \(C^y\) and the source terms \(\tilde{\boldsymbol{I}}^x\) and \(\tilde{\boldsymbol{I}}^y\) are specified in §3. The proposed LCD is based on the matrices \(C^x\) and \(C^y\), which are associated with the equilibrium-variable formulations (5 ). We then compute the matrices \(C^x\) and \(C^y\) at the grid points and use them to compute the local characteristic equilibrium variables, which are reconstructed (interpolated) to obtain high-order values of \(\boldsymbol{E}^x\) and \(\boldsymbol{E}^y\), which, in turn, give us the corresponding high-order values of \(\boldsymbol{U}\) (solving nonlinear systems of equations may be required).

We implement the new LCD technique in the framework of flux globalization based WB alternative weighted essentially non-oscillatory (A-WENO) FD schemes recently introduced in [17]. The local characteristic equilibrium variables are interpolated using the fifth-order affine-invariant WENO-Z (Ai-WENO-Z) interpolations [18]–[20]. The developed A-WENO scheme is applied to five systems of balance laws including the nozzle flow system, the one- and two-layer shallow water equations, the compressible Euler equations with gravitation, and the 2-D Ripa system. We conduct several numerical experiments to demonstrate that the proposed LCD of equilibrium variables reduces spurious oscillations while preserving the WB property of the underlying flux-globalization based scheme.

2 Flux Globalization Based WB A-WENO Schemes: An Overview↩︎

In this section, we give an overview of the flux globalization based WB A-WENO schemes introduced in[17] for general nonconservative systems.

2.1 1-D Scheme↩︎

The one-dimensional (1-D) hyperbolic systems of balance laws \[\boldsymbol{U}_t+\boldsymbol{F}(\boldsymbol{U})_x=B(\boldsymbol{U})\boldsymbol{U}_x+\boldsymbol{S}(\boldsymbol{U}) \label{1464b}\tag{6}\] can be written in an equivalent quasi-conservative form: \[\boldsymbol{U}_t+\boldsymbol{K}(\boldsymbol{U})_x=\boldsymbol{0},\] where \(\boldsymbol{K}(\boldsymbol{U})\) is a global flux \[\boldsymbol{K}(\boldsymbol{U})=\boldsymbol{F}(\boldsymbol{U})-\boldsymbol{R}(\boldsymbol{U}),\quad\boldsymbol{R}(\boldsymbol{U})= \int\limits_{\hat{x}}^x\Big[B(\boldsymbol{U}(\xi,t))\boldsymbol{U}_\xi(\xi,t)+\boldsymbol{S}(\boldsymbol{U}(\xi,t))\Big]\,{\rm d}\xi,\] and \({\hat{x}}\) is an arbitrary number.

We first introduce a uniform mesh consisting of the cells \([x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}]\) of size \(x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}\equiv\Delta x\) centered at \(x_j=(x_{j-\frac{1}{2}}+x_{j+\frac{1}{2}})/2\), \(j=1,\ldots,N\). We assume that at a certain time level \(t\), the approximate solution, realized in terms of its cell centered values \(\boldsymbol{U}_j\approx\boldsymbol{U}(x_j,t)\), is available (in the rest of the paper, we will suppress the time-dependence of all of the indexed quantities for the sake of brevity). The solution is then evolved in time by solving the following system of ODEs: \[\frac{{\rm d}\boldsymbol{U}_j}{{\rm d}t}=-\frac{\boldsymbol{{\cal K}}_{j+\frac{1}{2}}-\boldsymbol{{\cal K}}_{j-\frac{1}{2}}}{\Delta x}, \label{2463}\tag{7}\] where \(\boldsymbol{{\cal K}}_{j+\frac{1}{2}}\) are the fifth-order A-WENO numerical fluxes (see [17], [21]): \[\boldsymbol{{\cal K}}_{j+\frac{1}{2}}=\boldsymbol{{\cal K}}^{\rm FV}_{j+\frac{1}{2}}-\frac{(\Delta x)^2}{24}(\boldsymbol{K}_{xx})_{j+\frac{1}{2}}+\frac{7(\Delta x)^4}{5760}(\boldsymbol{K}_{xxxx})_{j+\frac{1}{2}}.\] Here, \(\boldsymbol{{\cal K}}^{\rm FV}_{j+\frac{1}{2}}\) is a FV numerical flux (in the numerical experiments reported in §4, we have used the second-order WB path-conservative central-upwind numerical flux introduced in [22]), and \((\boldsymbol{K}_{xx})_{j+\frac{1}{2}}\) and \(({\boldsymbol{K}_{xxxx}})_{j+\frac{1}{2}}\) are the high-order correction terms. The numerical fluxes \(\boldsymbol{{\cal K}}^{\rm FV}_{j+\frac{1}{2}}=\boldsymbol{{\cal K}}^{\rm FV}_{j+\frac{1}{2}}(\boldsymbol{U}_{j+\frac{1}{2}}^\pm,\widehat\boldsymbol{U}_{j+\frac{1}{2}}^\pm)\) are computed using the one-sided interpolated values of \(\boldsymbol{U}\), and to enforce the WB evolution, one needs to use two copies of those values denoted by \(\boldsymbol{U}_{j+\frac{1}{2}}^\pm\) and \(\widehat\boldsymbol{U}_{j+\frac{1}{2}}^\pm\); see [17], [22] for details. The correction terms \((\boldsymbol{K}_{xx})_{j+\frac{1}{2}}\) and \(({\boldsymbol{K}_{xxxx}})_{j+\frac{1}{2}}\) are computed using the numerical fluxes \(\boldsymbol{{\cal K}}^{\rm FV}_{j+\frac{1}{2}}\), which have been already obtained: \[\begin{align} &(\boldsymbol{K}_{xx})_{j+\frac{1}{2}}=\frac{1}{12(\Delta x)^2}\Big[-\boldsymbol{{\cal K}}^{\rm FV}_{j-\frac{3}{2}}+16\boldsymbol{{\cal K}}^{\rm FV}_{j-\frac{1}{2}}- 30\boldsymbol{{\cal K}}^{\rm FV}_{j+\frac{1}{2}}+16\boldsymbol{{\cal K}}^{\rm FV}_{j+\frac{3}{2}}-\boldsymbol{{\cal K}}^{\rm FV}_{j+\frac{5}{2}}\Big],\\ &(\boldsymbol{K}_{xxxx})_{j+\frac{1}{2}}=\frac{1}{(\Delta x)^4}\Big[\boldsymbol{{\cal K}}^{\rm FV}_{j-\frac{3}{2}}-4\boldsymbol{{\cal K}}^{\rm FV}_{j-\frac{1}{2}}+6\boldsymbol{{\cal K}}^{\rm FV}_{j+\frac{1}{2}}- 4\boldsymbol{{\cal K}}^{\rm FV}_{j+\frac{3}{2}}+\boldsymbol{{\cal K}}^{\rm FV}_{j+\frac{5}{2}}\Big]; \end{align}\] see [21] for details. We would like to stress that these correction terms are needed to increase the order of the resulting scheme to the fifth order and that adding these terms typically does not cause oscillations as long as the reconstruction is performed in the local characteristic variables; see, e.g., [2], [4], [8], [23]–[25].

2.2 2-D Scheme↩︎

The 2-D hyperbolic systems of balance laws (1 ) can be similarly written in an equivalent quasi-conservative form: \[\boldsymbol{U}_t+\boldsymbol{K}(\boldsymbol{U})_x+\boldsymbol{L}(\boldsymbol{U})_y=\boldsymbol{0},\] where \(\boldsymbol{K}(\boldsymbol{U})\) and \(\boldsymbol{L}(\boldsymbol{U})\) are the global fluxes \[\begin{align} &\boldsymbol{K}(\boldsymbol{U})=\boldsymbol{F}(\boldsymbol{U})-\boldsymbol{R}^x(\boldsymbol{U}),&&\boldsymbol{R}^x(\boldsymbol{U})= \int\limits_{\hat{x}}^x\Big[B^x(\boldsymbol{U}(\xi,y,t))\boldsymbol{U}_\xi(\xi,y,t)+\boldsymbol{S}^x(\boldsymbol{U}(\xi,y,t))\Big]\,{\rm d}\xi,\\ &\boldsymbol{L}(\boldsymbol{U})=\boldsymbol{G}(\boldsymbol{U})-\boldsymbol{R}^y(\boldsymbol{U}),&&\boldsymbol{R}^y(\boldsymbol{U})= \int\limits_{\hat{y}}^y\Big[B^y(\boldsymbol{U}(x,\eta,t))\boldsymbol{U}_\eta(x,\eta,t)+\boldsymbol{S}^y(\boldsymbol{U}(x,\eta,t))\Big]\,{\rm d}\eta, \end{align}\] and \({\hat{x}}\) and \({\hat{y}}\) are arbitrary numbers.

Let \([x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}]\times[y_{k-\frac{1}{2}},y_{k+\frac{1}{2}}]\), \(j=1,\cdots,N_x\), \(k=1,\cdots,N_y\) be the uniform 2-D cells centered at \((x_j,y_k)\) with \(x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}\equiv\Delta x\), \(y_{k+\frac{1}{2}}-y_{k-\frac{1}{2}}\equiv\Delta y\), \(x_j=(x_{j-\frac{1}{2}}+x_{j+\frac{1}{2}})/2\), and \(y_k=(y_{k-\frac{1}{2}}+y_{k+\frac{1}{2}})/2\). We denote by \(\boldsymbol{U}_{j,k}\approx\boldsymbol{U}(x_j,y_k,t)\) the computed cell centered values, which are assumed to be available at a certain time level \(t\). The point values \(\boldsymbol{U}_{j,k}\) are evolved in time by solving the following system of ODEs: \[\frac{{\rm d}\boldsymbol{U}_{j,k}}{{\rm d}t}=-\frac{\boldsymbol{{\cal K}}_{{j+\frac{1}{2}},k}-\boldsymbol{{\cal K}}_{{j-\frac{1}{2}},k}}{\Delta x}- \frac{\boldsymbol{{\cal L}}_{j,{k+\frac{1}{2}}}-\boldsymbol{{\cal L}}_{j,{k-\frac{1}{2}}}}{\Delta y}, \label{2468a}\tag{8}\] where \(\boldsymbol{{\cal K}}_{{j+\frac{1}{2}},k}\) and \(\boldsymbol{{\cal L}}_{j,{k+\frac{1}{2}}}\) are the fifth-order A-WENO numerical fluxes (see [17], [21]): \[\begin{align} &\boldsymbol{{\cal K}}_{{j+\frac{1}{2}},k}= \boldsymbol{{\cal K}}^{\rm FV}_{{j+\frac{1}{2}},k}-\frac{(\Delta x)^2}{24}(\boldsymbol{K}_{xx})_{{j+\frac{1}{2}},k}+\frac{7(\Delta x)^4}{5760}(\boldsymbol{K}_{xxxx})_{{j+\frac{1}{2}},k},\\ &\boldsymbol{{\cal L}}_{j,{k+\frac{1}{2}}}= \boldsymbol{{\cal L}}^{\rm FV}_{j,{k+\frac{1}{2}}}-\frac{(\Delta y)^2}{24}(\boldsymbol{L}_{yy})_{j,{k+\frac{1}{2}}}+\frac{7(\Delta y)^4}{5760}(\boldsymbol{L}_{yyyy})_{j,{k+\frac{1}{2}}}. \end{align}\] Here, \(\boldsymbol{{\cal K}}^{\rm FV}_{{j+\frac{1}{2}},k}\) and \(\boldsymbol{{\cal L}}^{\rm FV}_{j,{k+\frac{1}{2}}}\) are finite-volume numerical fluxes and \((\boldsymbol{K}_{xx})_{{j+\frac{1}{2}},k}\), \(({\boldsymbol{K}_{xxxx}})_{{j+\frac{1}{2}},k}\), \((\boldsymbol{L}_{yy})_{j,{k+\frac{1}{2}}}\), and \(({\boldsymbol{L}_{yyyy}})_{j,{k+\frac{1}{2}}}\) are the high-order correction terms. The numerical fluxes \(\boldsymbol{{\cal K}}^{\rm FV}_{{j+\frac{1}{2}},k}=\boldsymbol{{\cal K}}^{\rm FV}_{{j+\frac{1}{2}},k}\big(\boldsymbol{U}_{{j+\frac{1}{2}},k}^\pm\big)\) and \(\boldsymbol{{\cal L}}^{\rm FV}_{j,{k+\frac{1}{2}}}=\boldsymbol{{\cal L}}^{\rm FV}_{j,{k+\frac{1}{2}}}\big(\boldsymbol{U}_{j,{k+\frac{1}{2}}}^\pm\big)\) are computed using the one-sided interpolated values of \(\boldsymbol{U}\), and to enforce the WB evolution, one may need to use the numerical diffusion switch functions; see [26] for an example of such switch function used for the compressible Euler equations with gravitation. The correction terms \((\boldsymbol{K}_{xx})_{{j+\frac{1}{2}},k}\), \(({\boldsymbol{K}_{xxxx}})_{{j+\frac{1}{2}},k}\), \((\boldsymbol{L}_{yy})_{j,{k+\frac{1}{2}}}\), and \(({\boldsymbol{L}_{yyyy}})_{j,{k+\frac{1}{2}}}\) are computed using the numerical fluxes \(\boldsymbol{{\cal K}}^{\rm FV}_{{j+\frac{1}{2}},k}\) and \(\boldsymbol{{\cal L}}^{\rm FV}_{{j+\frac{1}{2}},k}\), which have been already obtained: \[\begin{align} &(\boldsymbol{K}_{xx})_{{j+\frac{1}{2}},k}=\frac{1}{12(\Delta x)^2}\Big[-\boldsymbol{{\cal K}}^{\rm FV}_{j-\frac{3}{2},k}+16\boldsymbol{{\cal K}}^{\rm FV}_{{j-\frac{1}{2}},k}- 30\boldsymbol{{\cal K}}^{\rm FV}_{{j+\frac{1}{2}},k}+16\boldsymbol{{\cal K}}^{\rm FV}_{j+\frac{3}{2},k}-\boldsymbol{{\cal K}}^{\rm FV}_{j+\frac{5}{2},k}\Big],\\ &(\boldsymbol{K}_{xxxx})_{{j+\frac{1}{2}},k}=\frac{1}{(\Delta x)^4}\Big[\boldsymbol{{\cal K}}^{\rm FV}_{j-\frac{3}{2},k}-4\boldsymbol{{\cal K}}^{\rm FV}_{{j-\frac{1}{2}},k}+ 6\boldsymbol{{\cal K}}^{\rm FV}_{{j+\frac{1}{2}},k}-4\boldsymbol{{\cal K}}^{\rm FV}_{j+\frac{3}{2},k}+\boldsymbol{{\cal K}}^{\rm FV}_{j+\frac{5}{2},k}\Big];\\ &(\boldsymbol{L}_{yy})_{j,{k+\frac{1}{2}}}=\frac{1}{12(\Delta y)^2}\Big[-\boldsymbol{{\cal L}}^{\rm FV}_{j,k-\frac{3}{2}}+16\boldsymbol{{\cal L}}^{\rm FV}_{j,{k-\frac{1}{2}}}- 30\boldsymbol{{\cal L}}^{\rm FV}_{j,{k+\frac{1}{2}}}+16\boldsymbol{{\cal L}}^{\rm FV}_{j,k+\frac{3}{2}}-\boldsymbol{{\cal L}}^{\rm FV}_{j,k+\frac{5}{2}}\Big],\\ &(\boldsymbol{L}_{yyyy})_{j,{k+\frac{1}{2}}}=\frac{1}{(\Delta y)^4}\Big[\boldsymbol{{\cal L}}^{\rm FV}_{j,k-\frac{3}{2}}-4\boldsymbol{{\cal L}}^{\rm FV}_{j,{k-\frac{1}{2}}}+ 6\boldsymbol{{\cal L}}^{\rm FV}_{j,{k+\frac{1}{2}}}- 4\boldsymbol{{\cal L}}^{\rm FV}_{j,k+\frac{3}{2}}+\boldsymbol{{\cal L}}^{\rm FV}_{j,k+\frac{5}{2}}\Big]; \end{align}\] see [21] for details.

3 Local Characteristic Decomposition of Equilibrium Variables↩︎

In this section, we apply the LCD approach introduced in [1]–[8] to the Ai-WENO-Z interpolation of the equilibrium variables.

Throughout this section, we use the following notation for the equilibrium variables: \[\label{3461a} \boldsymbol{E}^\nu(\boldsymbol{U})=\boldsymbol{W}^\nu(\boldsymbol{U})+\boldsymbol{I}^\nu(\boldsymbol{U}),\tag{9}\] where \(\nu\) stands for either \(x\) or \(y\), \(\boldsymbol{W}^\nu\) denotes the local part of the equilibrium variable \(\boldsymbol{E}^\nu\), and \(\boldsymbol{I}^\nu\) is its global (integral) part. When the equilibrium variables are local, \(\boldsymbol{I}^\nu\equiv0\), and thus \(\boldsymbol{E}^\nu=\boldsymbol{W}^\nu\). In the 1-D case, the directional superscript \(\nu\) is omitted so that we write \[\boldsymbol{E}(\boldsymbol{U})=\boldsymbol{W}(\boldsymbol{U})+\boldsymbol{I}(\boldsymbol{U}). \label{3461b}\tag{10}\]

3.1 1-D Case↩︎

We first rewrite the system (6 ) in the following equivalent (for smooth solution) form: \[\boldsymbol{U}_t+M(\boldsymbol{U})\boldsymbol{E}(\boldsymbol{U})_x=\boldsymbol{0}, \label{1465b}\tag{11}\] where \[M(\boldsymbol{U})\boldsymbol{E}(\boldsymbol{U})_x=\boldsymbol{F}(\boldsymbol{U})_x-B(\boldsymbol{U})\boldsymbol{U}_x-\boldsymbol{S}(\boldsymbol{U}), \label{1464c}\tag{12}\] which vanishes at steady states. In (12 ), \(M\in\mathbb{R}^{d\times d}\) and \(\boldsymbol{E}\) is the vector of equilibrium variables, which are constant at steady states. We then rewrite (11 ) again as \[\boldsymbol{E}(\boldsymbol{U})_t+C(\boldsymbol{U})\boldsymbol{E}(\boldsymbol{U})_x=\tilde{\boldsymbol{I}}(\boldsymbol{U}). \label{1465a}\tag{13}\]

We emphasize that (13 ) is not used as a time-evolution equation in the numerical method; it is introduced only to identify the matrix \(C(\boldsymbol{U})\) used in the LCD of equilibrium variables.

We now need to specify \(C\) and \(\tilde{\boldsymbol{I}}\) in (13 ). To this end, we consider two possible cases.

\(\bullet\) If \(\boldsymbol{E}(\boldsymbol{U})\) is a local equilibrium variable, that is \(\boldsymbol{I}(\boldsymbol{U})\equiv\boldsymbol{0}\) in (10 ), then we simply multiply (12 ) by \(\frac{\partial{\boldsymbol{E}}}{\partial{\boldsymbol{U}}}\) and obtain (13 ) with \(C(\boldsymbol{U})=\frac{\partial{\boldsymbol{E}}}{\partial{\boldsymbol{U}}}M(\boldsymbol{U})\) and \(\tilde{\boldsymbol{I}}(\boldsymbol{U})\equiv\boldsymbol{0}\).

\(\bullet\) If \(\boldsymbol{E}(\boldsymbol{U})\) is a global equilibrium variable, that is, \(\boldsymbol{I}(\boldsymbol{U})\ne\boldsymbol{0}\) and \[\boldsymbol{I}(\boldsymbol{U})=\int\limits^x\big[\boldsymbol{H}(\boldsymbol{U})+N(\boldsymbol{U})\boldsymbol{U}_x\big]\,{\rm d}x, \label{3461c}\tag{14}\] where \(\boldsymbol{H}(\boldsymbol{U})\) is a given vector function and \(N(\boldsymbol{U})\) is a given matrix, then we proceed in a different way. We first differentiate (10 ) and (14 ) with respect to \(t\) to obtain \[\boldsymbol{E}(\boldsymbol{U})_t+\frac{\partial\boldsymbol{W}}{\partial\boldsymbol{U}}(\boldsymbol{U})\,\boldsymbol{U}_t+\int\limits^x\bigg[\frac{\partial\boldsymbol{H}}{\partial\boldsymbol{U}}(\boldsymbol{U})\,\boldsymbol{U}_t+N(\boldsymbol{U})_t\,\boldsymbol{U}_x+ N(\boldsymbol{U})\boldsymbol{U}_{xt}\bigg]\,{\rm d}x. \label{3461d}\tag{15}\] We then use (11 ) and the identity \[N(\boldsymbol{U})\boldsymbol{U}_{xt}=(N(\boldsymbol{U})\boldsymbol{U}_t)_x-N(\boldsymbol{U})_x\,\boldsymbol{U}_t\] to rewrite (15 ) as \[\begin{align} \boldsymbol{E}(\boldsymbol{U})_t&+\bigg[\frac{\partial\boldsymbol{W}}{\partial\boldsymbol{U}}(\boldsymbol{U})+N(\boldsymbol{U})\bigg]M(\boldsymbol{U})\boldsymbol{E}(\boldsymbol{U})_x\\ =&-\int\limits^x\bigg[\frac{\partial\boldsymbol{H}}{\partial\boldsymbol{U}}(\boldsymbol{U})M(\boldsymbol{U})\boldsymbol{E}(\boldsymbol{U})_x-\bigg(\frac{\partial\boldsymbol{N}_1}{\partial\boldsymbol{U}}(\boldsymbol{U})\,\boldsymbol{U}_x\dots \frac{\partial\boldsymbol{N}_d}{\partial\boldsymbol{U}}(\boldsymbol{U})\,\boldsymbol{U}_x\bigg)M(\boldsymbol{U})\boldsymbol{E}(\boldsymbol{U})_x\\[0.5ex] &+\bigg(\frac{\partial\boldsymbol{N}_1}{\partial\boldsymbol{U}}(\boldsymbol{U})M(\boldsymbol{U})\boldsymbol{E}(\boldsymbol{U})_x\dots\frac{\partial\boldsymbol{N}_d}{\partial\boldsymbol{U}}(\boldsymbol{U})M(\boldsymbol{U})\boldsymbol{E}(\boldsymbol{U})_x \bigg)\bigg]\,{\rm d}x. \end{align} \label{3461e}\tag{16}\] Here, we have used the following notation: \(N=(\boldsymbol{N}_1\dots\boldsymbol{N}_d)\), where \(\boldsymbol{N}_i\) is the \(i^{\rm th}\) column of the matrix \(N\). Notice that (16 ) corresponds to (13 ) with \[C(\boldsymbol{U})=\bigg[\frac{\partial \boldsymbol{W}(\boldsymbol{U}) }{\partial \boldsymbol{U}}(\boldsymbol{U})+N(\boldsymbol{U}) \bigg]M(\boldsymbol{U})\] and \(\tilde{\boldsymbol{I}}(\boldsymbol{U})\) given by the right-hand side of (16 ).

In fact, the details of \(\tilde{\boldsymbol{I}}\) are not important as this source term does not influence the LCD of equilibrium variables, which is based on the matrix \(C\) only.

We then evaluate the matrix \(C\) at the grid points to obtain the constant matrices \(C_j:=C(\boldsymbol{U}_j)\), which can be diagonalized using the matrices \(Q_j\) and \(Q_j^{-1}\) to obtain \(\Lambda_j=Q_j^{-1}C_jQ_j\), where \(\Lambda_j\) is a diagonal matrix containing the eigenvalues of \(C_j\).

Next, we introduce the local characteristic equilibrium variables in the neighborhood of \(x=x_j\): \[\boldsymbol{\Gamma}_\ell=Q_j^{-1}\boldsymbol{E}_\ell,\quad\ell=j\pm2,j\pm1,j, \label{3463}\tag{17}\] apply the fifth-order Ai-WENO-Z (or any other fifth-order WENO-type) interpolation to evaluate the values \(\boldsymbol{\Gamma}^+_{j-\frac{1}{2}}\) and \(\boldsymbol{\Gamma}^-_{j+\frac{1}{2}}\), and finally obtain \[\boldsymbol{E}^\pm_{j\mp{\frac{1}{2}}}=Q_j\boldsymbol{\Gamma}^\pm_{j\mp{\frac{1}{2}}}. \label{3464}\tag{18}\] Equipped with these values, we proceed as in [17], [22] and solve the nonlinear equations (see [27]) to recover the values \(\boldsymbol{U}^\pm_{j\mp{\frac{1}{2}}}\) and \(\widehat\boldsymbol{U}^\pm_{j\mp{\frac{1}{2}}}\) needed to evaluate the numerical fluxes \(\boldsymbol{{\cal K}}^{\rm FV}_{j+\frac{1}{2}}\).

3.2 2-D Case↩︎

In the 2-D case, the LCD of the equilibrium variables is based on the formulation of the 2-D system (1 ) given in (5 ). In fact, to perform the LCD, we will only need to specify \(C^x\) and \(C^y\) since the terms with \(D^x\), \(D^y\), \(\tilde{\boldsymbol{I}}^x\), and \(\tilde{\boldsymbol{I}}^y\) do not influence the LCD process. The computation of \(C^x\) and \(C^y\) is similar to the computation of the matrix \(C\) in §3.1, and they are given by \[C^\nu(\boldsymbol{U})=\bigg[\frac{\partial\boldsymbol{W}^\nu}{\partial\boldsymbol{U}}(\boldsymbol{U})+N^\nu(\boldsymbol{U})\bigg]M^\nu(\boldsymbol{U}),\quad\nu\in\{x,y\}, \label{3464aa}\tag{19}\] where the matrices \(N^\nu(\boldsymbol{U})\) are from the expressions of the equilibrium variables (9 ) with \[\boldsymbol{I}^x=\int\limits^x\bigg[\boldsymbol{H}^x(\boldsymbol{U})+N^x(\boldsymbol{U})\,\boldsymbol{U}_x \bigg]\,{\rm d}x,\quad \boldsymbol{I}^y=\int\limits^y\bigg[\boldsymbol{H}^y(\boldsymbol{U})+N^y(\boldsymbol{U})\,\boldsymbol{U}_y \bigg]\,{\rm d}y.\] Here, \(\boldsymbol{H}^x(\boldsymbol{U})\) and \(\boldsymbol{H}^y(\boldsymbol{U})\) are given vector functions, and \(N^x(\boldsymbol{U})\) and \(N^y(\boldsymbol{U})\) are given matrices. In particular, if the equilibrium variables are local, then (19 ) reduces to \[C^\nu(\boldsymbol{U})=\frac{\partial\boldsymbol{E}^\nu}{\partial\boldsymbol{U}}(\boldsymbol{U})M^\nu(\boldsymbol{U}),\quad\nu\in\{x,y\}.\]

We then evaluate the matrices \(C^x\) and \(C^y\) at the grid points to obtain the constant matrices \(C^x_{j,k}:=C^x(\boldsymbol{U}_{j,k})\) and \(C^y_{j,k}:=C^y(\boldsymbol{U}_{j,k})\), which can be diagonalized using the matrices \(Q^x_{j,k}\) and \(Q^y_{j,k}\) to obtain \(\Lambda^x_{j,k}=(Q^x_{j,k})^{-1}C^x_{j,k}Q^x_{j,k}\) and \(\Lambda^y_{j,k}=(Q^y_{j,k})^{-1}C^y_{j,k}Q^y_{j,k}\), where \(\Lambda^x_{j,k}\) and \(\Lambda^y_{j,k}\) are the diagonal matrices containing the eigenvalues of \(C^x_{j,k}\) and \(C^y_{j,k}\), respectively.

Next, we introduce the local characteristic equilibrium variables in the neighborhood of \((x,y)=(x_j,y_k)\): \[\boldsymbol{\Gamma}^x_{\ell,k}=(Q^x_{j,k})^{-1}\boldsymbol{E}^x_{\ell,k},\quad\ell=j\pm2,j\pm1,j,\quadand\quad \boldsymbol{\Gamma}^y_{j,\ell}=(Q^y_{j,k})^{-1}\boldsymbol{E}^y_{j,\ell},\quad\ell=k\pm2,k\pm1,k, \label{3463a}\tag{20}\] apply the fifth-order Ai-WENO-Z (or any other fifth-order WENO-type) interpolation in the \(x\)- and \(y\)-directions to evaluate the values \((\boldsymbol{\Gamma}^x_{j\mp{\frac{1}{2}},k})^\pm\) and \((\boldsymbol{\Gamma}^y_{j,k\mp{\frac{1}{2}}})^\pm\), respectively, and finally obtain \[\big(\boldsymbol{E}^x_{j\mp{\frac{1}{2}},k}\big)^\pm=Q^x_{j,k}\big(\boldsymbol{\Gamma}^x_{j\mp{\frac{1}{2}},k}\big)^\pm,\quad \big(\boldsymbol{E}^y_{j,k\mp{\frac{1}{2}}}\big)^\pm=Q^y_{j,k}\big(\boldsymbol{\Gamma}^y_{j,k\mp{\frac{1}{2}}}\big)^\pm. \label{3464a}\tag{21}\]

Next, we show the application of the introduced LCD of equilibrium variables to five particular systems of balance laws.

3.3 Application to the Nozzle Flow System↩︎

In this section, we consider the 1-D nozzle flow system, which reads as (6 ) with \[\boldsymbol{U}=(\sigma\rho,\sigma\rho u)^\top,\quad\boldsymbol{F}(\boldsymbol{U})=(\sigma\rho u,\sigma\rho u^2+\sigma p)^\top,\quad B(\boldsymbol{U})=0,\quad\boldsymbol{S}(\boldsymbol{U})=(0,p\sigma_x)^\top,\] where \(\rho\) is the density, \(u\) is the velocity, \(p(\rho)=\kappa\rho^\gamma\) is the pressure, \(\kappa>0\) and \(1<\gamma<\frac{5}{3}\) are constants, and \(\sigma=\sigma(x)\) denotes the cross-section of the nozzle. The studied nozzle flow system admits steady-state solutions satisfying \(M(\boldsymbol{U})\boldsymbol{E}(\boldsymbol{U})_x=\boldsymbol{0}\) with \[\begin{align} &M(\boldsymbol{U})=\begin{pmatrix}1&0\\u&\sigma\rho\end{pmatrix},\quad\boldsymbol{E}(\boldsymbol{U})=\begin{pmatrix}q\\{\cal E}\end{pmatrix},\quad q=\sigma\rho u,\quad{\cal E}=\frac{u^2}{2}+\frac{\kappa\gamma}{\gamma-1}\rho^{\gamma-1}. \end{align}\]

In order to apply the fifth-order Ai-WENO-Z interpolation to the equilibrium variables, we first compute \[{\cal E}_j=\frac{u_j^2}{2}+\frac{\kappa\gamma}{\gamma-1}(\rho_j)^{\gamma-1},\] where \(u_j=q_j/(\sigma\rho)_j\), \(\rho_j=(\sigma\rho)_j/\sigma_j\), and \(\sigma_j=\sigma(x_j)\), and then evaluate the matrices \[\begin{align} &C_j=\begin{pmatrix}u_j&(\sigma\rho)_j\\\dfrac{\kappa\gamma(\rho_j)^{\gamma-1}}{(\sigma\rho)_j}&u_j\end{pmatrix},\quad Q_j=\begin{pmatrix}(\sigma\rho)_j&(\sigma\rho)_j\\ -\sqrt{\kappa\gamma}(\rho_j)^{\frac{\gamma-1}{2}}&\sqrt{\kappa\gamma}(\rho_j)^{\frac{\gamma-1}{2}}\end{pmatrix},\\ &Q_j^{-1}=\frac{1}{2\sqrt{\kappa\gamma}(\sigma\rho)_j(\rho_j)^{\frac{\gamma-1}{2}}} \begin{pmatrix}\sqrt{\kappa\gamma}(\rho_j)^{\frac{\gamma-1}{2}}&-(\sigma\rho)_j\\ \sqrt{\kappa\gamma}(\rho_j)^{\frac{\gamma-1}{2}}&(\sigma\rho)_j\end{pmatrix}. \end{align}\] We now implement the LCD of the equilibrium variables followed by the fifth-order Ai-WENO-Z interpolation giving \(\boldsymbol{\Gamma}^\pm_{j\mp{\frac{1}{2}}}\) and then \(q^\pm_{j\mp{\frac{1}{2}}}\) and \({\cal E}^\pm_{j\mp{\frac{1}{2}}}\). After that, we apply the same fifth-order Ai-WENO-Z interpolation to obtain the one-sided values of the cross-section of the nozzle \(\sigma^\pm_{j\mp{\frac{1}{2}}}\), and then solve the nonlinear equations as described in [17] to obtain \((\sigma\rho)^\pm_{j\mp{\frac{1}{2}}}\) and \((\widehat{\sigma\rho})^\pm_{j\mp{\frac{1}{2}}}\).

3.4 Application to the Saint-Venant System with Manning Friction↩︎

In this section, we consider the 1-D Saint-Venant system of shallow water equations with Manning friction, which reads as (6 ) with \[\boldsymbol{U}=(h,q)^\top,\quad\boldsymbol{F}(\boldsymbol{U})=\Big(q,hu^2+{\frac{1}{2}}gh^2\Big)^\top,\quad B(\boldsymbol{U})=0,\quad\boldsymbol{S}(\boldsymbol{U})=(0,-ghZ_x-ghS_f)^\top,\] where \(h\) is the water depth, \(u\) is the velocity, \(q=hu\) represents the discharge, \(Z(x)\) is a function describing the bottom topography, which can be discontinuous, \(g\) is the constant acceleration due to gravity, \(S_f\) is the Manning friction term (see, e.g., [29]) given by \(S_f=n^2q|q|h^{-\frac{10}{3}}\). The studied Saint-Venant system admits steady-state solutions satisfying \(M(\boldsymbol{U})\boldsymbol{E}(\boldsymbol{U})_x=\boldsymbol{0}\) with \[\begin{align} M(\boldsymbol{U})=\begin{pmatrix}1&0\\u&h\end{pmatrix},\quad\boldsymbol{E}(\boldsymbol{U})=\begin{pmatrix}q\\{\cal E}\end{pmatrix},\quad {\cal E}=\frac{u^2}{2}+g(h+Z)+\int\limits_{\hat{x}}^x{gS_f}\,{\rm d}\xi. \end{align}\]

In order to apply the fifth-order Ai-WENO-Z interpolation to the equilibrium variables, we first compute \[{\cal E}_j=\frac{u_j^2}{2}+g(h_j+Z_j)+{\cal I}_j,\] where \(u_j=q_j/h_j\), \(Z_j=Z(x_j)\), and \({\cal I}_j\) is a fifth-order approximation of the integral \(\int_{x_{-\frac{5}{2}}}^{x_j}gS_f\,{\rm d}x\), in which we have set \(\hat{x}=x_{-\frac{5}{2}}\). The values \({\cal I}_j\) are computed as follows. First, we evaluate \(f_j:=g(S_f)_j={gn^2q_j|q_j|}{h_j^{-\frac{10}{3}}}\) for \(j=1,\ldots,N\) and extend these values to \(j=-5,\ldots,0\) and \(j=N+1,\ldots,N+5\) using the prescribed boundary conditions (for \(h\) and \(q\)) implemented within the ghost cell framework. We then construct the interpolating polynomial for \(f\) using the five points \((x_{-\frac{5}{2}},f_{-\frac{5}{2}}^+)\), \((x_{-\frac{9}{4}},f_{-\frac{9}{4}})\), \((x_{-2},f_{-2})\), \((x_{-\frac{7}{4}},f_{-\frac{7}{4}})\), and \((x_{-\frac{3}{2}},f_{-\frac{3}{2}}^-)\) and integrate it over the interval \([x_{-\frac{5}{2}},x_{-2}]\) to obtain \[{\cal I}_{-2}=\frac{\Delta x}{360}\left[29f_{-\frac{5}{2}}^++124f_{-\frac{9}{4}}+24f_{-2}+4f_{-\frac{7}{4}}-f_{-\frac{3}{2}}^-\right], \label{3465a}\tag{22}\] where \(f_{j\pm\frac{1}{4}}\) are evaluated as in [28], but now we use the Ai-WENO-Z interpolant instead of the WENO-Z one implemented in [28].

We then proceed recursively: construct the interpolating polynomials for \(f\) using \((x_{j-1},f_{j-1})\), \((x_{j-\frac{3}{4}},f_{j-\frac{3}{4}})\), \((x_{j-\frac{1}{2}},f_{j-\frac{1}{2}})\), \((x_{j-\frac{1}{4}},f_{j-\frac{1}{4}})\), and \((x_j,f_j)\), integrate them over the corresponding interval \([x_{j-1},x_j]\), and end up with \[{\cal I}_j={\cal I}_{j-1}+\frac{\Delta x}{90}\left[7f_{j-1}+32f_{j-\frac{3}{4}}+12f_{j-\frac{1}{2}}+32f_{j-\frac{1}{4}}+7f_{j}\right],\quad j=1,\ldots,N+3, \label{3465b}\tag{23}\] where the point values \(f_{j-\frac{3}{4}}\) and \(f_{j-\frac{1}{4}}\) are computed as in [28], but using the Ai-WENO-Z interpolant and \(f_{{j-\frac{1}{2}}}:=\big(f_{j-\frac{1}{2}}^-+f_{j-\frac{1}{2}}^+\big)/2\).

We then evaluate the matrices \[\begin{align} C_j=\begin{pmatrix}u_j&h_j\\g&u_j\end{pmatrix},\quad Q_j=\begin{pmatrix}\sqrt{h_j}&\sqrt{h_j}\\-\sqrt{g}&\sqrt{g}\end{pmatrix},\quad Q_j^{-1}=\frac{1}{2\sqrt{gh_j}}\begin{pmatrix}\sqrt{g}&-\sqrt{h_j}\\\sqrt{g}&\sqrt{h_j}\end{pmatrix}, \end{align}\] and implement the LCD of the equilibrium variables followed by the fifth-order Ai-WENO-Z interpolation giving \(\boldsymbol{\Gamma}^\pm_{j\mp{\frac{1}{2}}}\) and then \(q^\pm_{j\mp{\frac{1}{2}}}\) and \({\cal E}^\pm_{j\mp{\frac{1}{2}}}\). After that, we apply the same fifth-order Ai-WENO-Z interpolation to obtain the one-sided values of the bottom topography \(Z^\pm_{j\mp{\frac{1}{2}}}\), and then solve the nonlinear equations \[{\cal E}^\pm_{j+\frac{1}{2}}=\frac{\big(q^\pm_{j+\frac{1}{2}}\big)^2}{2\big(h^\pm_{j+\frac{1}{2}}\big )^2}+g\left(h^\pm_{j+\frac{1}{2}}+Z^\pm_{j+\frac{1}{2}}\right)+{\cal I}_{j+\frac{1}{2}}, \label{3461}\tag{24}\] and \[{\cal E}^\pm_{j+\frac{1}{2}}=\frac{\big(q^\pm_{j+\frac{1}{2}}\big)^2}{2\big(\,\widehat h^\pm_{j+\frac{1}{2}}\big)^2}+g\left(\widehat h^\pm_{j+\frac{1}{2}}+Z_{j+\frac{1}{2}}\right)+ {\cal I}_{j+\frac{1}{2}},\quad Z_{j+\frac{1}{2}}={\frac{1}{2}}\left(Z^+_{j+\frac{1}{2}}+Z^-_{j+\frac{1}{2}}\right) \label{3462}\tag{25}\] to obtain \(h^\pm_{j+\frac{1}{2}}\) and \(\,\widehat h^\pm_{j+\frac{1}{2}}\), respectively. Here, the integrals \({\cal I}_{j+\frac{1}{2}}\) are evaluated recursively by the fifth-order quadrature: we first set \({\cal I}_{-\frac{5}{2}}=0\) and then compute \[{\cal I}_{j+\frac{1}{2}}={\cal I}_{j-\frac{1}{2}}+\frac{\Delta x}{90}\left[7f_{j-\frac{1}{2}}^++32f_{j-\frac{1}{4}}+12f_j+32f_{j+\frac{1}{4}}+7f_{j+\frac{1}{2}}^-\right] \label{3465c}\tag{26}\] for \(j=-2,\ldots,N+2\), where we have used the quadrature, which is obtained by constructing the interpolating polynomials for \(f\) using the five points \((x_{j-\frac{1}{2}},f_{j-\frac{1}{2}}^+)\), \((x_{j-\frac{1}{4}},f_{j-\frac{1}{4}})\), \((x_j,f_j)\), \((x_{j+\frac{1}{4}},f_{j+\frac{1}{4}})\), and \((x_{j+\frac{1}{2}},f_{j+\frac{1}{2}}^-)\) and integrating them over the corresponding intervals \([x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}]\). As before, the values \(f_{j\pm\frac{1}{4}}\) are computed as in [28], but using the Ai-WENO-Z interpolant. Notice that equations (24 ) and (25 ) are cubic and we solve them exactly as described in [27].

3.5 Application to the Two-Layer Shallow Water System↩︎

In this section, we consider the 1-D two-layer shallow water system, which reads as (6 ) with \[\begin{align} &\boldsymbol{U}=(h_1,q_1,h_2,q_2)^\top,\\ &\boldsymbol{F}(\boldsymbol{U})=\big(q_1,h_1u_1^2+\dfrac{g}{2}h_1^2,q_2,h_2u_2^2+\dfrac{g}{2}h_2^2\big)^\top,\\ &\boldsymbol{S}(\boldsymbol{U})=(0,-gh_1Z_x,0,-gh_2Z_x)^\top, \end{align} \qquad B(\boldsymbol{U})=\begin{pmatrix}0&0&0&0\\0&0&-gh_1&0\\0&0&0&0\\-rgh_2&0&0&0\end{pmatrix}.\] Here, \(h_1\) and \(h_2\) are the water depths in the upper and lower layers, respectively, \(u_1\) and \(u_2\) are the corresponding velocities, \(q_1=h_1u_1\) and \(q_2=h_2u_2\) represent the corresponding discharges, \(Z(x)\) and \(g\) are the same as in §3.4, and \(r=\frac{\rho_1}{\rho_2}<1\) is the ratio of the constant densities \(\rho_1\) (upper layer) and \(\rho_2\) (lower layer). The studied two-layer shallow water system admits steady-state solutions satisfying \(M(\boldsymbol{U})\boldsymbol{E}(\boldsymbol{U})_x=\boldsymbol{0}\) with \[\begin{align} M(\boldsymbol{U})=\begin{pmatrix}1&0&0&0\\u_1&h_1&0&0\\0&0&1&0\\0&0&u_2&h_2\end{pmatrix},\quad \boldsymbol{E}(\boldsymbol{U})=\begin{pmatrix}q_1\\{\cal E}_1\\q_2\\{\cal E}_2\end{pmatrix},\qquad \begin{aligned} &{\cal E}_1:=\frac{q_1^2}{2h_1^2}+g(h_1+h_2+Z),\\ &{\cal E}_2:=\frac{q_2^2}{2h_2^2}+g(rh_1+h_2+Z). \end{aligned} \end{align}\]

In order to apply the fifth-order Ai-WENO-Z interpolation to the equilibrium variables, we first compute \[\begin{align} ({\cal E}_1)_j&=\frac{(q_1)_j^2}{2(h_1)_j^2}+g\big[(h_1)_j+(h_2)_j+Z_j\big],\\ ({\cal E}_2)_j&=\frac{(q_2)_j^2}{2(h_2)_j^2}+g\big[r(h_1)_j+(h_2)_j+Z_j\big], \end{align}\] and evaluate the matrices \[\begin{align} C_j=\begin{pmatrix}(u_1)_j&(h_1)_j&0&0\\g&(u_1)_j&g&0\\0&0&(u_2)_j&(h_2)_j\\rg&0&g&(u_2)_j\end{pmatrix}, \end{align}\] where \((u_1)_j=(q_1)_j/(h_1)_j\) and \((u_2)_j=(q_2)_j/(h_2)_j\). We then compute the matrices \(Q_j\) and \(Q_j^{-1}\) numerically and implement the LCD of the equilibrium variables followed by the fifth-order Ai-WENO-Z interpolation giving \(\boldsymbol{\Gamma}^\pm_{j\mp{\frac{1}{2}}}\) and then \(\boldsymbol{E}^\pm_{j\mp{\frac{1}{2}}} =\big((q_1)^\pm_{j\mp{\frac{1}{2}}},({\cal E}_1)^\pm_{j\mp{\frac{1}{2}}},(q_2)^\pm_{j\mp{\frac{1}{2}}},({\cal E}_2)^\pm_{j\mp{\frac{1}{2}}}\big)^\top\). After that, we apply the same fifth-order Ai-WENO-Z interpolation to obtain the one-sided values of the bottom topography \(Z^\pm_{j\mp{\frac{1}{2}}}\), and then solve the nonlinear equations as described in [17] to obtain \((h_i)^\pm_{j+\frac{1}{2}}\) and \((\widehat h_i)^\pm_{j+\frac{1}{2}}\), \(i=1,2\).

3.6 Application to the 1-D Euler Equations with Gravitation↩︎

In this section, we consider the 1-D compressible Euler equations with gravitation, which can be written as (6 ) with \[\boldsymbol{U}=(\rho,m,{\cal E})^\top,\quad\boldsymbol{F}(\boldsymbol{U})=\left(m,\rho u^2+p,u({\cal E}+p)\right)^\top,\quad B(\boldsymbol{U})\equiv0,\quad \boldsymbol{S}(\boldsymbol{U})=(0,-\rho \phi_x,0)^\top,\] where \(\rho\) is the density, \(u\) is the velocity, \(m:=\rho u\) is the momentum, \({\cal E}:=E+\rho\phi\), \(E\) is the total energy, \(p\) is the pressure, and \(\phi(x)\) is the time-independent gravitational potential. The system is completed with an equation of state (EOS), which, in the case of ideal gas, reads as \(E=\frac{p}{\gamma-1}+{\frac{1}{2}}\rho u^2,\) where \(\gamma\) is a specific heat ratio. The studied 1-D Euler equations with gravitation admit steady-state solutions satisfying \(M(\boldsymbol{U})\boldsymbol{E}(\boldsymbol{U})_x=\boldsymbol{0}\) with \[M(\boldsymbol{U})=\begin{pmatrix}1&0&0\\0&1&0\\L&0&m\end{pmatrix},\quad\boldsymbol{E}(\boldsymbol{U})=\begin{pmatrix}m\\K\\L\end{pmatrix},\quad K:=\rho u^2+p+\int\limits_{\hat{x}}^x{\rho\phi_x}\,{\rm d}\xi,\quad L:=\frac{{\cal E}+p}{\rho}.\]

In order to apply the fifth-order Ai-WENO-Z interpolation to the equilibrium variables, we first compute \[K_j=\rho_ju_j^2+p_j+{\cal I}_j,\quad L_j=\frac{{\cal E}_j+p_j}{\rho_j}\] where \({\cal I}_j\) is a fifth-order approximation of the integral \(\int_{x_{-\frac{5}{2}}}^{x_j}\rho\phi_x\,{\rm d}x\), in which we have set \(\hat{x}=x_{-\frac{5}{2}}\). The values \({\cal I}_j\) are computed as follows. First, we evaluate \(f_j:=\rho_j\phi_x(x_j)\) for \(j=1,\ldots,N\) and extend these values to \(j=-5,\ldots,0\) and \(j=N+1,\ldots,N+5\) using the prescribed boundary conditions implemented within the ghost cell framework. We then compute the values \({\cal I}_j\), \(j=-2,\dots,N+3\) using (22 )–(23 ).

We then evaluate the matrices \[\begin{align} &C_j=\begin{pmatrix}0&1&0\\(\gamma-2)u_j^2+c_j^2&(3-\gamma)u_j&(\gamma-1)m_j\\[1.2ex] \dfrac{(\gamma-1)u_j^2+c_j^2}{\rho_j}&\dfrac{(1-\gamma)u_j}{\rho_j}&\gamma u_j\end{pmatrix},\quad Q_j=\begin{pmatrix}\dfrac{(1-\gamma)m_j}{c_j^2}&-\dfrac{\rho_j}{c_j}&\dfrac{\rho_j}{c_j}\\[1.5ex] \dfrac{(1-\gamma)m_j^2}{c_j^2\rho_j}&\rho_j-\dfrac{m_j}{c_j}&\rho_j+\dfrac{m_j}{c_j}\\[1.5ex]1&1&1\end{pmatrix},\\ &Q_j^{-1}=\begin{pmatrix}\dfrac{u_j}{\rho_j}&-\dfrac{1}{\rho_j}&1\\[1.5ex] \dfrac{(1-\gamma)u_j^2}{2c_j\rho_j}-\dfrac{1}{2\rho_j}(u_j+c_j)&\dfrac{(\gamma-1)u_j}{2c_j\rho_j}+\dfrac{1}{2\rho_j}& \dfrac{(1-\gamma)u_j}{2 c_j}\\[1.8ex] \dfrac{(\gamma-1)u_j^2}{2c_j\rho_j}-\dfrac{1}{2\rho_j}(u_j-c_j)&\dfrac{(1-\gamma)u_j}{2c_j\rho_j}+\dfrac{1}{2\rho_j}& \dfrac{(\gamma-1)u_j}{2c_j}\end{pmatrix}, \end{align}\] where \(c_j:=\sqrt{\gamma p_j/\rho_j}\), and implement the LCD of the equilibrium variables followed by the fifth-order Ai-WENO-Z interpolation giving \(\boldsymbol{\Gamma}^\pm_{j\mp{\frac{1}{2}}}\) and then \(m^\pm_{j\mp{\frac{1}{2}}}\), \(K^\pm_{j\mp{\frac{1}{2}}}\), and \(L^\pm_{j\mp{\frac{1}{2}}}\). After that, we solve the quadratic equations \[(\gamma-1)\big(L_{j+\frac{1}{2}}^\pm-\phi(x_{j+\frac{1}{2}})\big)\big(\rho_{j+\frac{1}{2}}^\pm\big)^2-\gamma\big(K_{j+\frac{1}{2}}^\pm-{\cal I}_{j+\frac{1}{2}}\big)\rho_{j+\frac{1}{2}}^\pm+ \frac{(\gamma+1)\big(m_{j+\frac{1}{2}}^\pm\big)^2}{2}=0 \label{34613}\tag{27}\] using the method described in [26] to obtain \(\rho^\pm_{j\mp{\frac{1}{2}}}\) and then \({\cal E}^\pm_{j\mp{\frac{1}{2}}}\). In (27 ), the integrals \({\cal I}_{j+\frac{1}{2}}\) are evaluated recursively using (26 ). For the numerical flux \(\boldsymbol{{\cal K}}^{\rm FV}_{j+\frac{1}{2}}\), we use the one given in [26].

3.7 Application to the 2-D Euler Equations with Gravitation↩︎

In this section, we consider the 2-D Euler equations with gravitation, which read as (1 ) with \[\begin{align} &\boldsymbol{U}=(\rho,m,n,{\cal E})^\top\!,~\boldsymbol{F}(\boldsymbol{U})=\Big(m,\rho u^2+p,\frac{mn}{\rho},u({\cal E}+p)\Big)^\top\!,~ \boldsymbol{G}(\boldsymbol{U})=\Big(n,\frac{mn}{\rho},\rho v^2+p,v({\cal E}+p)\Big)^\top\!,\\ &B^x(\boldsymbol{U})=B^y(\boldsymbol{U})\equiv0,\quad\boldsymbol{S}(\boldsymbol{U})=(0,-\rho\phi_x,-\rho\phi_y,0)^\top, \end{align}\] where \(v\) is the \(y\)-velocity, \(n:=\rho v\) is the \(y\)-momentum, and the rest of the notation is the same as in the 1-D case. The system is completed with the EOS for an ideal gas: \(E=\frac{p}{\gamma-1}+{\frac{1}{2}}\rho(u^2+v^2)\).

The studied 2-D Euler equations with gravitation admit steady-state solutions satisfying
\(M^x(\boldsymbol{U})\boldsymbol{E}^x(\boldsymbol{U})_x=\boldsymbol{0}\) and \(M^y(\boldsymbol{U})\boldsymbol{E}^y(\boldsymbol{U})_y=\boldsymbol{0}\) with \[\begin{align} &M^x(\boldsymbol{U})= \begin{pmatrix}1&0&0&0\\0&0&1&0\\v+(\gamma+1)u\psi^x&u+(1-\gamma)v\psi^x&-\gamma\psi^x&(\gamma-1)\rho\psi^x\\L&0&0&m\end{pmatrix},\\[0.5ex] &M^y(\boldsymbol{U})= \begin{pmatrix}0&1&0&0\\v+(1-\gamma)u\psi^y&u+(\gamma+1)v\psi^y&-\gamma\psi^y&(\gamma-1)\rho\psi^y\\0&0&1&0\\0&L&0&n\end{pmatrix},\\[0.5ex] &\boldsymbol{E}^x(\boldsymbol{U})=(m,n,K^x,L)^\top,\quad\boldsymbol{E}^y(\boldsymbol{U})=(m,n,K^y,L)^\top, \end{align}\] where \[\begin{align} &K^x:=\rho u^2+p+\int\limits_{\hat{x}}^x{\rho\phi_x}\,{\rm d}\xi,\quad K^y:=\rho v^2+p+\int\limits_{\hat{y}}^y{\rho\phi_y}\,{\rm d}\eta,\quad L:=\frac{{\cal E}+p}{\rho},\\ &\psi^x:=\frac{2mn}{(\gamma-1)(2\rho^2L+n^2)-(\gamma+1)m^2},\quad\psi^y:=\frac{2mn}{(\gamma-1)(2\rho^2L+m^2)-(\gamma+1)n^2}. \end{align}\]

In order to apply the fifth-order Ai-WENO-Z interpolation to the equilibrium variables, we first compute \[K_{j,k}^x=\rho_{j,k}u_{j,k}^2+p_{j,k}+{\cal I}_{j,k}^x,\quad K_{j,k}^y=\rho_{j,k}v_{j,k}^2+p_{j,k}+{\cal I}_{j,k}^y,\quad L_{j,k}=\frac{{\cal E}_{j,k}+p_{j,k}}{\rho_{j,k}},\] where \(u_{j,k}=m_{j,k}/\rho_{j,k}\), \(v_{j,k}=n_{j,k}/\rho_{j,k}\), and \({\cal I}_{j,k}^x\) and \({\cal I}_{j,k}^y\) are fifth-order approximations of the integral \(\int_{x_{-\frac{5}{2}}}^{x_j}\rho\phi_x\,{\rm d}x\) and \(\int_{y_{-\frac{5}{2}}}^{y_k}\rho\phi_y\,{\rm d}y\), respectively. The values \({\cal I}_{j,k}^x\) are computed as follows. For each \(k=1,\dots,N_y\), we evaluate \(f_{j,k}:=\rho_{j,k}\phi_x(x_j,y_k)\) for \(j=1,\ldots,N_x\), and extend these values to \(j=-5,\ldots,0\) and \(j=N+1,\ldots,N+5\) using the prescribed boundary conditions implemented within the ghost cell framework. We then compute \({\cal I}_{j,k}^x\), \(j=-2,\dots,N_x+3\) using (22 )–(23 ). The values of \({\cal I}_{j,k}^y\), \(j=1,\ldots,N_x\), \(k=-2,\dots,N_y+3\) can be computed similarly.

We then take the local parts of \(\boldsymbol{E}^x(\boldsymbol{U})\) and \(\boldsymbol{E}^y(\boldsymbol{U})\), which are \(W^x(\boldsymbol{U})=(m,n,\rho u^2+p,L)\) and \(W^y(\boldsymbol{U})=(m,n,\rho v^2+p,L)\), and evaluate their Jacobians \[\begin{align} &\frac{\partial W^x}{\partial U}(\boldsymbol{U})=\begin{pmatrix}0&1&0&0\\0&0&1&0\\[0.5ex] (1-\gamma)\phi+\dfrac{\gamma-3}{2}u^2+\dfrac{\gamma-1}{2}v^2 &(3-\gamma)u&(1-\gamma)v&\gamma-1\\[1.5ex] (\gamma-1)\dfrac{u^2+v^2}{\rho}-\dfrac{\gamma{\cal E}}{\rho^2}&(1-\gamma)\dfrac{u}{\rho}&(1-\gamma)\frac{v}{\rho}&\dfrac{\gamma}{\rho} \end{pmatrix},\\[1.5ex] &\frac{\partial W^y}{\partial U}(\boldsymbol{U})=\begin{pmatrix}0&1&0&0\\0&0&1&0\\[0.5ex] (1-\gamma)\phi+\dfrac{\gamma-1}{2}u^2+\dfrac{\gamma-3}{2}v^2&(1-\gamma)u&(3-\gamma)v&\gamma-1\\[1.5ex] (\gamma-1)\dfrac{u^2+v^2}{\rho}-\dfrac{\gamma{\cal E}}{\rho^2}&(1-\gamma)\dfrac{u}{\rho}&(1-\gamma)\frac{v}{\rho}&\dfrac{\gamma}{\rho} \end{pmatrix}. \end{align}\]

Notice that computing the matrices \(C^x_{j,k}=\frac{\partial W^x}{\partial U}(\boldsymbol{U}_{j,k})M^x(\boldsymbol{U}_{j,k})\) and \(C^y_{j,k}=\frac{\partial W^y}{\partial U}(\boldsymbol{U}_{j,k})M^y(\boldsymbol{U}_{j,k})\) analytically is quite cumbersome and their evaluation will be computationally expensive. Moreover, their analytical eigenstructures are unavailable. We therefore compute \(C^x_{j,k}\) and \(C^y_{j,k}\) together with the matrices \(Q^x_{j,k}\), \((Q^x_{j,k})^{-1}\), \(Q^y_{j,k}\), and \((Q^y_{j,k})^{-1}\) numerically. We then implement the LCD of the equilibrium variables followed by the fifth-order Ai-WENO-Z interpolation in the \(x\)- and \(y\)-directions separately. In the \(x\)-direction, this gives \(\boldsymbol{\Gamma}^\pm_{j\mp{\frac{1}{2}},k}\) and then \(m^\pm_{j\mp{\frac{1}{2}},k}\), \(n^\pm_{j\mp{\frac{1}{2}},k}\), \((K^x)^\pm_{j\mp{\frac{1}{2}},k}\), and \(L^\pm_{j\mp{\frac{1}{2}},k}\). After that, we solve the quadratic equations \[\begin{align} (\gamma-1)\big(L^\pm_{{j+\frac{1}{2}},k}-\phi(x_{j+\frac{1}{2}},y)\big)\big(\rho^\pm_{{j+\frac{1}{2}},k}\big)^2-\gamma\big((K^x)^\pm_{{j+\frac{1}{2}},k}-I^x_{{j+\frac{1}{2}},k}\big) \rho^\pm_{{j+\frac{1}{2}},k}&\\ +\frac{(\gamma+1)\big(m^\pm_{{j+\frac{1}{2}},k}\big)^2}{2}-\frac{(\gamma-1)\big(n^\pm_{{j+\frac{1}{2}},k}\big)^2}{2}&=0 \end{align} \label{34614}\tag{28}\] using the method described in [26] to obtain \(\rho^\pm_{j\mp{\frac{1}{2}},k}\) and then \({\cal E}^\pm_{j\mp{\frac{1}{2}},k}\). In (28 ), the integrals \(I^x_{{j+\frac{1}{2}},k}\) are evaluated recursively using (26 ). The point values \(m^\pm_{j,k\mp{\frac{1}{2}}}\), \(n^\pm_{j,k\mp{\frac{1}{2}}}\), \((K^y)^\pm_{j,k\mp{\frac{1}{2}}}\), \(L^\pm_{j,k\mp{\frac{1}{2}}}\), \(\rho^\pm_{j,k\mp{\frac{1}{2}}}\), and \({\cal E}^\pm_{j,k\mp{\frac{1}{2}}}\) can be computed using a similar interpolation procedure carried out in the \(y\)-direction. Finally, we use the numerical fluxes described in [26] for \(\boldsymbol{{\cal K}}_{{j+\frac{1}{2}},k}\) and \(\boldsymbol{{\cal L}}_{j,{k+\frac{1}{2}}}\).

3.8 Application to the 2-D Ripa System↩︎

In this section, we consider the 2-D Ripa system, which was introduced in [30], [31] and reads as (1 ) with \[\begin{align} &\boldsymbol{U}=(h,q^x,q^y,h\theta,Z)^\top\!,~\boldsymbol{F}(\boldsymbol{U})=\Big(q^x,q^xu+P,\frac{q^xq^y}{h},q^x\theta,0\Big)^\top\!,~ \boldsymbol{G}(\boldsymbol{U})=\Big(q^y,\frac{q^xq^y}{h},q^yv+P,q^y\theta, 0\Big)^\top\!,\\ &B^x(\boldsymbol{U})=\begin{pmatrix}0&0&0&0&0\\0&0&0&0&-h\theta\\0&0&0&0&0\\0&0&0&0&0\\0&0&0&0&0\end{pmatrix},\quad B^y(\boldsymbol{U})=\begin{pmatrix}0&0&0&0&0\\0&0&0&0&0\\0&0&0&0&-h\theta\\0&0&0&0&0\\0&0&0&0&0\end{pmatrix},\quad\boldsymbol{S}(\boldsymbol{U})\equiv\boldsymbol{0}, \end{align}\] where \(h\) is the water depth, \(u\) and \(v\) are the velocities in \(x\)- and \(y\)-directions, respectively, \(q^x=hu\), \(q^y=hv\), \(Z(x,y)\) is the bottom topography, \(\theta\) is the potential temperature, and \(P={\frac{1}{2}}h^2\theta\) is the pressure.

We restrict our consideration to the quasi 1-D moving-water equilibria: the \(x\)-directional, \[\begin{align} &q^x_x=q^y=\theta_x=({\cal E}^x)_x\equiv0,\quad{\cal E}^x=\frac{u^2}{2}+\theta(h+Z)+{\cal I}^x,\\ &{\cal I}^x=-\int\limits_{\widehat x}^x\left[\sqrt{2P(\xi,y,t)}\big(\sqrt{\theta(\xi,y,t)}\big)_\xi+Z(\xi,y)\theta_\xi(\xi,y,t)\right] {\rm d}\xi, \end{align} \label{4467}\tag{29}\] and the \(y\)-directional, \[\begin{align} &q^x=q^y_y=\theta_y=({\cal E}^y)_y\equiv0,\quad{\cal E}^y=\frac{v^2}{2}+\theta(h+Z)+{\cal I}^y,\\ &{\cal I}^y=-\int\limits_{\widehat y}^y\left[\sqrt{2P(x,\eta,t)}\big(\sqrt{\theta(x,\eta,t)}\big)_\eta+Z(x,\eta)\theta_\eta(x,\eta,t) \right]{\rm d}\eta, \end{align} \label{4468}\tag{30}\] ones. The steady states (29 ) and (30 ) satisfy \(M^x(\boldsymbol{U})\boldsymbol{E}^x(\boldsymbol{U})_x=\boldsymbol{0}\) and \(M^y(\boldsymbol{U})\boldsymbol{E}^y(\boldsymbol{U})_y=\boldsymbol{0}\), respectively, with \[\begin{align} &M^x(\boldsymbol{U})=\begin{pmatrix}1&0&0&0&0\\u&h&0&0&0\\v&0&q^x&0&0\\\theta&0&0&q^x&0\\0&0&0&0&0\end{pmatrix},\quad M^y(\boldsymbol{U})=\begin{pmatrix}1&0&0&0&0\\u&0&q^y&0&0\\v&h&0&0&0\\\theta&0&0&q^y&0\\0&0&0&0&0\end{pmatrix},\quad \begin{matrix}\boldsymbol{E}^x(\boldsymbol{U})=\Big(q^x,{\cal E}^x,v,\theta,Z\Big)^\top\!,\\ \boldsymbol{E}^y(\boldsymbol{U})=\Big(q^y,{\cal E}^y,u,\theta,Z\Big)^\top\!. \end{matrix} \end{align}\] For the 2-D Ripa system, the matrices \(C^x(\boldsymbol{U})\) and \(C^y(\boldsymbol{U})\) can be explicitly computed: \[C^x(\boldsymbol{U})=\begin{pmatrix}u&h&0&0&0\\[0.5ex]\theta&u&0&\dfrac{q^x}{2}&0\\[1.5ex]0&0&u&0&0\\0&0&0&u&0\\0&0&0&0&0\end{pmatrix},\quad C^y(\boldsymbol{U})=\begin{pmatrix}v&h&0&0&0\\\theta&v&0&\dfrac{q^y}{2}&0\\0&0&v&0&0\\0&0&0&v&0\\0&0&0&0&0\end{pmatrix}, \label{44614}\tag{31}\] and they have complete eigensystems so that the corresponding matrices \(Q^x(\boldsymbol{U})\) and \(Q^y(\boldsymbol{U})\) along with their inverses are given by \[\begin{align} &Q^x(\boldsymbol{U})=\begin{pmatrix}0&0&-\dfrac{q^x}{2\theta}&\dfrac{c}{\theta}&-\dfrac{c}{\theta}\\[1.5ex]0&0&0&1&1\\0&1&0&0&0\\0&0&1&0&0\\1&0&0&0&0 \end{pmatrix},\quad[Q^x(\boldsymbol{U})]^{-1}=\begin{pmatrix}0&0&0&0&1\\0&0&1&0&0\\0&0&0&1&0\\[0.5ex] \dfrac{\theta}{2c}&\dfrac{1}{2}&0&-\dfrac{q^x}{4c}&0\\[2.0ex]-\dfrac{\theta}{2c}&\dfrac{1}{2}&0&-\dfrac{q^x}{4c}&0\end{pmatrix},\\[1.5ex] &Q^y(\boldsymbol{U})=\begin{pmatrix}0&0&-\dfrac{q^y}{2\theta}&\dfrac{c}{\theta}&-\dfrac{c}{\theta}\\[1.5ex]0&0&0&1&1\\0&1&0&0&0\\0&0&1&0&0\\1&0&0&0&0 \end{pmatrix},\quad[Q^y(\boldsymbol{U})]^{-1}=\begin{pmatrix}0&0&0&0&1\\0&0&1&0&0\\0&0&0&1&0\\[0.5ex] \dfrac{\theta}{2c}&\dfrac{1}{2}&0&-\dfrac{q^y}{4c}&0\\[2.0ex]-\dfrac{\theta}{2c}&\dfrac{1}{2}&0&-\dfrac{q^y}{4c}&0\end{pmatrix}. \end{align}\] In order to apply the fifth-order Ai-WENO-Z interpolation to the equilibrium variables, we first compute \[{\cal E}^x_{j,k}=\frac{u_{j,k}^2}{2}+\theta_{j,k}(h_{j,k}+Z_{j,k})+{\cal I}^x_{j,k},\quad {\cal E}^y_{j,k}=\frac{v_{j,k}^2}{2}+\theta_{j,k}(h_{j,k}+Z_{j,k})+{\cal I}^y_{j,k},\] where \({\cal I}^x_{j,k}\) and \({\cal I}^y_{j,k}\) are computed using the path-conservative integration as in [16]. We then implement the LCD of the equilibrium variables followed by the fifth-order Ai-WENO-Z interpolation in the \(x\)- and \(y\)-directions separately. In the \(x\)-direction, this gives \(\boldsymbol{\Gamma}^\pm_{j\mp{\frac{1}{2}},k}\) and then \(({\cal E}^x)_{{j+\frac{1}{2}},k}^\pm\), \((q^x)^\pm_{j\mp{\frac{1}{2}},k}\), \(v^\pm_{j\mp{\frac{1}{2}},k}\), \(\theta^\pm_{j\mp{\frac{1}{2}},k}\), and \(Z^\pm_{j\mp{\frac{1}{2}},k}\). After that, we solve the quadratic equations \[{\frac{1}{2}}\left(\frac{(q^x)_{{j+\frac{1}{2}},k}^\pm}{h_{{j+\frac{1}{2}},k}^\pm}\right)^2+\theta_{{j+\frac{1}{2}},k}^\pm\big(h_{{j+\frac{1}{2}},k}^\pm+Z_{{j+\frac{1}{2}},k}^\pm\big)+ ({\cal I}^x)_{{j+\frac{1}{2}},k}^\pm-({\cal E}^x)_{{j+\frac{1}{2}},k}^\pm=0 \label{34624f}\tag{32}\] to obtain \(h^\pm_{j\mp{\frac{1}{2}},k}\). In (32 ), the values \(({\cal I}^x)_{{j+\frac{1}{2}},k}^\pm\) are evaluated recursively also using the path-conservative integration. A similar interpolation procedure is carried out in the \(y\)-direction. Finally, we use the numerical fluxes for \(\boldsymbol{{\cal K}}_{{j+\frac{1}{2}},k}\) and \(\boldsymbol{{\cal L}}_{j,{k+\frac{1}{2}}}\). We refer the reader to [32] for details.

4 Numerical Examples↩︎

In this section, we test the proposed fifth-order WB A-WENO scheme based on the LCD of equilibrium variables on several numerical examples for the nozzle flow system, one- and two-layer shallow water equations, compressible Euler equations with gravitation, and the 2-D Ripa system. For the sake of brevity, this scheme will be referred to as Scheme 1 and its performance will be compared with the following two schemes:

\(\bullet\) Scheme 2: The WB A-WENO scheme from [17], in which the Ai-WENO-Z interpolation is applied to the equilibrium variables \(\boldsymbol{E}\) without any LCD;

\(\bullet\) Scheme 3: The A-WENO scheme from [28], which can only preserve the simplest “lake-at-rest” steady states, but uses the LCD applied to the conservative variables \(\boldsymbol{U}\) to reduce the spurious oscillations.

In all of the examples, we solve the ODE systems (7 ) and (8 ) using the three-stage third-order SSP Runge–Kutta method (see, e.g., [33], [34]) with the time-step restricted by the CFL number \(0.5\).

When the eigendecomposition of matrices \(C\), \(C^x\), or \(C^y\) is required, we perform it using the FORTRAN LAPACK routine GEEV.

4.1 Nozzle Flow System↩︎

Example 1—Flow in Continuous Divergent Nozzle↩︎

In the first example taken from [17], [22], we consider the divergent nozzle described using the smooth cross-section \[\sigma(x)=0.976+0.748\tanh(0.8x-4).\] We first take the steady states with \(q_{\rm eq}(x)\equiv8\), \({\cal E}_{\rm eq}(x)\equiv21.9230562619897\) and compute the discrete values of \(\rho_{\rm eq}(x)\) by solving the corresponding nonlinear equations; see [17], [22]. We then obtain \(u_{\rm eq}(x)=q_{\rm eq}(x)/(\sigma(x)\rho_{\rm eq}(x))\).

Equipped with these steady states, we add a small perturbation to the density field and consider the initial data \[\rho(x,0)=\rho_{\rm eq}(x)+\begin{cases}10^{-2},&x\in[0.5,1.5],\\0,&otherwise,\end{cases}\qquad q(x,0)=\sigma(x)\rho(x,0)u_{\rm eq}(x),\] which are prescribed in the computational domain \([0,10]\) subject to the homogeneous Neumann boundary conditions.

We compute the numerical solutions until the final time \(t=0.8\) by Schemes 1–3 on a uniform mesh with \(\Delta x=1/20\) and plot the differences \(\rho(x,0.8)-\rho_{\rm eq}(x)\) and \(q(x,0.8)-q_{\rm eq}(x)\) together with the reference solution computed by Scheme 1 on a much finer mesh with \(\Delta x=1/400\) in Figure 1. As one can see, Scheme 1 clearly outperforms Schemes 2 and 3 as there are no oscillations in the results computed by Scheme 1. We also stress that in this example, Scheme 3, which is not WB, produces the largest oscillations as we capture small perturbations of the steady state.

Example 2—Flow in Continuous Convergent Nozzle↩︎

In the second example also taken from [17], [22], we consider the convergent nozzle described using the smooth cross-section \[\sigma(x)=0.976-0.748\tanh(0.8x-4).\] We first take the steady states with \(q_{\rm eq}(x)\equiv8\), \({\cal E}_{\rm eq}(x)\equiv58.3367745090349\), compute the discrete values of \(\rho_{\rm eq}(x)\) by solving the corresponding nonlinear equations, and then obtain \(u_{\rm eq}(x)=q_{\rm eq}(x)/(\sigma(x)\rho_{\rm eq}(x))\); see [17], [22]. We then consider the initial data containing a substantially larger perturbation than the one studied in Example 1: \[\rho(x,0)=\rho_{\rm eq}(x)+\begin{cases}0.3,&x\in[0.5,1.5],\\0,&otherwise,\end{cases}\qquad q(x,0)=\sigma(x)\rho(x,0)u_{\rm eq}(x).\] As in Example 1, the computational domain is \([0,10]\) and the homogeneous Neumann boundary conditions are imposed.

We compute the numerical solutions until the final time \(t=0.5\) by Schemes 1–3 on a uniform mesh with \(\Delta x=1/20\) together with the reference solution computed by Scheme 1 on a much finer mesh with \(\Delta x=1/400\) and plot the differences \(\rho(x,0.5)-\rho_{\rm eq}(x)\) and \(q(x,0.5)-q_{\rm eq}(x)\) in Figure 2. As one can see, Scheme 3 does not produce visible oscillations as the magnitude of the perturbation is apparently larger than the size of the truncation errors. On the contrary, Scheme 2, which does not use any LCD, now produces larger oscillations than in Example 1, where the size of perturbation was much smaller.

4.2 Saint-Venant System with Manning Friction↩︎

Example 3—Riemann Problem (\(n=0.4\))↩︎

In this example, we test the performance of Schemes 1–3 on a Riemann problem with the following initial data (prescribed in the computational domain \([-0.1,0.3]\) subject to the homogeneous Neumann boundary conditions): \[h(x,0)=\left\{\begin{align}&1,&&x<0,\\&0.8,&&x>0,\end{align}\right.\qquad u(x,0)=\left\{\begin{align}&2,&&x<0,\\&4,&&x>0,\end{align}\right.\] and the bottom topography also containing a jump at \(x=0\): \[Z(x)=\left\{\begin{align}&1,&&x<0,\\&1.9,&&x>0.\end{align}\right.\]

We compute the numerical solutions by the three studied schemes until the final time \(t=0.03\) on a uniform mesh with \(\Delta x=1/200\) together with the reference solution computed by Scheme 1 on a much finer mesh with \(\Delta x=1/10000\). The obtained results (\(h\), \(q\), and \({\cal E}\)) are shown in Figure 3, where one can observe that the solutions computed by both Schemes 2 and 3 are oscillatory, whereas Scheme 1 solution is oscillation-free. In order to further study the performance of Schemes 1–3, we refine the mesh to \(\Delta x=1/2500\) and compare the high-resolution solutions of the three studied schemes together with the reference solution computed by Scheme 1 on a much finer mesh with \(\Delta x=1/10000\); see Figure 4. One can observe that the high-frequency oscillations produced by Scheme 3 on a coarse mesh were almost suppressed when the mesh was refined, whereas Scheme 2, which does not employ any LCD, still produces oscillations, which can be clearly seen in the zoom views shown in the right column.

Example 4—Convergence to a Steady State (\(n=0.15\))↩︎

In this example, we study the convergence of the solutions computed by Schemes 1–3 towards the steady flow over a hump. We consider the continuous bottom topography given by \[Z(x)=\left\{\begin{align} &0.2&&if~8\le x\le12,\\ &0&&otherwise, \end{align}\right.\] and the initial and boundary data that correspond to a subcritical flow: \[h(x,0)\equiv2-Z(x),\quad q(x,0)\equiv0,\quad q(0,t)=4.42,\quad h(25,t)=2,\] with the boundary conditions for \(h\) at \(x=0\) and \(q\) at \(x=25\) set to be homogeneous Neumann.

We compute the numerical solutions until the final time \(t=500\) on the computational domain \([0,25]\) covered by a uniform mesh with \(\Delta x=1/4\) together with the reference solution computed by Scheme 1 on a much finer mesh with \(\Delta x=1/80\). The obtained numerical solutions \(h+Z\), \(q\), and \({\cal E}-I_{\frac{1}{2}}\) (note that we subtract \(I_{\frac{1}{2}}\) from \(\cal E\) to obtain the quantity, which is supposed to be constant disregarding what values of \(\widehat x\) were used in the evaluation of \({\cal E}_j\)) are plotted in Figure 5. One can clearly see that the WB Schemes 1 and 2 converge to the constant \(q\) and \({\cal E}\), whereas the non-WB Scheme 3 generates spurious oscillations.

We then test the ability of the studied schemes to capture the propagation of a small perturbation of the obtained moving-water equilibria. To this end, we denote the obtained steady states by \(h_{\rm eq}(x)\) and \(q_{\rm eq}(x)\) (notice that each scheme has its own discrete equilibrium), and then consider the following initial data: \[h(x,0)=h_{\rm eq}(x)+\left\{\begin{align}&10^{-4},&&9.5\le x\le10.5,\\&0,&&otherwise,\end{align}\right.\qquad q(x,0)=q_{\rm eq}(x).\] We compute the solutions by the three studied schemes until the final time \(t=1.5\) on the same uniform mesh with \(\Delta x=1/4\). The obtained differences \(h(x,1.5)-h_{\rm eq}(x)\) are plotted in Figure 6. One can clearly see that both the WB Schemes 1 and 2 can capture the time evolution of the perturbation quite accurately, whereas the non-WB Scheme 3 generates spurious oscillations as it can preserve still-water equilibria only. Even though Scheme 2 is WB, it still generates some oscillations, and thus Scheme 1 clearly outperforms both of its counterparts in this example.

Example 5—Efficiency vs. Accuracy (\(n=0\))↩︎

In this example, we assess both the WB property and the computational efficiency of Schemes 1–3 on a nontrivial moving-water steady state. The bottom topography is \[Z(x)=-0.2e^{-40(x-10)^2},\] the equilibrium state is \[E_{\rm eq}(x)\equiv32,\quad q_{\rm eq}(x)\equiv2,\] and the discrete values of the corresponding equilibrium water depth \(h_{\rm eq}(x_j)\) are obtained by solving the nonlinear equations \[\frac{4}{2h_{\rm eq}^2(x_j)}+g\left(h_{\rm eq}(x_j)+Z(x_j)\right)=32.\] These data are prescribed in the computational domain \([0,25]\) subject to the homogeneous Neumann boundary conditions.

To test the ability of the three schemes to preserve this steady state and capture small perturbations, we take the following initial data: \[h(x,0)=h_{\rm eq}(x)+10^{-4}e^{-4(x-12)^2},\quad q(x,0)=q_{\rm eq}(x)\equiv2,\] in which a small Gaussian-shaped perturbation is added to the equilibrium state.

We compute the numerical solutions up to the final time \(t=1\) on a sequence of uniform meshes with \(\Delta x=1/2\), \(1/4\), \(1/8\), \(1/16\), \(1/32\), and \(1/64\). We measure the \(L^1\)-errors in \(h\) using the Runge formula, which is based on the solutions computed on the three consecutive uniform grids with the mesh sizes \(\Delta x\), \(2\Delta x\), and \(4\Delta x\) and denoted by \((\cdot)^{\Delta x}\), \((\cdot)^{2\Delta x}\), and \((\cdot)^{4\Delta x}\), respectively: \[{\rm Error}(\Delta x)\approx\frac{\delta_{12}^2}{|\delta_{12}-\delta_{24}|}, \label{runge}\tag{33}\] where \(\delta_{12}:=\|(\cdot)^{\Delta x}-(\cdot)^{2\Delta x}\|_{L^1}\) and \(\delta_{24}:=\|(\cdot)^{2\Delta x}-(\cdot)^{4\Delta x}\|_{L^1}\). We report the \(L^1\)-errors along the CPU times (in seconds) in Table 1, where one can see that Scheme 3 is the least computationally expensive, but it produces errors several orders of magnitude larger than those produced by Schemes 1 and 2. Concerning the WB property, it is also evident that Schemes 1 and 2 capture the small perturbation of the steady state very accurately, while Scheme 3 fails to suppress spurious oscillations, resulting in substantially larger errors. The latter can be clearly seen in Figure 7, where we plot \(h(x,1)-h_{eq}(x)\) computed by the three studied schemes on coarse meshes with \(\Delta x=1/16\) and \(1/32\) (here, the finest-mesh solution computed by Scheme 1 is being considered as the reference solution).

4.3 Two-Layer Shallow Water System↩︎

Example 6—Experimental Order of Accuracy↩︎

In this example taken from [35], we consider the following initial data: \[h_1(x,0)=5+e^{\cos(2\pi x)},\quad h_2(x,0)=5-e^{\cos(2\pi x)}-\sin^2(\pi x),\quad q_1(x,0)=q_2(x,0)\equiv0,\] and a continuous bottom topography \[Z(x)=\sin^2(\pi x)-10,\] both prescribed in the computational domain \([0,1]\) subject to the periodic boundary conditions.

We compute the numerical solution until the final time \(t=0.1\) by the studied Schemes 1–3 on a sequence of uniform meshes with \(\Delta x=1/40\), \(1/80\), \(1/160\), \(1/320\), and \(1/640\). We measure the \(L^1\)-errors in \(h_1\) using the Runge formula (33 ) and estimate the experimental convergence rates as follows: \[{\rm Rate}(\Delta x)\approx\log_2\left(\frac{\delta_{24}}{\delta_{12}}\right).\] The obtained results are reported in Table 2, where one can clearly see that the expected orders of accuracy are achieved by all the three studied schemes. Note that in order to achieve the fifth order of accuracy, we have used smaller time steps with \(\Delta t\sim(\Delta x)^\frac{5}{3}\).

Example 7—Small Perturbation of a Discontinuous Steady State↩︎

In this example taken from [17], [22], we consider a discontinuous steady state given by \[\begin{align} (h_1)_{\rm eq}(x):&=\begin{cases}1.22373355048230,&x<0,\\1.44970064153589,&x>0,\end{cases}\qquad(q_1)_{\rm eq}(x)\equiv12,\\ (h_2)_{\rm eq}(x):&=\begin{cases}0.968329515483846,&x<0,\\1.12439026921484,&x>0,\end{cases}\qquad(q_2)_{\rm eq}(x)\equiv10,\\ \end{align}\] and a discontinuous bottom topography \[Z(x)=\begin{cases}-2,&x<0,\\-1,&x>0.\end{cases}\] In order to test the ability of the studied schemes to capture quasi-steady solutions, we add a small perturbation to the upper layer depth and take the following initial data: \[\begin{align} h_1(x,0)&=(h_1)_{\rm eq}(x)+\begin{cases}0.04,&x\in[-0.9,-0.8],\\0,&otherwise,\end{cases}\\ h_2(x,0)&=(h_2)_{\rm eq}(x),\quad q_1(x,0)=(q_1)_{\rm eq}(x),\quad q_2(x,0)=(q_2)_{\rm eq}(x), \end{align}\] prescribed in the computational domain \([-1,1]\) subject to the homogeneous Neumann boundary conditions.

We compute the numerical solutions until the final time \(t=0.1\) by Schemes 1–3 on a uniform mesh with \(\Delta x=1/100\) and obtain the reference solution using Scheme 1 on a much finer mesh with \(\Delta x=1/1000\). The differences \(h_1(x,0.1)-(h_1)_{\rm eq}(x)\) and \(h_2(x,0.1)-(h_2)_{\rm eq}(x)\) are plotted in Figure 8, where one can clearly see that, unlike the proposed Scheme 1, Schemes 2 and 3 produce oscillatory numerical results.

Example 8—Riemann Problem↩︎

In this example taken from [17], [22], we numerically solve a Riemann problem with the following initial data: \[(h_1,q_1,h_2,q_2)(x,0)=\begin{cases}(1,1.5,1,1),&x<0,\\(0.8,1.2,1.2,1.8),&otherwise,\end{cases}\] and discontinuous bottom topography: \[Z(x)=\begin{cases}-2,&x<0,\\-1.5,&otherwise,\end{cases}\] prescribed in the computational domain \([-1,1]\) subject to the homogeneous Neumann boundary conditions.

We compute the numerical solutions until the final time \(t=0.1\) by Schemes 1–3 on a uniform mesh with \(\Delta x=1/50\). The obtained upper layer depth \(h_1\) and lower layer depth \(h_2\) are plotted in Figure 9 together with the reference solution computed by Scheme 1 on a much finer mesh with \(\Delta x=1/2000\). As one can see, Scheme 1 clearly outperforms Schemes 2 and 3 as there are no oscillations in the results computed by Scheme 1.

4.4 1-D Euler Equations with Gravitation↩︎

Example 9—Shock Tube Problem↩︎

In this example, which is a modification of the example studied in [36], [37], we set a nonlinear gravitational potential \(\phi(x)=\frac{1}{10x+1}\) in the computational domain \([0,1]\), where we prescribe the following initial data: \[(\rho(x,0),u(x,0),p(x,0))=\left\{\begin{align}&(1,0,1),&&x\le0.5,\\&(0.125,0,0.1),&&x>0.5,\end{align}\right.\] and the reflecting (solid wall) boundary conditions, which are implemented using the ghost cell technique: we set the same values of \(\rho\) and \(p\) in the ghost cells while for \(u\) the sign is negated.

We compute the numerical solutions until the final time \(t=0.2\) by Schemes 1 and 2 on a uniform grid with \(\Delta x=1/30\) together with the reference solution computed by Scheme 1 on a much finer mesh with \(\Delta x=1/8000\). The numerical results are shown in Figure 10, where one can clearly see that Scheme 1 outperforms Scheme 2 as there are no oscillations in the results computed by Scheme 1.

4.5 2-D Euler Equations with Gravitation↩︎

Example 10—Discontinuous Perturbation of a Steady State↩︎

In this example, which is a modification of the example studied in [36], [38], we consider the following hydrostatic equilibrium: \[\rho_{\rm eq}(x,y)=1.21e^{-1.21\phi(x,y)},\quad u_{\rm eq}(x,y)=v_{\rm eq}(x,y)\equiv0,\quad p_{\rm eq}(x,y)=e^{-1.21\phi(x,y)},\] with the gravitational potential \(\phi(x,y)=x+y\). We take the computational domain \([0,1]\times[0,1]\) and impose the homogeneous Neumann boundary conditions.

We construct the discrete steady state as it was described in [26], but with the integrals in \(K^x\) and \(K^y\) evaluated within the fifth order of accuracy.

Equipped with the discrete steady state, we first numerically verify that Schemes 1 and 2 can preserve it within the machine accuracy, and then introduce a discontinuous pressure perturbation and consider the perturbed initial data: \[(\rho,u,v)(x,y,0)=(\rho_{\rm eq},u_{\rm eq},v_{\rm eq})(x,y),\quad p(x,y,0)=p_{\rm eq}(x,y)+ \left\{\begin{align}&0.5,&&(x-0.5)^2+(y-0.5)^2\le0.2^2,\\&0,&&otherwise.\end{align}\right.\]

We compute the numerical solutions by Schemes 1 and 2 on a uniform grid with \(\Delta x=\Delta y=1/80\) together with the reference solution computed by Scheme 1 on a much finer mesh with \(\Delta x=\Delta y=1/400\) at time \(t=0.12\) and plot the obtained results in Figure 11. As one can see, Scheme 1 clearly outperforms Scheme 2 as there are almost no oscillations in the results computed by Scheme 1.

4.6 2-D Ripa System↩︎

Example 11—Quasi 1-D Moving-Water Equilibrium and its Circular Perturbation↩︎

In this example, we take a continuous bottom topography, which varies in the \(x\)-direction only: \[Z(x)=\left\{\begin{align} &0.2-0.05(x-10)^2,&&8\le x\le12,\\ &0,&&otherwise, \end{align}\right.\] and the following subcritical initial and boundary conditions: \[\begin{align} &(\mathcal{E}^x,q^x,q^y,\theta)\Big|_{(x,y,0)}= \big(\mathcal{E}^x_{\rm eq},q^x_{\rm eq},q^y_{\rm eq},\theta_{\rm eq}\big)\Big|_{(x,y)}\equiv(110.33025,4.42\sqrt{5},0,49.06),\\ &h(25,y,t)=2,~~q^x(0,y,t)=4.42\sqrt{5}, \end{align}\] prescribed in the computational domain \([0,25]\times[0,10]\). The rest of the boundary conditions in the \(x\)-direction are homogeneous Neumann, and the reflecting (solid wall) boundary conditions are imposed in the \(y\)-direction. The latter are implemented using the ghost cell technique: we set the same values of \(h\), \(u\), and \(\theta\) in the ghost cells while for \(v\) the sign is negated.

We compute the solution by Schemes 1–3 until the final time \(t=20\) on a uniform mesh with \(\Delta x=\Delta y=1/4\) and present the equilibrium errors in Table 3. As one can see, the errors of Schemes 1 and 2 are close to the machine errors, while Scheme 3 generates large \(L^1\)-errors. We also note that the errors of Scheme 1 are smaller than those of Scheme 2, which implies that Scheme 1 performs slightly better than Scheme 2 in maintaining the moving-water equilibrium.

We then add a small genuinely 2-D circular perturbation to the water depth without changing the initial \(q^x\), \(q^y\), and \(h\theta\). The initially perturbed \(h\) is \[h(x,y,0)=h_{\rm eq}(x,y)+ \left\{\begin{align} &0.01,&&(x-6)^2+(y-5)^2<0.25,\\ &0,&&otherwise. \end{align}\right.\]

We compute the numerical solutions by Schemes 1–3 until the final time \(t=0.4\) on a uniform mesh with \(\Delta x=\Delta y=1/4\) together with the reference solution computed by Scheme 1 on a much finer mesh with \(\Delta x=\Delta y=1/40\) and plot the obtained results in Figure 16. As one can see, Scheme 3, which cannot preserve the moving-water equilibrium, produces large spurious oscillations near \(x=8\) and \(x=12\) due to the effect of the bottom topography. On the other hand, the solutions computed by Schemes 1 and 2 are almost oscillation-free.

Example 12—2-D Dam-Break Problem↩︎

In the final example, we consider the initial data \[(h,u,v,\theta)\Big|_{(x,y,0)}= \left\{\begin{align} &(5,0.5,0,9.812),&&\max\{|x|,|y|\}\le0.5,\\ &(3,2.75,0,15.2086),&&\quad~~~otherwise, \end{align}\right.\] prescribed in the computational domain \([-1,1]\times[-1,1]\) subject to the homogeneous Neumann boundary conditions. The bottom topography is flat (\(Z(x,y)\equiv0\)).

We compute the solutions by Schemes 1–3 until the final time \(t=0.075\) on a uniform mesh with \(\Delta x=\Delta y=1/100\) and plot the obtained results in Figure 17. The reference solution is computed by Scheme 1 on a much finer mesh with \(\Delta x=\Delta y=1/400\). As one can see, Scheme 2, which does not use the LCD, generates spurious oscillations, which can be clearly seen in the 1-D slices along the line \(y=0.125\).

Acknowledgment↩︎

The work of S. Chu was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - SPP 2410 Hyperbolic Balance Laws in Fluid Mechanics: Complexity, Scales, Randomness (CoScaRa) within the Project(s) HE5386/26-1 (Numerische Verfahren für gekoppelte Mehrskalenprobleme,525842915) and (Zufällige kompressible Euler Gleichungen: Numerik und ihre Analysis, 525853336) HE5386/27-1, and the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - SPP 2183: Eigenschaftsgeregelte Umformprozesse with the Project(s) HE5386/19-2,19-3 Entwicklung eines flexiblen isothermen Reckschmiedeprozesses für die eigenschaftsgeregelte Herstellung von Turbinenschaufeln aus Hochtemperaturwerkstoffen (424334423). The work of A. Kurganov was supported in part by NSFC grants 12171226 and W2431004. The work of B.-S. Wang was supported in part by NSFC grant 12301530 and the startup funding provided by the Ocean University of China.

References↩︎