Ocean Modelling 27 (2009) 198–214
Contents lists available at ScienceDirect
Ocean Modelling journal homepage: www.elsevier.com/locate/ocemod
A depth-integrated model for weakly dispersive, turbulent, and rotational fluid flows Dae-Hong Kim, Patrick J. Lynett *, Scott A. Socolofsky Zachry Department of Civil Engineering, Texas A&M University
a r t i c l e
i n f o
Article history: Received 13 June 2008 Received in revised form 8 January 2009 Accepted 17 January 2009 Available online 30 January 2009 Keywords: Boussinesq equations Finite volume method Vorticity Turbulent flow
a b s t r a c t A set of weakly dispersive Boussinesq-type equations, derived to include viscosity and vorticity terms in a physically consistent manner, is presented in conservative form. The model includes the approximate effects of bottom-induced turbulence, in a depth-integrated sense, as a second-order correction. Associated with this turbulence, vertical and horizontal rotational effects are captured. While the turbulence and horizontal vorticity models are simplified, a model with known physical limitations has been derived that includes the quadratic bottom friction term commonly added in an ad hoc manner to the inviscid equations. An interesting result of this derivation is that one should take care when adding such ad hoc models; it is clear from this exercise that (1) it is not necessary to do so – the terms can be included through a consistent derivation from the viscous primitive equations – and (2) one cannot properly add the quadratic bottom friction term without also adding a number of additional terms in the integrated governing equations. To solve these equations numerically, a highly accurate and stable model is developed. The numerical method uses a fourth-order MUSCL-TVD scheme to solve the leading order (shallow water) terms. For the dispersive terms, a cell averaged finite volume method is implemented. To verify the derived equations and the numerical model, four cases of verifications are given. First, solitary wave propagation is examined as a basic, yet fundamental, test of the models ability to predict dispersive and nonlinear wave propagation with minimal numerical error. Vertical velocity distributions of spatially uniform flows are compared with existing theory to investigate the effects of the newly included horizontal vorticity terms. Other test cases include comparisons with experiments that generate strong vorticity by the change of bottom bathymetry as well as by tidal jets through inlet structures. Very reasonable agreements are observed for the four cases, and the results provide some information as to the importance of dispersion and horizontal vorticity. Ó 2009 Elsevier Ltd. All rights reserved.
1. Introduction Boussinesq equation models are a popular choice for the simulation of weakly dispersive free-surface phenomena, such as wind waves in the nearshore area. Significant advances have been made in deriving the Boussinesq-type equations and in developing numerical models to solve them in recent years. Nearly all of this progress is based on the inviscid flow assumption; the Boussinesq equations are typically derived using the three-dimensional Euler or potential flow equations as a starting point. One of the common derivation methods is the perturbation approach, with the final equation model cast in depth-averaged form (e.g. Peregrine, 1967). Nwogu (1993) derived a new set of Boussinesq equations by using an arbitrary elevation za and extended the applicable water depth to the intermediate water regime. While Nwogu’s equation model is depth-integrated, the velocity variables are not * Corresponding author. E-mail address:
[email protected] (P.J. Lynett). 1463-5003/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ocemod.2009.01.005
in a depth-averaged form, and provide information about the varying vertical flow structure. Wei et al. (1995) extended Nwogu’s approach to capture nonlinear-dispersive effects, allowing for the simulation of very large amplitude, weakly dispersive waves. Here, we extend these models to incorporate friction effects by bottomgenerated turbulence through the common quadratic bottom stress. The effects represent a second-order correction to the leading order, shallow-water hydrodynamics, which when consistently applied, allows the model to capture vertical and horizontal vorticity. These model enhancements are important to simulate energy and constituent transport by large two-dimensional coherent structures in coastal flows. Inclusion of rotationality in the standard Boussinesq model has been the subject of some recent studies. As mentioned above, the Boussinesq model is usually founded on either potential flow or Euler’s equations. However, many derivations, including the Euler’s-based approach, enforce zero horizontal vorticity. From a strict physical standpoint, if the two horizontal vorticity components are zero, the third must be constant (or zero for practical
199
D.-H. Kim et al. / Ocean Modelling 27 (2009) 198–214
consideration). Thus, it is difficult to consider an equation model with zero horizontal vorticity, whether derived from Eulers or Navier–Stokes equations, as capable of predicting vertical vorticity. This physical inconsistency has been approached in a number of ways. In the derivation of Hsiao et al. (2002), which is based on Eulers equations, a number of z-dependent terms resulted in the final Boussinesq model. As the model is depth-integrated, these residual z-dependent terms are mathematically nonsense, implying there is not a unique solution to the Hsiao et al. model. However, it was found that by enforcing zero vertical vorticity, all the z-dependent terms would disappear, resulting in a solvable, irrotational model. Chen (2006), faced with a similar issue of z-dependent terms, used a different approach. He eliminated the z-dependent terms by double-integrating the Boussinesq model, providing a model that included vertical vorticity, although without explicitly including horizontal vorticity. Others have included horizontal vorticity directly, such as Musumeci et al. (2005), who derived a set of Boussinesq-type equations with horizontal vorticity, where the vorticity was solved with a separate transport equation. To include the effect of the turbulent mixing in the Boussinesq equation model, Chen et al. (1999) used a quadratic bottom friction dissipation term with a current-based subgrid dissipation model. Both terms are added in an ad hoc manner to the derived, inviscid Boussinesq equations. Hinterberger et al. (2007) presented a ‘‘depth-averaged Large Eddy Simulation (LES)” closure, used within a shallow water wave equation model. They showed that the depth-averaged LES is considerably more economic and can produce results that are generally of sufficient accuracy for practical purposes, compared with the full three-dimensional LES, in shallow flows. One of the most important characteristics that a numerical model must have is stability, because, in many cases, numerical models are to be applied to real or complex conditions. Recently, approximate Riemann solvers have been widely used with finite volume methods to provide a stable and accurate solution for the analysis of the Euler equations and the shallow water equations. For the Boussinesq equations, Erduran et al. (2005) used the fourth-order MUSCL-TVD scheme and used an approximate Riemann solver to solve leading order terms in one-dimensional space. They used a finite difference scheme to solve the dispersive terms in their numerical model. The object of this paper is to present Boussinesq-type equations that include the depth-integrated viscosity and associated horizontal and vertical vorticity terms. Hence, the derivation starts from the Navier–Stokes equations. In order to make the numerical model stable, the Boussinesq equations are converted to the conservative form and a fourth-order MUSCL-TVD (Monotone Upstream-centered Scheme for Conservation Laws-Total Variation Diminishing) scheme is used to solve the leading order terms. For the dispersive terms, a cell averaged finite volume method is used. To verify the derived equations and the implemented numerical model, four tests are conducted. First, a solitary wave is simulated. Second, the velocity of uniform and steady flow is compared with analytical theory. The other cases center on comparisons with experiments that generate strong vortices and eddies by the change of the bottom bathymetry and by a tidal jet through an inlet.
derivation of the approximate, depth-integrated model, a nondimensionalization, or scaling, of the primitive equations is the first step. Consistent with previous Boussinesq-type approaches, it is expected that the leading order solution will be shallow water, and thus a long wave scaling is used. A spatial region is characterhorizontal length scale ‘o , wave ized by a typical water depth ho , a p ffiffiffiffiffiffiffiffi amplitude ao , and a time scale ‘o = gho . With these variables, the following dimensionless variables and a parameter can be introduced:
pffiffiffiffiffiffiffiffi 0 ðx0 ; y0 Þ z0 t0 gho h f0 ; z¼ ; t¼ ; h¼ ; f¼ ; ‘o ‘ ho ho ho o0 0 U ;V W0 p0 ho ðU; V Þ ¼ pffiffiffiffiffiffiffiffi ; W ¼ pffiffiffiffiffiffiffiffi ; p ¼ ; l¼ qgho ‘o gho l gho
ðx; yÞ ¼
ð1Þ
where ðx0 ; y0 Þ denotes horizontal axes, z0 is a vertical axis, t 0 is time, 0 h is water depth, f0 is water surface elevation, ðU 0 ; V 0 Þ are horizontal velocities, W 0 is the vertical direction velocity, and p0 is pressure. The g and q are a gravitational acceleration and density, respectively. All these variables are dimensional. The l is a standard parameter for a scale analysis of long waves. For this study, due to the depth integration and resulting loss of flow details in the vertical plane, it will be reasonable to divide the turbulent eddy viscosity into horizontal and vertical components, as is commonly done for shallow mixing studies. The Smagorinsky 0 model (1963) will be used for the horizontal eddy viscosity mht , that qffiffiffiffiffiffiffiffiffiffiffiffi 2 0 0 0 h0 is, mt ¼ ðC s D0Þ 2Sij Sij , where C s is a constant, the Sij is a strain rate tensor and D0 is the grid size. By applying the above scalings to the 0 horizontal eddy viscosity, mht can be expressed as h0 t
m ¼
C 2s D2 ho
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 qffiffiffiffiffiffiffiffi 2 @U @U @W gho þ 2l 2 þ 2l2 þ @z @x @z
ð2Þ
Eq. (2) is rewritten in the compact form 0
qffiffiffiffiffiffiffiffi
mht ¼ aho gho mht ;
ð3Þ
where a ¼ C 2s D2 . For the vertical eddy viscosity, we presume a shallow flow formulation, where the vertical turbulence is driven by the 0 bottom shear only. The expression mvt ¼ C h H0u0 is used where H0 is 0 the total water depth and u is the friction velocity. The vertical eddy viscosity can be non-dimensionalized as 0
qffiffiffiffiffiffiffiffi
qffiffiffiffiffiffiffiffi
mvt ¼ bho gho Hub ¼ bho gho mvt ;
ð4Þ
where ub represents a near-bottom free stream velocity such that u ¼ C ub and b is equal to C h C . The typical magnitude of b is Oð0:01Þ as C h Oð0:1Þ by Elder’s model (1959) and C Oð0:1Þ when Reynolds number is 5000–10,000. Finally, the continuity equation and the Navier–Stokes equations can be scaled with Eqs. (1), (3) and (4):
rUþ
@W ¼ 0; @z
ð5Þ
@U @U b @ @U mvt þ U rU þ W þ rp ¼ alr mht rU þ ; @t @z @z l @z @W @W @p l2 þ l2 U rW þ l2 W þ þ1 @t @z @z @ @W mv ¼ al3 r mht rW þ bl @z t @z
ð6Þ
ð7Þ
2. Boussinesq equations with viscosity terms 2.2. Derivation of the depth-integrated momentum equations 2.1. Dimensionless governing equations The basic approach for including viscous effects into the Boussinesq equations is to derive the governing equations not from Eulers equations but from the Navier–Stokes equations. For the
This derivation will be of the perturbation type, and a small parameter assumption must be made. Looking to the vertical momentum equation (7), it is assumed that Oðl2 Þ ¼ OðlbÞ 1, yielding
200
D.-H. Kim et al. / Ocean Modelling 27 (2009) 198–214
@p þ 1 ¼ Oðl2 ; lbÞ: @z
ð8Þ
The above indicates that to leading order, the pressure is hydrostatic, which will permit the standard depth integration to obtain a long wave model. Thus the derived model will be restricted to weakly dispersive waves and flow with weak vertical turbulence and rotation. This step provides significant physical insight into this class of problem, and indicates that in order for the flow to be assumed hydrostatic to leading order, both dispersive and turbulent effects must be weak. Any model that assumes hydrostatic pressure implicity assumes weak turbulence. Typically, the perturbation of the inviscid primitive equations is performed using l2 as the small parameter. In these inviscid cases, where of course a ¼ b ¼ 0, the small parameter choice essentially required by (8) is clear. This would be the choice when deriving the typical (inviscid) shallow water or Boussinesq-type equations. It is noted that the ‘‘true” Boussinesq equations are derived assuming a balance between nonlinearity, or wave amplitude to depth ratio, and frequency dispersion, l2 , where both effects are considered to be small. The weak nonlinearity assumption is often violated by nearshore wind waves, and can be discarded from the derivation (e.g. Wei et al., 1995). These new fully nonlinear equations are still referred to as Boussinesq or Boussinesq-type equations, despite the fact that they no longer employ the scaling assumptions associated with their namesake. Back to the scaling found here, with viscosity, the choice for the expansion parameter is not clear, as either l2 or lb could be used as the small parameter. Mathematically, there is no reason to choose one over the other, as in fact both would result in the same final dimensional equations. For the derivation presentation, l2 will be used, and this issue of ambiguity will be addressed later. Physical values are expanded with power series following:
f ¼
N X
lð2nÞ fn ;
ð9Þ
n¼1
where f ¼ p; U; V; W and l2 assumed to be small. Substituting this expansion into (7) or (8) gives po as hydrostatic. It follows that rpo is independent of z. This implies that in the horizontal momentum equation, all the other leading order should also be z-independent functions (Dellar and Salmon, 2005). Consequently, U o becomes U o ðx; y; tÞ. At the water surface and at the bottom, the following boundary conditions W f ¼ @f=@t þ U f rf at z ¼ f and W h þ U h rh ¼ 0 at z ¼ h can be applied. The vertical velocity can be expressed with the horizontal velocity terms by integrating the continuity equation, yielding
W o ¼ zS T;
ð10Þ
where S ¼ r U o and T ¼ r ðhU o Þ. With the perturbation analysis, the horizontal vorticity is expressed as
@U 0 @U 1 co 0 2 co r W ¼ l r W þ O l4 o @z0 ho @z ho co co ¼ l2 x1 þ O l4 ; ho ho
ð11Þ
pffiffiffiffiffiffiffiffi where co ¼ gho . A vertical profile of U 1 can be derived from Eq. (11) through a vertical integration:
1 1 2 U 1 ¼ z2 rS zrT þ h rS hrT þ 2 2 þ O l2
Z
Z z 1 2 1 2 x1 dz z rS þ zrT h rS þ hrT þ l2 2 2 h þ l2 U 1 ðhÞ þ O l4
U ¼ U o l2
As this derivation will make use of Nwogu’s (1993) approach, the horizontal velocity is evaluated at an arbitrary elevation z ¼ za ,
Z za 1 2 1 2 x1 dz za rS þ za rT h rS þ hrT þ l2 2 2 h ð14Þ þ l2 U 1 ðhÞ þ O l4
U a ¼ U o l2
Subtracting Eq. (14) from Eq. (13), U can be expressed in terms of Ua.
U ¼ U a þ l2
1 2 za z2 rS þ ðza zÞrT þ l2 X þ O l4 ; 2
x1 dz þ U 1 ðhÞ
h
ð12Þ
such that the horizontal velocity, including up to second-order terms, becomes
ð15Þ
Rz
where X ¼ za x1 dz. For later use, the horizontal velocity can be ex / 4 r 2 pressed in which U r1 ¼ X and /as /U ¼ Ua þ l U1 þ U1 þ O l / U 1 ¼ U 1 ; V 1 is defined as
U /1 ¼
1 2 z z2 rS þ ðza zÞrT: 2 a
ð16Þ
The vertical profile of pressure is found through integration of the vertical momentum equation. Noting that the vertical distribution of mvt is independent of z as shown in Eq. (4), the pressure can be expressed as
@S 1 2 @T z f2 þ l2 ðz f Þ 2 @t @t 2 2 21 2 z f U o rS þ l ðz fÞU o rT þl 2 1 þ l2 f2 z2 S2 þ l2 ðf zÞTS þ O l4 ; al3 ; bl3 2
p ¼ f z þ l2
ð17Þ
The next step in deriving the horizontal depth-integrated momentum equation is to express each term of the horizontal momentum equations through U a . These terms, included to elucidate how vorticity and viscosity terms appear, become
@U @U a @ 1 2 @X ¼ za z2 rS þ ðza zÞrT þ l2 þ O l4 þ l2 @t @t 2 @t @t ð18Þ
1 2 U rU ¼ U a rU a þ l2 r U a za z2 rS þ ðza zÞrT 2 2 2 ð19Þ þ l rðU a XÞ þ l n þ O l4 @U ¼ l2 z2 SrS þ zT rS þ zSrT þ T rT þ W o x1 þ O l4 @z 1 @S @T 1 @S l2 r f þ l2 r z2 rp ¼ rf l2 r f2 2 @t @t 2 @t @T 1 þ l2 r z l2 r f2 U a rS l2 rðfU a rT Þ @t 2 1 1 þ l2 r f2 S2 þ l2 r z2 U a rS þ l2 rðzU a rT Þ 2 2 1 l2 r z2 S2 þ l2 rðfTSÞ l2 rðzTSÞ þ O l4 2 W
alr ðmht rUÞ ¼ alr mht rU a þ O al3
z
ð13Þ
b @ l @z
mvt
@U @z
¼ bl
@ mvt x1 blmvt rS þ O bl3 @z
In Eq. (19), n ¼ ðnx ; ny Þ is defined as
ð20Þ
ð21Þ ð22Þ ð23Þ
D.-H. Kim et al. / Ocean Modelling 27 (2009) 198–214
8