Two-dimensional Boussinesq equations applied to channel flows: deducing and applying the equations

Equações de Boussinesq em duas dimensões aplicada a escoamentos em canais: dedução e aplicação das equações

ABSTRACT

A basic hypothesis adopted for theoretical formulation of fluid flows is the hydrostatic pressure distribution. However, many researchers have pointed out that this simplification can lead to errors, in cases such as dam break flow. Discrepancy between computational solution and the experiment is attributed to the pressure distribution. These findings are not new, but it is not presented any formulation in the literature that considers the non-hydrostatic pressure distribution in 2D flow. This article deduces the Boussinesq Equations as an evolution of the Shallow Water Equations with the hypothesis of non-hydrostatic pressure distribution in the vertical direction. XYZ Orthogonal Cartesian System is used, considering the influence of channel bed slope and head losses of flow. It is presented the non-hydrostatical correction in the Boussinesq equation in two dimension using Fourier series. The solution uses Runge-Kutta Discontinuous Galerkin Method and the formulation is applied to a cylindrical dam-break.

Keywords:
Boussinesq equations; Channel flow; Mathematical modeling; Non-hydrostatic pressure distribution; Two-dimensional model

RESUMO

Uma das hipóteses básicas adotada para a formulação teórica dos escoamentos de fluidos é a distribuição da pressão hidrostática na vertical. No entanto, muitos pesquisadores apontaram que essa simplificação pode levar a erros, em casos como o escoamento na ruptura de barragens. A discrepância entre a solução computacional e o experimento é atribuída à distribuição de pressão. Essas observações não são novas, mas não é apresentada qualquer formulação na literatura que considere a distribuição da pressão não-hidrostática no fluxo 2D. Este artigo deduz as equações de Boussinesq como uma evolução das equações de águas rasas com a hipótese de distribuição não hidrostática de pressão na direção vertical. Utiliza-se o sistema cartesiano ortogonal XYZ, considerando a influência da inclinação do fundo do canal e das perdas de carga no escoamento. É apresentada a correção não-hidrostática na equação de Boussinesq em duas dimensões usando as séries de Fourier. A solução utiliza o método Runge-Kutta Galerkin Descontínuo, e a formulação é aplicada ao rompimento de uma barragem cilíndrica.

Palavras-chave:
Equações de Boussinesq; Escoamento em canais; Modelagem matemática; Distribuição não-hidrostática de pressões; Modelo bidimensional

INTRODUCTION

The study of non-steady flow in channels is very important, considering the possible effects of floods or transient waves on population downstream. This flow is calculated from Navier-Stokes equations, which are a three-dimensional representation of the relationship between pressure and velocities at each point in space and time. The Navier-Stokes equations also do not allow a general solution for any type of flow. For practical applications, these equations must be approached through the Reynolds Averaged Navier-Stokes Equations (RANS), with the use of some closing equations for turbulence. The presentation of these equations is not part of this article but can be viewed for instance in Schilichting (1979) and Lai, Rubin and Krempl (1993).

These additional equations, associated with the greater number of discretization points in three-dimensional modeling, greatly increase the size of the systems of equations to be solved by computations, which makes the analysis of engineering works more complex or even impossible.

Models to simplify the flows of one or even two main dimensions of flow are often used to solve this problem.

This work presents the deduction of these simplified forms, called 2D Boussinesq Equations, an evolution from the Shallow Water Equations, and the non-hydrostatic pressure distribution is focused upon. Thus, we can study flows in a simplified way in two dimensions but considering possible effects of non-hydrostatic pressure distribution.

The influence of the non-hydrostatic distribution of pressures in a flow versus the traditional consideration of hydrostatic pressure distribution is analyzed.

This formulation is applied to a case study of rupture of a hypothetical cylindrical dam. The discretization of the correction term of non-hydrostatic pressure distribution is presented, and the influence of each part analyzed.

The generic name of Boussinesq equations include many partial differential equations with a common characteristic: they are treated phenomena that traditionally are simplified to linear, without simplification. We can find equations with this name in several fields of knowledge: in channel hydraulics or marine hydraulics – our case – in groundwater flow, atmospheric boundary layer, celestial mechanics, to stay in some areas of knowledge. Traditionally in Channel Hydraulics and Marine Hydraulics there are Boussinesq equations with consideration of non-hydrostatic pressure distribution, but do not consider the roughness of the bottom nor its slope (deep-water models). It is also observed that these models, because they are of maritime hydraulics, use spherical coordinates, not Cartesian coordinates.

LITERATURE REVIEW – UNSTEADY FLOW EQUATIONS

The Saint–Venant equations approximate the Navier –Stokes Equations to a one dimension where only the longitudinal variations of flow are represented, neglecting the variations of velocity and pressure in the transversal direction, or even the vertical one. They can be deduced from the Navier-Stokes Equation and were first presented by Adémas Jean-Claude Barrè, Count of Saint-Venant, in 1871. One of the many ways of writing the Saint-Venant equations is shown in Equation 1. These equations can also be called Shallow Water Equations in one dimension.

