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

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.


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

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.
( ) ( ) ( ) Where: h : depth of flow (m); u: water velocity in x direction (m s -1 ); g: acceleration of gravity (9.81 m s -2 ); 0 S : bottom slope (m/m); e f S : energy line slope (m/m), calculated, for instance, by the Manning equation (Equation 2): (2) Where: n -Manning coefficient (m -1/3 s) and h R 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: (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) Where: v: flow velocity in y direction (m s -1 ); 0x S and 0 y S : bottom slopes in x and y directions (m/m); and f x S and f y S : 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: 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. 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. (2009Roger et al. ( , 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.

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: With w as the velocity (m s -1 ) in the z direction.Multiplying by dz and integrating over the entire depth h -between = (atmospheric pressure), we can write the vertical velocity at the free surface as Equation 10: Returning this equation in Equation 9and integrating once again with the Leibniz rule, and recalling that s b z z h = + , Equation 11 is obtained: 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 :

Conservation of the Momentum in z direction
The equation of Momentum in z direction may be written as Equation 13: With: p is fluid pressure (Pa) and ρ 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 p ρ can be written as in Equation 14: Integrating Equation 14 in the vertical direction z , Equation 15 is obtained.Equation 15shows 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.

Conservation of Momentum in x direction
The equation of Momentum in x direction may be written as shown in Equation 17: Where term x 0 gS represents the component of weight of fluid in the direction of movement and x f gS 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: Returning to Equation 14, multiplying it by dz and integrating it throughout its depth, Equation 19 is obtained: Applying the Leibniz rule in the remaining integral in Equation 18 and rewriting the terms we have Equation 20: since we have atmospheric pressure at the top of runoff.Substituting Equation 19 in Equation 20, the result is Equation 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: The sum of the previous equations contains the continuity equation which is null, from which Equation 26 is obtained: Designating as B -Boussinesq Correction Coefficientthe function (shown in Equation 27) that is being derived in x direction in Equation 26: In this way, Equation 21 is written in the form of Equation 28:

Conservation of Momentum in y direction
The equation of Momentum in y direction can be written as in Equation 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 Fabiani and Ota 5/10 x for y.Hence, following the same steps as the previous item, we can describe the correction term as shown in Equation 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.
The Boussinesq equations in two horizontal dimensions finally result in the form presented in Equation 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.

;
. min 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.
Where each term is described by Equation 35to Equation 38and ( ) . ; , , The corrective term B is given by Equation 27.The products of depth by velocity represent the flows in the directions x and y ( ; ), 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.
In Equation 39 one must change the unknown U by its approximation h U , given by Equation 40, which uses the nodal values ( i U ) and form functions ( ( ) ), usually polynomials.

( )
In the RKDG method the spatial and temporal discretizations can be totally separated, as this is an explicit process.In this way, Two-dimensional Boussinesq equations applied to channel flows: deducing and applying the equations 6/10 we can write Equation 39 in the form of Equation 41-where The unknowns appear explicitly and are the depth (h) and the streams in the directions x and y ( x q and x q ).The discretization in time by the second-order Runge-Kutta method is performed by applying Equation 42 to the system of equations Equation 4).
Where ( ) 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.
Where x L and y L 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 0 a to 4 a -we derive the derivatives contained in Equation 27, according to equations Equation 44 to Equation 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.
From Equation 27and the derivatives calculated by Equation 44to Equation 48, we can write the three portions of the Boussinesq correction, as in Equation 49to Equation 51, with Equation 52showing 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.
Figure 2. Representation of a generic element, its neighbors and surroundings, and a generic function to be approximated (in vertical axis).
Fabiani and Ota 7/10 More details about deduction and applications were shown in Fabiani (2016).
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.

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.
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 49to 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.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.
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.
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    Fabiani and Ota 9/10 (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.
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).
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.

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

Figure 3 .
Figure 3. Initial condition in cylindrical dam-break study.

Figure 5 .
Figure 5. Results obtained using Shallow Water equation.

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

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

Figure 9 .
Figure 9. Detail of the result on the positive x -axis.
Two-dimensional Boussinesq equations applied to channel flows: deducing and applying the equations 10/10 S 0 = bottom slope [dimensionless, or m/m]; S f = energy line slope [dimensionless, or m/m]; 0x S and 0 y S = bottom slopes, in x and y directions [dimensionless, or m/m]; f x S and f y S = 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]; ρ = fluid density [kg m -3 ].
, shows this formulation.