$\left\{\begin{array}{c}\frac{\partial h}{\partial t}+\frac{\partial }{\partial x}\left(uh\right)=0\\ \frac{\partial }{\partial t}\left(uh\right)+\frac{\partial }{\partial x}\left\{u^{2}h+g\frac{h^2}{2}\right\}=gh\left(S_0-S_f\right)\end{array}\right.$ (1)

Where: $h$: depth of flow (m); $u$: water velocity in $x$ direction (m s^-1); $g$: acceleration of gravity (9.81 m s^-2); $S_0$: bottom slope (m/m); e $S_f$: energy line slope (m/m), calculated, for instance, by the Manning equation (Equation 2):

$u=\frac{1}{n}{R}_h^{\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}{S}_f^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$ (2)

Where: $n$ – Manning coefficient (m^-1/3s) and $R_h$ is the hydraulic radius of the section (m) – calculated by the division of the area of flow by the wetted perimeter. Conserving the signal (i.e., the direction) of flow, we can describe Equation 3:

$S_f=\frac{u^2{n}^2}{R_h^{\raisebox{1ex}{$4$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}}=\frac{u^2{n}^2}{h^{\raisebox{1ex}{$4$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}}=\frac{u\left|u\right|n^2}{h^{\raisebox{1ex}{$4$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}}$ (3)

Equation 3 also considers a very wide channel, where it can be accepted that the hydraulic radius of the section is equal to the depth of flow.

In the case of the Shallow Water Equations – usually the two-dimensional equations are given this name – where we have the expansion of the Saint-Venant equations to two horizontal dimensions, variations in the transversal direction are considered, but not velocities or vertical accelerations. Equation 4, presented by Chaudhry (2008), shows this formulation.

$\left\{\begin{array}{c}\frac{\partial h}{\partial t}+\frac{\partial }{\partial x}\left(uh\right)+\frac{\partial }{\partial y}\left(vh\right)=0\\ \frac{\partial }{\partial t}\left(uh\right)+\frac{\partial }{\partial x}\left(u^{2}h+g\frac{h^2}{2}\right)+\frac{\partial }{\partial y}\left(vuh\right)=gh\left(S_0{}_x-S_f{}_x\right)\\ \frac{\partial }{\partial t}\left(vh\right)+\frac{\partial }{\partial x}\left(uvh\right)+\frac{\partial }{\partial y}\left(v^{2}h+g\frac{h^2}{2}\right)=gh\left(S_0{}_y-S_f{}_y\right)\end{array}\right.$ (4)

Where: $v$: flow velocity in $y$ direction (m s^-1); $S_0{}_x$and $S_0{}_y$: bottom slopes in x and y directions (m/m); and $S_f{}_x$and $S_f{}_y$: energy line slopes in the same directions (m/m) calculated, for instance, by the Manning equation – Equation 5 and 6, using total velocity V (m s^-1) given by Equation 7:

$S_f{}_x=\frac{u\left|V\right|n^2}{h^{\raisebox{1ex}{$4$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}}$ (5)

$S_f{}_y=\frac{v\left|V\right|n^2}{h^{\raisebox{1ex}{$4$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}}$ (6)

$V=\sqrt{u^2+v^2}$ (7)

The one- or two-dimensional Boussinesq equations are different from the equations mentioned above – Saint-Venant and Shallow Water – because they consider a variation of flow velocity in the vertical direction (adopted linear), which causes the acceleration of flow in this direction. This invalidates the hypothesis of hydrostatic distribution of pressures.

The introduction of corrective terms enables a more precise analysis without calibration parameters. Chaudhry (2008), in his chapter 12, presents the deduction of these equations for the one-dimensional case, here reproduced in Equation 8.

$$\left\{\begin{array}{l}\frac{\partial h}{\partial t}+\frac{\partial }{\partial x}\left(uh\right)=0\hfill \\ \begin{array}{r}\hfill \frac{\partial }{\partial x}\left\{u^{2}h+g\frac{h^2}{2}-\underset{Boussinesq’scorrectionfactor=B}{\underbrace{\frac{h^3}{3}\left[\frac{{\partial }^{2}u}{\partial x\partial t}+u\frac{{\partial }^{2}u}{\partial x^2}-{\left(\frac{\partial u}{\partial x}\right)}^2\right]}}\right\}+\dots \\ \hfill \dots +\frac{\partial }{\partial t}\left(uh\right)=gh\left(S_0-S_f\right)\end{array}\hfill \end{array}\right.$$ (8)

Equation 8 is mentioned in Mahmood and Yevjevich (1975, p. 62) as a higher order formulation of the Saint-Venant equations – which means that $w$ (velocity in the vertical direction (m s^-1)) is no longer null as assumed in the Saint-Venant Equations.

DEDUCTION OF BOUSSINESQ EQUATIONS IN TWO DIMENSIONS IN THE HORIZONTAL PLANE

According to what has been presented in the previous item, we can expand the Boussinesq assumption – existence of vertical acceleration based on the curvature of the current lines causing a non-hydrostatic pressure distribution – to two horizontal dimensions, creating a change like that previously presented for the Saint-Venant Equations also for the Shallow Water Equations.

The study of flows in two spatial dimensions is traditionally carried out by the equations called the Shallow Water Equations. These equations show the variation of water depth ($h$) and velocity components ($u$ and $v$) in the two spatial directions ($x$ and $y$). One of the basic hypotheses taken from the deduction of these equations is that the pressure distribution in vertical direction is hydrostatic. However, researchers have verified that in some flows the pressure distribution can result non-hydrostatic, due to downward flows generated by strong gradients. The studies presented in Roger et al. (2009, 2010), performed with the Shallow Water Equations, conclude that the observed difference between the values experimentally measured in the laboratory and those obtained in the computational modeling can be credited to the non-uniform velocity distribution and the non-hydrostatic pressure distribution.

We shall follow the path proposed by Chaudhry (2008), performing the necessary changes for two dimensions.

Figure 1 presents the flow characteristics in two dimensions used in the following derivation.

Figure 1
Presentation of the Characteristics used to Deduce the Two-Dimensional Boussinesq Equations.

Continuity equation

The Continuity equation in two horizontal dimensions and one vertical dimension (which will be simplified) can be written in the form of Equation 9:

$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0$ (9)

With $w$ as the velocity (m s^-1) in the $z$ direction. Multiplying by $dz$ and integrating over the entire depth $h$ – between $(z=z_b$ and $z=z_s=z_b+h)$ – using the Leibniz rule and imposing bottom boundary conditions – where ($z=z_b)$ results in $w_b=0$, because $w=0$ at the bottom (impervious bottom) – and on the free surface, where ($z=z_s)$ results in $p_s=0$ (atmospheric pressure), we can write the vertical velocity at the free surface as Equation 10:

$w_s=\frac{dh}{dt}=\frac{\partial h}{\partial t}+u\frac{\partial h}{\partial x}+v\frac{\partial h}{\partial y}$ (10)

Returning this equation in Equation 9 and integrating once again with the Leibniz rule, and recalling that $z_s=z_b+h$, Equation 11 is obtained:

$\frac{\partial h}{\partial t}+\frac{\partial }{\partial x}\left(uh\right)+\frac{\partial }{\partial y}\left(vh\right)=0$ (11)

This equation is identical to the Shallow Water equation, Equation 4a, as expected. Expanding Equation 11, and with a few additional mathematical treatments, Equation 12 is found, showing the equation for vertical velocity $w$:

$w=\frac{1}{h}\left(\frac{\partial h}{\partial t}+u\frac{\partial h}{\partial x}+v\frac{\partial h}{\partial y}\right)$ (12)

Conservation of the Momentum in $z$ direction

The equation of Momentum in $z$ direction may be written as Equation 13:

$\frac{\partial w}{\partial t}+u\frac{\partial w}{\partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial z}-g$ (13)

With: p is fluid pressure (Pa) and $\rho$ represents fluid density (kg m^-3).

Ignoring the head losses that occurs in the vertical flow and multiplying Equation 13 by $z$ and rearranging the terms, the result is that the term of pressure $\raisebox{1ex}{$p$}\!\left/ \!\raisebox{-1ex}{$\rho $}\right.$ can be written as in Equation 14:

$\frac{p}{\rho }=\frac{1}{\rho }\frac{\partial }{\partial z}\left(pz\right)+gz+\frac{\partial }{\partial t}\left(wz\right)+\frac{\partial }{\partial x}\left(uzw\right)+\frac{\partial }{\partial y}\left(vzw\right)+\frac{\partial }{\partial z}\left(w^{2}z\right)-w^2$ (14)

Integrating Equation 14 in the vertical direction $z$, Equation 15 is obtained. Equation 15 shows that if $w=0$ (without velocity in $z$ direction), the pressure distribution is hydrostatic and obtained directly, as shown by Equation 16, which is used in the Saint-Venant and Shallow Water equations.

${{\displaystyle \int }}_0^h\frac{p}{\rho }dz=\frac{\partial }{\partial t}{{\displaystyle \int }}_0^{h}wzdz+\frac{\partial }{\partial x}{{\displaystyle \int }}_0^{h}uwzdz+\frac{\partial }{\partial y}{{\displaystyle \int }}_0^{h}vwzdz-{{\displaystyle \int }}_0^h{w}^{2}dz+\frac{1}{2}g{h}^2$ (15)

${{\displaystyle \int }}_0^h\frac{p}{\rho }dz=g\frac{h^2}{2}$ (16)

Conservation of Momentum in $x$ direction

The equation of Momentum in $x$ direction may be written as shown in Equation 17:

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial x}+g{S}_{0_x}-g{S}_{f_x}$ (17)

Where term $g{S}_{0_x}$ represents the component of weight of fluid in the direction of movement and $g{S}_{f_x}$ represents the head losses in the $x$ direction (because of fluid viscosity and contour roughness). Based on Equation 9 integrated in the vertical and applying a few mathematical operations, Equation 18 is obtained:

$\left[\frac{\partial }{\partial t}\left(uh\right)+\frac{\partial }{\partial x}\left(u^{2}h\right)+\frac{\partial }{\partial y}\left(uvh\right)+{{\displaystyle \int }}_0^h\frac{1}{\rho }\frac{\partial p}{\partial x}dz\right]=gh\left(S_{0_x}-S_{f_x}\right)$ (18)

Returning to Equation 14, multiplying it by $dz$ and integrating it throughout its depth, Equation 19 is obtained:

${{\displaystyle \int }}_0^h\frac{p}{\rho }dz=\frac{\partial }{\partial t}{{\displaystyle \int }}_0^{h}wzdz+\frac{\partial }{\partial x}{{\displaystyle \int }}_0^h\left(wuz\right)dz+\frac{\partial }{\partial y}{{\displaystyle \int }}_0^h\left(wvz\right)dz-{{\displaystyle \int }}_0^h{w}^{2}dz+\frac{1}{2}g{h}^2$ (19)

Applying the Leibniz rule in the remaining integral in Equation 18 and rewriting the terms we have Equation 20:

$\frac{\partial }{\partial t}\left(uh\right)+\frac{\partial }{\partial x}\left(u^{2}h\right)+\frac{\partial }{\partial y}\left(uvh\right)+\frac{\partial }{\partial x}{{\displaystyle \int }}_0^h\frac{p}{\rho }dz=gh\left(S_{0_x}-S_{f_x}\right)$ (20)

since we have atmospheric pressure at the top of runoff. Substituting Equation 19 in Equation 20, the result is Equation 21:

$\begin{array}{l}\frac{\partial }{\partial t}\left(uh\right)+\frac{\partial }{\partial x}\left(u^{2}h+\frac{1}{2}g{h}^2\right)+\frac{\partial }{\partial y}\left(uvh\right)+\\ \frac{\partial }{\partial x}\left[\frac{\partial }{\partial t}{{\displaystyle \int }}_0^{h}wzdz+\frac{\partial }{\partial x}{{\displaystyle \int }}_0^{h}wuzdz+\frac{\partial }{\partial y}{{\displaystyle \int }}_0^{h}wvzdz-{{\displaystyle \int }}_0^h{w}^{2}dz\right]=gh\left(S_{0_x}-S_{f_x}\right)\end{array}$ (21)

The last term on the left side of Equation 21 can be determined by the sum of Equations 22 to 25. Focusing on this term (between brackets), there is:

$\begin{array}{l}\frac{\partial }{\partial t}{{\displaystyle \int }}_0^{h}wzdz=\frac{\partial }{\partial t}{{\displaystyle \int }}_0^h\left[-\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)\right]z^{2}dz=-\frac{\partial }{\partial t}\left[\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)\frac{h^3}{3}\right]=\\ -\frac{h^3}{3}\left(\frac{{\partial }^{2}u}{\partial x\partial t}+\frac{{\partial }^{2}v}{\partial y\partial t}\right)-\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)h^2\frac{\partial h}{\partial t}\end{array}$ (22)

$\begin{array}{l}\frac{\partial }{\partial x}{{\displaystyle \int }}_0^{h}wuzdz=\frac{\partial }{\partial x}{{\displaystyle \int }}_0^h\left[-\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)\right]z^{2}udz=-\frac{\partial }{\partial x}\left[\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)u\frac{h^3}{3}\right]=\\ -\left[u\frac{h^3}{3}\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y\partial x}\right)+\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)h^{2}u\frac{\partial h}{\partial x}+\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)\frac{h^3}{3}\frac{\partial u}{\partial x}\right]\end{array}$ (23)

$\begin{array}{l}\frac{\partial }{\partial y}{{\displaystyle \int }}_0^{h}wvzdz=\frac{\partial }{\partial y}{{\displaystyle \int }}_0^h\left[-\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)\right]z^{2}vdz=-\frac{\partial }{\partial y}\left[\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)v\frac{h^3}{3}\right]=\\ -\left[v\frac{h^3}{3}\left(\frac{{\partial }^{2}u}{\partial y\partial x}+\frac{{\partial }^{2}v}{\partial y^2}\right)+\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)h^{2}v\frac{\partial h}{\partial y}+\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)\frac{h^3}{3}\frac{\partial v}{\partial y}\right]\end{array}$ (24)

$\underset{0}{\overset{h}{{\displaystyle \int }}}{w}^{2}dz=\underset{0}{\overset{h}{{\displaystyle \int }}}{\left[-\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)\right]}^2{z}^{2}dz=\left[\frac{h^3}{3}{\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)}^2\right]$ (25)

The sum of the previous equations contains the continuity equation which is null, from which Equation 26 is obtained:

$\dots =-\frac{\partial }{\partial x}\left\{\frac{h^3}{3}\left\{\left(\frac{{\partial }^{2}u}{\partial x\partial t}+\frac{{\partial }^{2}v}{\partial y\partial t}\right)+\left[u\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y\partial x}\right)+v\left(\frac{{\partial }^{2}u}{\partial y\partial x}+\frac{{\partial }^{2}v}{\partial y^2}\right)\right]-{\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)}^2\right\}\right\}$ (26)

Designating as $B$ – Boussinesq Correction Coefficient – the function (shown in Equation 27) that is being derived in $x$ direction in Equation 26:

$B=\left\{\frac{h^3}{3}\left\{\left(\frac{{\partial }^{2}u}{\partial x\partial t}+\frac{{\partial }^{2}v}{\partial y\partial t}\right)+\left[u\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y\partial x}\right)+v\left(\frac{{\partial }^{2}u}{\partial y\partial x}+\frac{{\partial }^{2}v}{\partial y^2}\right)\right]-{\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)}^2\right\}\right\}$ (27)

In this way, Equation 21 is written in the form of Equation 28:

$\frac{\partial }{\partial t}\left(uh\right)+\frac{\partial }{\partial x}\left(u^{2}h+\frac{1}{2}g{h}^2\right)+\frac{\partial }{\partial y}\left(uvh\right)-\frac{\partial }{\partial x}B=gh\left(S_{0_x}-S_{f_x}\right)$ (28)

Conservation of Momentum in$y$ direction

The equation of Momentum in$y$ direction can be written as in Equation 29:

$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial y}+g{S}_{0_y}-g{S}_{f_y}$ (29)

In the same way as presented in item 3.3, when the variation in the momentum in $x$ direction was analyzed, in $y$ direction we have the exchange of velocity $u$ for velocity $v$ and the directions $x$ for $y$. Hence, following the same steps as the previous item, we can describe the correction term as shown in Equation 30:

$\dots =-\frac{\partial }{\partial y}\left\{\frac{h^3}{3}\left\{\left(\frac{{\partial }^{2}u}{\partial x\partial t}+\frac{{\partial }^{2}v}{\partial y\partial t}\right)+\left[u\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y\partial x}\right)+v\left(\frac{{\partial }^{2}u}{\partial y\partial x}+\frac{{\partial }^{2}v}{\partial y^2}\right)\right]-{\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)}^2\right\}\right\}$ (30)

Which is like Equation 26 with the change of the spatial derivate of term $B$. The Equation of Momentum in $y$ direction is like Equation 28, but with the appropriate adaptations due to the differences in direction shown in Equation 4b and Equation 4c, as shown in Equation 31.

$\frac{\partial }{\partial t}\left(uh\right)+\frac{\partial }{\partial x}\left(uvh\right)+\frac{\partial }{\partial y}\left(v^{2}h+\frac{1}{2}g{h}^2\right)-\frac{\partial }{\partial y}B=gh\left(S_{0_y}-S_{f_y}\right)$ (31)

The Boussinesq equations in two horizontal dimensions finally result in the form presented in Equation 32:

$\left\{\begin{array}{c}\frac{\partial h}{\partial t}+\frac{\partial }{\partial x}\left(uh\right)+\frac{\partial }{\partial y}\left(vh\right)=0\\ \frac{\partial }{\partial t}\left(uh\right)+\frac{\partial }{\partial x}\left(u^{2}h+g\frac{h^2}{2}-B\right)+\frac{\partial }{\partial y}\left(uvh\right)=\dots \\ \dots =gh\left(S_{0_y}-S_{f_y}\right)\\ \frac{\partial }{\partial t}\left(vh\right)+\frac{\partial }{\partial x}\left(uvh\right)+\frac{\partial }{\partial y}\left(v^{2}h+g\frac{h^2}{2}-B\right)=\dots \\ \dots =gh\left(S_{0_y}-S_{f_y}\right)\end{array}\right.$ (32)

The correction Boussinesq term is calculated by Equation 27.

MODELING THE EQUATIONS

To model Equation 32 the Runge-Kutta Discontinuous Galerkin method (RKDG), described in Schwanenberg (2003) and Schwanenberg and Harms (2004), was adopted. This method allows the interpolated functions to be discontinuous at the element boundaries, unlike the traditional Finite Element Methods in which the functions must be continuous between the elements.

The Runge-Kutta method is an explicit method in time, that is, the conditions in the future time are calculated as a function only of the known values in the previous instant. To ensure convergence in this kind of explicit method it is necessary to respect the Courant-Fiedrichs-Lewis condition, given generically by Equation 33, which must be calculated point by point.

$$CFL=\frac{V_{máx}}{\min \left(\frac{\Delta x}{\Delta t};\frac{\Delta y}{\Delta t}\right)}<1.0$$ (33)

In the case of this RKDG method, using a 2^nd order method, the condition is more restrictive. The upper limit becomes 0.35 instead of 1.0. Further details can be obtained from Schwanenberg (2003).

This method allows the resulting functions to be discontinuous between the elements, allowing the introduction of fluxes in the contours. In this way, the method allows the analysis of discontinuities such as, for example, abrupt wave fronts.

In the RKDG Method it is necessary to adopt a limiter to avoid strong numerical oscillations. There are several methods cited in the literature; Cockburn and Shu (2001) present a set of these limiters. In this work it is used the HLL method (Harten-Lax & van Leer), also used in the works presented by Schwanenberg (2003). This limiter is based in the method of the characteristics and seeks to limit the transfer of information to the future by the speed and celerity of the flow at the studied point. This type of approach seeks to avoid strong discontinuities - of numerical and non-physical origin - that can lead to simulation to present great instabilities. Schwanenberg (2003) presents the application of the Discontinuous Galerkin method to the Saint-Venant Equations (in one dimension) and Shallow Water (2D Equations). Such formulation is in the sequence, already adapted to contain the term referring to the correction of Boussinesq, object of this work.

The Shallow Water equations can be written in vector and conservative manner as Equation 34.

${\partial }_{t}U+{\partial }_{x}F\left(U\right)+{\partial }_{y}G\left(U\right)=S\left(U\right)$ (34)

Where each term is described by Equation 35 to Equation 38 – and ${\partial }_i=\raisebox{1ex}{$\partial \left(.\right)$}\!\left/ \!\raisebox{-1ex}{$\partial i$}\right.;i=x,y,t$.

$U=\left\{\begin{array}{c}h\\ uh\\ vh\end{array}\right\}$ (35)

$F\left(U\right)=\left\{\begin{array}{c}uh\\ u^{2}h+\frac{1}{2}g{h}^2\\ uvh\end{array}\right\}$ (36)

$G\left(U\right)=\left\{\begin{array}{c}vh\\ uvh\\ v^{2}h+\frac{1}{2}g{h}^2\end{array}\right\}$ (37)

$S\left(U\right)=\left\{\begin{array}{c}0\\ gh\left(S_{0_x}-S_{f_x}\right)+\frac{\partial B}{\partial x}\\ gh\left(S_{0_y}-S_{f_y}\right)+\frac{\partial B}{\partial y}\end{array}\right\}$ (38)

The corrective term B is given by Equation 27. The products of depth by velocity represent the flows in the directions$x$ and $y$ ($q_x=uh;q_y=vh$), respectively. The term of Boussinesq already appears here being moved to the source term, because it is small (since it is a correction).

The Finite Element Method in its solution multiplies the original equation by a weight function (ψ) and integrates it into each element K, as described by Equation 39, to minimize the approximation error between the calculated and observed values.

$\underset{K}{\overset{}{{\displaystyle \int }}}\left({\partial }_{t}U+{\partial }_{x}F\left(U\right)+{\partial }_{y}G\left(U\right)-S\left(U\right)\right)\psi dxdy=0$ (39)

In Equation 39 one must change the unknown $U$ by its approximation $U_h$, given by Equation 40, which uses the nodal values ($U^i$) and form functions (${\psi }_i\left(x\right)$), usually polynomials.

$U_h={\displaystyle \sum }_{i=1}^n{U}^i{\psi }_i\left(x\right)$ (40)

In the RKDG method the spatial and temporal discretizations can be totally separated, as this is an explicit process. In this way, we can write Equation 39 in the form of Equation 41 – where $\nabla \left(a,b\right)$ is a gradient tensor representation.

${\partial }_{t}U+\nabla \left(F\left(U,\nabla U\right),G\left(U,\nabla U\right)\right)=S\left(U\right)$ (41)

The unknowns appear explicitly and are the depth ($h$) and the streams in the directions $x$ and $y$ ($q_x$ and $q_x$).

The discretization in time by the second-order Runge-Kutta method is performed by applying Equation 42 to the system of equations Equation 4).

$U\left(t+\text{Δ}t\right)=U\left(t\right)+\frac{\text{Δ}t}{2}{L}_h\left(U\left(t\right),\frac{\partial U\left(t\right)}{\partial x},\frac{\partial U\left(t\right)}{\partial y}\right)$ (42)

Where $L_h\left(.\right)$ represents the linear operator that allows calculating the variation in time from the original variables and their spatial derivatives.

The modeling performed by Schwanenberg (2003) and Schwanenberg and Harms (2004) is also adopted in this work, using mesh with triangular finite elements.

DISCRETIZATION OF THE BOUSSINESQ TERMS

The term of Boussinesq was modeled as a source term, in order to be a term of correction for the original equations. The study of flows in two spatial dimensions is carried out traditionally by the equations denominated Shallow Water. These equations show the variation of levels and flows in each horizontal direction. In the case of the correction term of Boussinesq, its determination is made from the speeds and their derivatives in the spatial and temporal directions, besides the own depth. The values ​​referring to this element, its three direct neighbors and the neighbors of these elements, called surrounding ones, were selected for the calculation in each element. In this way, in each element, we will have 10 elements being considered for the calculations of the derivatives. To represent these functions the Fourier Series was adopted. Such a representation is interesting because the Fourier functions are well adapted to represent functions with strong gradients or even discontinuities, and always have derivatives of any order. We selected the Fourier function given by Equation 43, where we have sine and cosine in the two spatial directions, but only with one harmonic in each direction. The generic function φ represents the velocity components at $x$ and $y$ ($u$ and $v$) at the past (known) and future (calculation) instants.

$\phi =a_0+a_{1}cos\left(\frac{\pi x}{L_x}\right)+a_{2}sen\left(\frac{\pi x}{L_x}\right)+a_{3}cos\left(\frac{\pi y}{L_y}\right)+a_{4}sen\left(\frac{\pi y}{L_y}\right)$ (43)

Where $L_x$ and $L_y$ are the periodicity of the function, adopted as the maximum amplitude of the coordinates in the elements used for the calculation in the respective direction.

The use of a Fourier series is proposed in this work, since initial tests with polynomial functions did not converge. In the same way, Fourier series with more harmonics also did not converge (or more points were necessary).

From this formulation and adopting the least squares method to determine the coefficients for each approximation – coefficients $a_0$ to $a_4$ – we derive the derivatives contained in Equation 27, according to equations Equation 44 to Equation 48.

$\frac{\partial \phi }{\partial x}=-\frac{\pi }{L_x}{a}_{1}sen\left(\frac{\pi x}{L_x}\right)+\frac{\pi }{L_x}{a}_{2}cos\left(\frac{\pi x}{L_x}\right)$ (44)

$\frac{\partial \phi }{\partial y}=-\frac{\pi }{L_y}{a}_{3}sen\left(\frac{\pi y}{L_y}\right)+\frac{\pi }{L_y}{a}_{4}cos\left(\frac{\pi y}{L_y}\right)$ (45)

$\frac{{\partial }^2\phi }{\partial x^2}=-{\left(\frac{\pi }{L_x}\right)}^2{a}_{1}cos\left(\frac{\pi x}{L_x}\right)-{\left(\frac{\pi }{L_x}\right)}^2{a}_{2}sen\left(\frac{\pi x}{L_x}\right)$ (46)

$\frac{{\partial }^2\phi }{\partial y^2}=-{\left(\frac{\pi }{L_y}\right)}^2{a}_{3}cos\left(\frac{\pi y}{L_y}\right)-{\left(\frac{\pi }{L_y}\right)}^2{a}_{4}sen\left(\frac{\pi y}{L_y}\right)$ (47)

$\frac{{\partial }^2\phi }{\partial x\partial y}=0$ (48)

Figure 2 shows a generic element of the mesh, with a generic function represented in each centroid and the interpolation surface generated by the Least Squares method with the adopted Fourier Series.

Figure 2
Representation of a generic element, its neighbors and surroundings, and a generic function to be approximated (in vertical axis).

From Equation 27 and the derivatives calculated by Equation 44 to Equation 48, we can write the three portions of the Boussinesq correction, as in Equation 49 to Equation 51, with Equation 52 showing the form of the complete correction.

The neighboring and surrounding elements influence the derivatives, but the values are determined only for the element in question, which is in a position approximately central to the set.

$B_1=\left\{\frac{h^3}{3}\left[\left(\frac{{\partial }^{2}u}{\partial x\partial t}+\frac{{\partial }^{2}v}{\partial y\partial t}\right)\right]\right\}$ (49)

$B_2=\left\{\frac{h^3}{3}\left[u\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}v}{\partial x\partial y}\right)+v\left(\frac{{\partial }^{2}u}{\partial x\partial y}+\frac{{\partial }^{2}v}{\partial y^2}\right)\right]\right\}$ (50)

$B_3=-\left\{\frac{h^3}{3}\left[{\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)}^2\right]\right\}$ (51)

$B=B_4=B_1+B_2+B_3$ (52)

More details about deduction and applications were shown in Fabiani (2016).

APPLICATION TO A CYLINDRICAL DAM-BREAK

The method described above was applied to the case study of cylindrical dam-break proposed by Alcrudo and Garcia-Navarro (1993), and analyzed by many other researchers: Anastasiou and Chan (1997), Louaked and Hanich (1998), Mingham and Causon (1998), Ming and Chu (2000), Fujihara and Borthwick (2000), Schwanenberg (2003), Schwanenberg and Harms (2004), Gottardi and Venutelli (2004) and Fabiani (2016). It is a theoretical case, presented to test the behavior of numerical models, still without experimental results.

The case consists of a square tank with 40 m of side, filled with water at rest until a height of 1 m. Centered in this square there is a cylinder with 11 m radius, filled with water up to a height of 10 m. There is an instantaneous disruption of this cylinder, and the water profile is analyzed after 0.69 s of flow. Figure 3 graphically shows this condition, while Figure 4 shows the mesh adopted in the simulation, with 5,422 triangular elements.

Figure 3
Initial condition in cylindrical dam-break study.

Figure 4
Grid adopted in cylindrical dam-break simulation (with 5422 elements).

Applying the conventional Shallow Water Equations

The first simulation performed sought to reproduce with the proposed model the results presented by Toro (2001), which has a semi-analytical formulation. In Figure 5 the results obtained against the semi-analytical ones in the radial form (i.e., the distance to the center, $r$, and not the $x$ and $y$ coordinates) are plotted. In this way the results are coupled.

Figure 5
Results obtained using Shallow Water equation.

It is observed that the obtained results are consistent, presenting a damping of the wave fronts, characteristic of low order models (in this case, Runge-Kutta of 2^nd order).

Applying the Boussinesq Equations

The following application consisted of analyzing the influence of each of the three parts of the non-hydrostatic correction term. The flow was simulated with each of the plots given by Equation 49 to Equation 52, with the results being confronted with the results obtained with the Shallow Water Equations, shown in Figure 5. The selection in the software is performed by a flag called Bou, which varies from 0 to 4 - with zero being the conventional Shallow Water Equation, without the Boussinesq correction, and 4 represents the Fully Boussinesq correction. Figure 6 presents these results of the complete correction, in the radial format. A consideration of the presented results is that the inclusion of the non-hydrostatic term increases the dispersion of the responses - which is observed by the amplitude of the points around the result for Shallow Waters. This characteristic is observed mainly in the region near the wave front.

Figure 6
Results obtained with full Boussinesq correction (a) perspective; (b) radial plot.

The analysis of the influence of each part in the total correction can be done through transversal sections by the axes and, also, separating the parts, as shown in Figure 7 and Figure 8 (section by the x and y axes, respectively). Figure 9 shows a detail of the results on the positive wave front on the $x$-axis – highlighted in Figure 7 – to better illustrate the results obtained.

Figure 7
Results with Boussinesq correction – Section by $x$-axis.

Figure 8
Results with Boussinesq correction – Section by $y$-axis.

Figure 9
Detail of the result on the positive $x$-axis.

From the analysis of these figures it is highlighted that the correction of the non-hydrostatic term is generating an asymmetry in the flow for all plots, small in the third term, Equation 51, and large in the first term, Equation 49.

The computation of the Boussinesq correction has a high computational cost. Table 1 presents the processing time values for the case of the use of the conventional Shallow Water Equations, considering each of the correction terms and the complete set. An increase of at least six times in processing time is observed.

Table 1
Processing times and computational effort.

Simulation option	Processing time (s) a	Compu-tational Effort b
Shallow Water (Bou = 0)	21.04	1.00
First term of correction (Bou = 1)	172.93	8.22
Second term of correction (Bou = 2)	139.73	6.64
Third term of correction (Bou = 3)	140.30	6.67
Fully Boussinesq correction (Bou = 4)	175.51	8.34

• ^a
  Obtained in i7 computer, with 16 GB RAM memory;
• ^b
  The effort is defined by the relationship between times of simulations with Boussinesq correction and with the Shallow Water Equations.

The results presented were obtained by limiting the maximum variation at 10% - at each point and at each time step - of the original variation by the conventional Shallow Water Equations at the same instant. Such a limiter was set since the consideration of non-hydrostatic pressure distribution is not the main effect of the flow variations, but a correction to one of the terms. This procedure explains the fact that several points (see Figure 9) show the same result considering only the first term (Bou = 1) or all of them are equal (Bou = 4). This indicates that the first term - which presents cross-derivatives in time and space - has reached this limit at this point, usually at the points of the wave front. In average terms, the Boussinesq correction is well below this limit, with the average values observed in the simulations presented in Table 2.

Table 2
Mean corrections.

Simulation option	Mean variation (%)
Shallow Water (Bou = 0)	-
First term of correction (Bou = 1)	3.41
Second term of correction (Bou = 2)	2.35
Third term of correction (Bou = 3)	1.99
Fully Boussinesq correction (Bou = 4)	3.41

We can also compare the results with the values expected by the semi-analytic theory – Toro (2001) – and with the results obtained with the use of the Shallow Water Equations by other researchers; Figure 10 presents this analysis graphically (in radial view).

Figure 10
Comparison of the results of this study with those obtained by Toro (2001), Tseng and Chu (2000), Schwanenberg (2003) and Gottardi and Venutelli (2004) for t=0.69 s.

The results obtained with the proposed methodology are shown to be better than those obtained by Tseng and Chu (2000) and Gottardi and Venutelli (2004), approaching the best results obtained in the literature – by Schwanenberg (2003). The difference for the best results can be credited to the various calculation of derivatives from velocities, in time and space (including cross-derivatives).

CONCLUSIONS

The Boussinesq Equations presented here are an evolution of the Shallow Water Equations in two horizontal dimensions where it is possible to consider a non-hydrostatic pressure distribution. This variation in the pressure distribution is generated by vertical flow, which causes an acceleration to appear in this spatial direction. Also considered are the channel bottom slopes – which include the part of weight in the equations – and the head losses in the two horizontal directions ($x$ and $y$), while the head loss in vertical direction is ignored.

The application of the Runge-Kutta Discontinuous Galerkin method to the Boussinesq equations proved to be feasible, but with additional computational effort compared to the use of conventional Shallow Water Equations. This work is more complex (with more calculation steps) and has created a slight asymmetry. However, an advance in the sequence seems to be feasible through the following evidences: when the portion containing only the cross-derivative in time (Bou = 1) is included, the results tend to be more asymmetrical, compared to the results obtained using only spatial derivatives terms (Bou = 2 and Bou = 3); the results obtained with the complete analysis (Bou = 4) lead to results that are very similar, and even identical at several points, to those obtained with only the first term, due to the importance of this term and the limiters adopted in the solution method that guarantees the stability of the method.

It is proposed to continue the research with the use of other approximation functions to calculate spatial velocity derivatives and to allow a better representation of the correction terms due to non-hydrostatic pressure distribution. Other cases of non-permanent flows or strong gradients should be analyzed in future research, seeking to generalize the methodology.

Notation

b = sub-index indicating that the variable is taken at the bottom of the channel;

B = Boussinesq correction factor [m³ s^-2];

g = gravity acceleration [9.81 m s^-2];

h = water depth [m];

n = Manning roughness coefficient [m^-1/3 s];

p = fluid pressure [Pa, N m^-2];

$R_h$ = Hydraulic radius [m];

$s$ = sub-index indicating that the variable is taken at the free surface;

S₀ = bottom slope [dimensionless, or m/m];

S_f = energy line slope [dimensionless, or m/m];

$S_0{}_x$and $S_0{}_y$= bottom slopes, in x and y directions [dimensionless, or m/m];

$S_f{}_x$and $S_f{}_y$= energy line slopes, in x and y directions [dimensionless, or m/m];

$t$ = temporal dimension [s];

$u,v$ and $w$= velocities [m s^-1], respectively in$x,y$and $z$ directions;

V = total velocity [m s^-1];

$x,y$and $z$ = coordinated directions [m];

$\rho$ = fluid density [kg m^-3].

ACKNOWLEDGEMENTS

We thank Professor Dirk Schwanenberg for the section of his shallow water source code, used as the basis for this work.

REFERENCES

• ALCRUDO, F.; GARCIA-NAVARRO, P. A high-resolution godunov-type scheme in finite volumes for 2D shallow-water equations. International Journal for Numerical Methods in Fluids, v. 16, n. 6, p. 489-505, 1993. http://dx.doi.org/10.1002/fld.1650160604
  » http://dx.doi.org/10.1002/fld.1650160604
• ANASTASIOU, K.; CHAN, C. T. Solution of the 2D shallow water equations using the finite volume method on unstructured triangular meshes. International Journal for Numerical Methods in Fluids, v. 24, n. 11, p. 1225-1245, 1997. http://dx.doi.org/10.1002/(SICI)1097-0363(19970615)24:11<1225::AID-FLD540>3.0.CO;2-D
  » http://dx.doi.org/10.1002/(SICI)1097-0363(19970615)24:11<1225::AID-FLD540>3.0.CO;2-D
• CHAUDHRY, M. H . Open channel flow New York: Springer, 2008. 523 p.
• COCKBURN, B.; SHU, C.-W. Runge-Kutta discontinuous galerkin methods for convection-dominated problems. Journal of Scientific Computing, v. 16, n. 3, p. 173-261, 2001. http://dx.doi.org/10.1023/A:1012873910884
  » http://dx.doi.org/10.1023/A:1012873910884
• FABIANI, A. L. T. Estudo da Equação de Boussinesq em duas dimensões horizontais. 2016. 152 f. Tese (Doutorado em Engenharia) - Universidade Federal do Paraná, Curitiba, 2016.
• FUJIHARA, M.; BORTHWICK, A. G. L. Godunov-type solution of curvilinear shallow-water equations. Journal of Hydraulic Engineering, v. 126, n. 11, p. 827-836, 2000. http://dx.doi.org/10.1061/(ASCE)0733-9429(2000)126:11(827)
  » http://dx.doi.org/10.1061/(ASCE)0733-9429(2000)126:11(827)
• GOTTARDI, G.; VENUTELLI, M. Central scheme for two-dimensional dam-break flow simulation. Advances in Water Resources, v. 27, n. 3, p. 259-268, 2004. http://dx.doi.org/10.1016/j.advwatres.2003.12.006
  » http://dx.doi.org/10.1016/j.advwatres.2003.12.006
• LAI, W. M.; RUBIN, D.; KREMPL, E. Introduction to continuum mechanics 3rd ed. Butterworth & Heinemann: Exeter, 1993. 556 p.
• LOUAKED, M.; HANICH, L. TVD scheme for the shallow water equations. Journal of Hydraulic Research, v. 36, n. 3, p. 363-378, 1998. http://dx.doi.org/10.1080/00221689809498624
  » http://dx.doi.org/10.1080/00221689809498624
• MAHMOOD, K.; YEVJEVICH, V. Unsteady flow in open channels USA: Water Resources Pub., 1975. 923 p. 3 v.
• MING, H. T.; CHU, C. R. Two-dimensional shallow water flows simulation using TVD-MacCormack Scheme. Journal of Hydraulic Research, v. 38, n. 2, p. 123-131, 2000. http://dx.doi.org/10.1080/00221680009498347
  » http://dx.doi.org/10.1080/00221680009498347
• MINGHAM, C. G.; CAUSON, D. M. High-resolution finite-volume method for shallow water flows. Journal of Hydraulic Engineering, v. 124, n. 6, p. 605-614, 1998. http://dx.doi.org/10.1061/(ASCE)0733-9429(1998)124:6(605)
  » http://dx.doi.org/10.1061/(ASCE)0733-9429(1998)124:6(605)
• ROGER, S.; DEWALS, B. J.; ERPICUM, S.; SCHWANENBERG, D.; SCHÜTTRUMPF, H.; KÖNGETER, J.; PIROTTON, M. Experimental and numerical investigations of dike-break induced flows. Journal of Hydraulic Research, v. 47, n. 3, p. 349-359, 2009. http://dx.doi.org/10.1080/00221686.2009.9522006
  » http://dx.doi.org/10.1080/00221686.2009.9522006
• ROGER, S.; KÖNGETER, J.; SCHÜTTRUMPF, H.; ERPICUM, S.; ARCHAMBEAU, P.; PIROTTON, M.; SCHWANENBERG, D.; DEWALDS, B. J. Hybrid modeling of dike-break induced flows. In: RIVER FLOW 2010: INTERNATIONAL CONFERENCE ON FLUVIAL HYDRAULICS, 2010, Braunschweig. Proceedings… Karlsruhe: Bundesanstalt für Wasserbau, 2010. p. 523-531.
• SCHILICHTING, H. Boundary layer theory. 7th ed. USA: McGrawHill, 1979. 817 p.
• SCHWANENBERG, D. Die Runge-Kutta-Discontinuous-Galerkin-Methode zur Lösung konvections-dominierter tiefengemittelter Flachwasserprobleme 2003. 161 f. Dissertation (Doktors der Ingenieurwissenschaften) – Aachen University, Aachen, 2003.
• SCHWANENBERG, D.; HARMS, M. Discontinuous galerkin finite-element method for transcritical two-dimensional shallow water flows. Journal of Hydraulic Engineering, v. 130, n. 5, p. 412-421, 2004. http://dx.doi.org/10.1061/(ASCE)0733-9429(2004)130:5(412)
  » http://dx.doi.org/10.1061/(ASCE)0733-9429(2004)130:5(412)
• TORO, E. F. Shock-capturing methods for free-surface shallow flows Wiltshire: John Wiley and Sons, 2001. 309 p.
• TSENG, M.-H.; CHU, C. Two-dimensional shallow water flows simulation using TVD-MacCormack Scheme. Journal of Hydraulic Research, v. 38, n. 2, p. 123-131, 2000. http://dx.doi.org/10.1080/00221680009498347
  » http://dx.doi.org/10.1080/00221680009498347

Publication Dates

• Publication in this collection
  25 Apr 2019
• Date of issue
  2019

History

• Received
  27 Sept 2018
• Reviewed
  24 Feb 2019
• Accepted
  23 Mar 2019