Baixe o app para aproveitar ainda mais
Prévia do material em texto
Elsevier Oceanography Series, 55 SHALLOW WATER HYDRO DY NAM I CS Mathematical Theory and Numerical Solution for a Two-dimensional System of Shallow Water Equations Tan Weiyan Nanjing Research Institute of Hydrology and Water Resources Nanjing 2 10024. China Water & Power Press Elsevier Beijing, China Amsterdam 1992 Responsible editors Yan Cunli. Tang Xinhua Published by Water 8. Power Press, Beijing The distribution of this book is being handled by the following publishers For the U. S . A. and Canada Elsevier Science Publishing Company, Inc. 52, Vanderbilt Avenue New York, NY 10017,U.S.A. For the People’s Republic of China Water &. Power Press 6 ,Sanlihe Road Beijing IO00.14, China For all remaining areas Elsevier Science Publishers P. 0. Box 2 1 I , 1000 AE Amsterdam, The Netherlands ISBN O-A11-98751-7(V01.55) ISBN 0 - 1 1 1 - 1 1 623 - 1 (Series ) Copyright 0 1992 Water 8. Power Press and Elsevier Science Publishers B. V. All rights reserved. No part of this publication may be reproduced, stored in a retrieval system or transmitted i n any form or by any means, electronic, mechanical, photo- copying, recording or otherwise, without the prior written permission of the publish- ers. Printed in Hong Kong 1 CHAPTER 1 BACKGROUND IN MECHANICS 1.1 PHYSICAL OEJECTS I . THE GEOMETRY OF WATER BODIES UNDER STUDY The geometry is characterized by : 1. A free surface. 2. A gentle bottom slope: the inclined angle a is such that tan a- sin a , and there is no abrupt change in underwater topography. 3. Shallow water : water depth h is much smaller than wave length or the char- acteristic length of the water body L , h<<L. Generally, it is required that h / L < lo-'. If the amplitude of h is of the same order of magnitude as depth itself, it is referred to as a very shallow flow. 4. The value of the horizontal space scale is often between 1 m and 1000 km. I I . PROPERTIES OF FUiIDS 1. Continuum. The physical and mechanical properties of the fluid should not approach infinity or suffer a jump at any isolated points. 2. Viscosity. For laminar flows, the fluid can be described approximately as a Newtonian fluid, in which molecular viscosity plays an important role. For turbulent flows, it is actually Enon-Newtonian fluid featured .+its turbulent viscosity. 3. Incompressibility. The density of a fluid element does not change during its motion. The effect on the flow of volumetric variation can be neglected. 4 . Homogeneity. The fluid, as a medium in mass transportation and heat con- duction, is well mixed, or the spatial distribution of its density has no influence on the flow. Density-driven flow and layered flow due to sediment, salinity, pollution, temperature, etc. , are not taken into consideration. Flow and transportation can be calculated separately, by taking the resulting water level and flow velocity as inputs of transport equations. Homogeneity together with incompressibility means that the density p is constant. We often set the value of the density of fresh water at 1000 kg/m3 and that of seawater at 1025 kg/m3. 5. Isotropy. Parameters of material properties, such as the viscosity coefficient ,u, do not vary with direction. 1II. FLOW BEHAVIOR 1. Unsteady (nonstationary) flow. A steady flow can be considered as the limit 2 of an unsteady flow under fixed external conditions as time t increases infinitely. We take an inertial coordinate system as a reference frame for differentiating the two types of flows, since a steady flow in that system would become unsteady in another noninertial system. In shallow-water flow computations, it is common practice to take a particular horizontal plane as a coordinate plane, ignoring the Earth's curva- ture. 2. Due to shallowness, a relatively uniform distribution of the horizontal veloc- ity over a vertical is obtained, so that depth-averaging can be justified. An intrinsi- cally three-dimensional flow can be simplified into a plane flow by integrating the horizontal velocity over a vertical to obtain a depth-averaged value, and by ignoring the effects of vertical velocities. 3. The flow is generally rotational, where production, transportation, diffusion and dissipation of vortices take place simultaneously and continuously. Large-scale flows exist in horizontal planes, whereas only small ones occur in vertical planes. Large vortices can be regarded as a type of secondary flow superposed to the main flow. 4. Temperature is often treated as a constant, since heat production due to fric- tional dissipation and heat transfer is negligible. If there is a difference in tempera- ture, we do not consider the variation in density, viscosity and thermal conductivity, so that the flow field is decoupled from the temperature field. They can be calculated separately. 5. Free surface elevation varies gradually with a small curvature, so that com- pared with gravitational acceleration, the acceleration in the vertical direction can be ignored. This may be formulated as the assumption that a hydrostatic (linear) pres- sure distribution over depth is a good approximation in regions of continuous flow. Special measures should be taken in cases of steep bed slopes and discontinuities such as hydraulic jumps. When vertical accelerations cannot be ignored, an assumption may occasionally be introduced that vertical velocity is zero at the bottom of a water body, varies lin- early, and reaches a maximum value at the water surface. The momentum equation thus obtained is the Boussinesq equation. 6. The effect of surface tension can be ignored. 7. The flow has a time scale between one second and several days. N . EXTERNAL FORCES 1. Gravity is the chief force initiating and governing the flow, and gives rise to the name of 'gravity wave'. Underwater topography mainly manifests its influence through gravitational effects. 2. Coriolis inertial force, due to the earth revolving around its own axis. 3. Tide-raising forces exerted by the moon and the sun. 4. Frictional forces between the flow and bed. Dissipation of mechanical energy due to molecular and turbulent viscosity can also be included in this term. With such a technique, the fluid can be viewed as inviscid, but such an assumption does not hold at discontinuities of the fluid flow. 5. Wind stress generated by the wind field over the water surface. 3 6. Pressure gradient force generated by the atmospheric-pressure field on the water surface. The first three are body forces, whose values are related to the water depth, while the last three surface forces depend on the horizontal area of a fluid column an- alyzed. 1 . 2 SYSTEM OF FLUID-DYNAMICS EQUATIONS FOR HOMOGENEOUS ISOTROPIC INCO- MPRESSIBLE FLOWS 1. FOUR FUNDAMENTAL PHYSICAL LAWS FOR FLUID F W W S In any inertial coordinate system hold the following fundamental physical laws : 1. Mass conservation law d -(m) = 0 d t (1. 2. 1) where m=mass of the fluid element under study. The equation derived from this law is usually called the continuity equation or equation of mass conservation. 2. Momentum conservation law (Newton's second law in dynamics) Conservation of linear momentum d dt - ( m V ) = F (1. 2. 2) where V=velocity vector, P=external force. The equation derived from this law is usually called the equation of motion or equation of momentum conservation. Conservation of angular momentum (1. 2. 3) d - ( ( m T X V ) = F X T dt where T = position vector of the given mass. Only on the microscopic scale is Eq. (1. 2. 3) independent of Eq. (1. 2. 2). For a macroscopic continuum, a reasonable assumption says that the directions of the mi- croscopic angular momenta are randomly distributed. It has been proved in the me- chanics of continuous media that due to Eq. (1. 2. 3) stress appears to be symmetric at any point, just as in static equilibrium. Conversely,if only stress is symmetric, conservation of angular momentum must be satisfied. At that time, Eq. (1. 2. 3) can be derived from Eq. (1. 2. 2 ) , and does not bring about new restrictions to the flow. The equation derived from the second type of the law is usually called the equation of vorticity , which is used chiefly in vortex analyses of incompressible fluid flows. 3. Energy conservation law (the first law of thermodynamics) d dQ dW -(I?) = - + - dt dt dt (1. 2. 4) 4 where E=energy of mass m , consisting of internal energy and mechanical energy. Internal energy denotes thermodynamic energy associated with microscopic random motions of molecules. Mechanical energy denotes both the kinetic energy and the po- tential energy associated with observable macroscopic motions. & = heat, both that exchanged with the environment and that produced by frictional dissipation and, if any, variation of its volume. Heat flux represents total transportation of microscopic energy. W=work done on the mass by external forces. The equation derived from this law is usually called the energy equation or equation of energy conservation. If the heating process is independent of mechanical movement, then the one overall equation of energy can be decomposed into two, related to mechanical energy and thermal energy respectively. The former, which is no more than a dot product of the linear momentum equation and the velocity vector, describes the variation of me- chanical energy due to nonequilibrium between the forces. The latter describes the change in internal energy due to heat transfer and fluid deformation. On the other hand, if the heating process makes a notable impact on the mechanical movement, then the equation of energy is a unified and independent relation that must be satis- fied. As the effect of temperature on fluid flow is assumed to be small, this equation will not be utilized hereafter. Here, we just recall that heating or cooling due to com- pression or expansion of the fluid is a reversible process, whereas viscosity dissipation (which always increases internal energy) and heat conduction are irreversible pro- cesses. 4. The second law of thermodynamics dQ d S - - > O T (1. 2. 5 ) where S = entropy, and T = absolute temperature. Lower-case symbols s and e are used for thermodynamic parameters of the fluid per unit mass, called specific entropy and specific internal energy. They characterize uniquely a thermodynamic state. Through the equation of state (one type of constitutive equation showing properties of the material) , e = e (s , u ) , where specific volume ZI = 1 /p, the two thermodynamic laws can be connected together. Here are some special thermodynamic processes. When d&= 0, it is adiabatic and we have dS>O, i. e. , entropy must be non-decreasing. When the equality dS= dQ/T holds, it is reversible ; conversely, for an irreversible process, dS>dQ/T. Un- der both adiabatic and reversibility conditions, we have dS = 0 and call the process isentropic. If only the adiabatic condition holds, a steady, irrotational flow is also isen tropic. We can make an analysis for inviscid compressible flows. In a smooth region of the flow, the equality symbol holds identically under both adiabatic and isentropy conditions. Then, to close our problem needs to use merely the energy conservation law (of course, depending on whether or not heating is coupled with motion) and the equation of state. If discontinuities appear in the flow, where production of entropy causes the isentropy condition to fail, we must additionally use the second law. These conclusions will be discussed in detail in Section 3. 3. 5 I I . TWO ALTERNATIVE VIEWPOINTS IN FUIID DYNAMICS ANALYSIS Two alternative approaches can be distinguished according to what kind of space coordinate system we adopt. 1. The Lagrangian viewpoint, also called material description of fluid flow. The object under study is the motion of a given fluid particle (or element) occupying a space region. The surface of the region moves at a local velocity, so that it is al- ways covered by the same material. Independent variables used in the formulation are time and the initial positions of all the particles (the latter are the marks of different particles). For the convenience of analysis, we often establish a reference coordinate system fixed at the fluid element and moving together with it. New space coordinates are called Lagrangian coordinates, which are functions of time and the initial posi- tions. 2. The Eulerian viewpoint, also called space description of fluid flow. The ob- ject in this approach is the flow at a given position or in a fixed region. Independent variables are time and physical space position. Due to the different roles played by space and time, they can be dealt with separately. Time is often viewed as a parame- ter-like independent variable, so that a flow can be described by the evolution of a transient flow field. The Eulerian viewpoint is perfectly suited to steady flows, because the variable 1 does not appear in governing equations. With its use, most unsteady flows can also easily be dealt with due to the adoption of a fixed space coordinate system. The La- grangian viewpoint is used chiefly in those one-dimensional problems with moving in- ternal boundaries, such as discontinuities and interfaces between two layers of fluids. A moving coordinate system fixed at such a point on the internal boundary would be convenient. In summary, the Lagrangian viewpoint is associated with the particle approach, and the Eulerian viewpoint with the field formulation. III. INTEGRAL FORMS OF GOVERNING EQUATIONS IN FUIID DYNAMICS FOR A FINITE CONTROL VOUIME The above physical laws are described from the Lagrangian terms, when a coordinate system is fixed at the fluid element with mass m. Now we reformulate them at instant t from the Eulerian terms, by selecting a fixed finite control volume in the flow re- gion. 1. Equation of continuity (1. 2. 6) where Q =domain of the control volume with element h; &'=surface of the control volume, whose area element is ds represented by a vector in the outward-normal di- rection. 6 2. Equation of motion (1.2. 7) where Fs=sum of surface-force vectors exerted on S ; F8=sum of body forces exert- ed on the fluid per unit mass. Note that formal similarity exists between the above two equations. The first term on the left-hand side expresses a time rate of increment of some physical quanti- ty (mass or momentum) in the control volume. The second term expresses a flux of that physical quantity flowing out through the surface of the control volume. The terms on the right-hand side express external influences. In the integral forms of conservation laws there appear both volume and area in- tegrals. It is common practice to use the Green's theorem in transforming an area in- tegral into a volume integral. Define a convex regular region as one whose boundary surfaces are composed of several pieces such that outward-normals to each piece form a continuous vector field. For any convex regular region or any region that can be decomposed into a finite set of such regions, the Green's theorem can be expressed by JD g d o = J n,Ads ( i = 1,2,3) (1. 2. 8) where A=a continuously differentiable function, vector or Cartesian tensor (cf. IV) in Q n,=x,-component of the unit vector n directed along an outward-normal to S. Forms of this theorem in common use are as follows: Gradient theorem Gauss' divergence theorem JQv do = v nds L integral theorem JQv x Vdm = n x Vds L where p=scalar function and operator and denoted by grad. (1. 2. 9 ) (1. 2.9a) (1. 2. 9b) V =gradient operator, also called space Hamiltonian In a rectangular (also called Cartesian or Descartes) a a a ax ag az coordinate system, V = i - + j - + k - ; i , j , k are unit basis vectors in coordi- natedirections. V may be viewed as a vector entering into operations formally. V V is a dot (scalar) product of Vand V , called divergence and also denoted by div V . V X V is a cross (vector) product of V and V , called vorticity and also denoted by rot V. If V = (uI , u2 ,u3 IT is in terms of the Cartesian coordinates ( x l , z2, z 3 ) , then 7 (1. 2. 10) In this case the above theorems can be reduced to relations between integrals over a plane region and its boundary curve, and still retain its name, Gauss' theo- rem. All the integral laws stated above can be considered as specific forms of the gen- eral transport equation (1.2.11) where J, is the density of a transported quantity, f is the flux density of that quanti- t y , and cp is its source density. By using Green's formula, the second term can be re- duced to a volume integral so that under some differentiability condition we have the differential form 3 + , . f = y at (1. 2. 12) N . BRIEF INTROWCTION TO CARTESIAN TENSOR A tensor defined in a Cartesian coordinate system is called a Cartesian tensor. A scalar is zeroth order (order-0) tensor. A vector is an order-1 tensor. In a three-di- mensional space, an order-2 tensor T is composed of 32 = 9 components { t,, } ( z , j = 1 ,2,3) . The definitions of scalar, vector and tensor are based on how their compo- nents change under a rotation of the coordinate system. A scalar does not change at all. Any set of three scalars forms a vector only if they are subject to the well-known coordinate rotation formula. An order-2 tensor can be regarded as a matrix such that, when coordinates ( r l , zz, r 3 ) are changed into ( G I , i z , &) under a rotational transformation, the relation between its old and new components { t,, } and { t,, } (z, j = 1 , 2 , 3 ) is given by t,, = ~mnB,mRB,8 (1. 2. 13) where B,, is the direction cosine of the angle made by the $,-axis with the r,-axis. For brevity of notation, the Einstein summation convention is adopted in the above equa- tion and also hereafter. If a pair of common indices appears in a certain term, then it - " a u z a3 a, axl az ax3 denotes a summation over the index, e. g. , 2 = - + - + - . The notation of tensor is similar to that of vector. A tensor is denoted either by an uppercase only or by a lowercase with subscripts. To the former, called a symbol notation, no coordinate system is assigned. It has the merit of conciseness, but we should follow some special rules in writing a tensor equation. The latter, named an index notation, can show clearly the order of tensor, the arrangement of its compo- nents, and the role of the coordinate system. When writing component equations for the convenience of numerical solution, there is no special requirement as each compo- 8 nent is a scalar. One reason for using the tensor as our tool is that if some physical law has been formulated in the form of a Cartesian tensor equation, then it must hold in any or- thogonal coordinate system. Partial derivatives of a Cartesian tensor constitute a new higher-order tensor. For example, those of a first-order tensor {ti} satisfy the definition of second-order tensor (1. 2. 1 4 ) An order-2 tensor whose components can be arranged to form a symmetric ma- trix is called a symmetric tensor. A coordinate system can be found to annul its non- diagonal elements. Such coordinates are called major axes, while its diagonal ele- ments are called major values. Suppose that under an arbitrary orthogonal transformation of rectangular coordi- nates, the values of all components of a given tensor do not change; then it is an isotropic tensor. There exists no order-1 isotropic tensor. Any order-2 isotropic ten- sor can be written in the form @$, , where p is a scalar, while 4, denotes the Kroneck- er delta s,, = O(2 # .7), 6, = l ( 2 = 3 ) (1. 2. 15) Moreover, a symmetric and isotropic order-4 tensor can be expressed in a form con- taining only two scalars : t r j k l = nsb]sU + p(ackdjL + 6'16jk) (1. 2. 16) Chief operations defined for order-2 tensors include : a = t,,s,, (denoted as LI = T: S ) T,k = S,,t,k(denOted aS R = S T) (1 ) scalar product of two tensors (2) inner product of two tensors ( 3 ) vector product of a tensor and a vector ( 4 ) tensor product of two vectors (diad) (1. 2. 17) (1. 2. 18) (1. 2. 19) (1.2. 20) u, = v,t,,(denoted as u = V T) t,, = u,v,(denoted as T = uv) V. DIFFERENTIAL FORMS OF GOVERNING EQUATIONS IN FLUID DYNAMICS Assume that: (i) the integral conservation laws hold for any bounded control volume selected within some region; (ii) there is no singularity in the solution for that region. Applying the integral forms of governing equations to an infinitesimal fluid element, and using the Gauss theorem for changing area integrals into volume integrals, we arrive at the desired differential forms, in which the left-hand sides are just the integrands under volume integral symbol (the right-hand sides are zero). Such an example has been given in Eq. (1. 2. 12) . The derivation is based on the theorem that if an integral of a given function over any control volume vanishes, then that function is identically equal to zero at all of its points of continuity. There- fore, when some derivative or nonhomogeneous term in a differential equation suffer a discontinuity or approaches infinity, we must still use the integral form. 9 1. Equation of continuity (1.2.21) Defining the material derivative (also called the substantial derivative) as D / D t = a/at + V V (note that in Eqs. (1. 2. 1 ) - (1. 2. 4 ) the derivative d / d t should be understood as D / L & ) , then the above equation can be written as b + p v * v = o Dt (1.2. 22) It should be noted that, when we use the integral conservation laws, the order of the two operations, material derivative and integral, generally cannot be inversed. For an incompressible fluid, Eq. (1. 2. 21) reduces to v.v = 0 (1. 2. 23) A velocity vector field satisfying this equation is called a solenoidal vector field. In- deed, we can use the equation as a definition of generalized incompressible flow, in- stead of constancy of density. 2. Equation of motion P E = Dv v u + pFf l (1. 2. 24) where u denotes a symmetric order-2 Cartesian tensor, called the stress tensor, acting on an infinitesimal fluid element. Its nine components {a,, } fully describe the stress behavior at a point. Its diagonal elements correspond to normal stresses, while the nondiagonal ones correspond to shear stresses. The right-hand side of the above equation represents the net external force exert- ed on the fluid element per unit volume. According to the d'Alember Principle, the left-hand side can be moved to the right-hand side, and considered as a part of the ex- ternal force, called the inertial force. As for the integral conservation laws, the re- sultant of the inertial forces should be exerted on the centre of the mass. Then the moving element can be viewed as if it were in equilibrium. The left-hand side of the above equation can be expanded to yield some different forms as follows: Eulerian form (also called convective form or nonconservative form) ] ~7 + P F B av P[X + (v v>v = v which can be applied to Cartesian coordinate systems only. Conservation form (1. 2. 25) (1. 2. 26) where VV denotes a diad having the associative property that A (BC) = ( A B ) C . Lamb-Gromeko or Bernoulli form 10 (1. 2. 27) The above forms are suitable for any fluid. For an incompressible fluid, the constant p can be moved outside the derivation symbol. As a stress tensor appears in the equation of motion, the system, Eqs. (1. 2. 21 ) and (1. 2. 24) , is open, and cannot be solved. It is thus necessary to introduce a constitutive equation, a relation between stress tensor and strain tensor. If that e- quation does not itself change under any orthogonaltransformation of the coordinate system, the fluid is isotropic. V I . CONSTITUTIVE EQUATIONS FOR FLUIDS In fluid dynamics, a constitutive equation expresses the relation between the stress tensor and the deformation-rate tensor of a fluid. We will now explain the ba- sic logics involved. 1. The stress tensors in incompressible and compressible fluids are somewhat dif- ferent, expecially the two pressures are radically different. For an incompressible fluid, the stress tensor can be decomposed into two parts u =- P I f z ( I . 2. 28) where p=scalar , I=identity matrix, and t=order-2 bias (or viscosity) stress ten- sor. The first term on the right-hand side is a diagonal matrix expressing the isotropic part of the stress tensor. The diagonal element p is determined by p = - tru/3 = - (uli + uZ2 + u3,>/3 (1. 2. 29) where tru=sum of diagonal elements of the matrix u , called trace. When the fluid is motionless, t = 0 and p denotes the static pressure. For a moving fluid, p is called the dynamic pressure, which is intrinsically a mechanical quantity which does not depend on other thermodynamic state variables. Hence, after the movement has stopped, stress tensor will approach that in static state. If the flow is compressible, p is generally not equal to the value given by the above equation. It is necessary to consider the influence of the additional viscosity stemming from dilatation of the fluid. Difference between the two values of p often is a linear function of the inflation rate of the volume. However, we usually adopt Stokes' hypothesis, which ignores this factor, leading to a Stokes' fluid. Another feature of the compressible flow is that p not only satisfies balance rela- tions among stresses, but also depends on other thermodynamic parameters, thus gaining its name of thermodynamic pressure. Since we neglect the variation of tem- perature, the dependency can be expressed by an equation of state f (p , p ) = 0. A flu- id for which p is only a function of p is called a barotropic fluid. Specifically, for an isentropic flow (adiabatic and reversible) of a perfect polytropic gas, the equation of state can be formulated as - = const (1. 2. 30) PY where y = C,/C,, is the ratio of specific heats under constant pressure C, and under constant volume C,,. 11 Under isentropy conditions, another thermodynamic parameter, c , related to p , is defined as (1. 2. 31) For perfect gases, if Eq. (1. 2. 30) is introduced into the above equation, we obtain c = &&I. As the right-hand side is the coefficient of order-1 term in a Taylor's series expansion of p , c denotes the speed of propagation of a small distur- bance, the speed of sound in gas dynamics. However, the application of Eq. ( I. 2. 31) is conditional. The condition of isentropy also requires that there is no strong shock wave (discontinuity) occurring in the flow. The second term on the right-hand side of Eq. (1. 2. 28) represents the non- isotropic part of the stress tensor. Its diagonal elements are called (viscosity) normal stress, while nondiagonal ones are (viscosity) shear stress. The total normal stress in any coordinate direction equals the sum of negative pressure - p and viscosity normal stress, while the total shear stress is just viscosity shear stress. The above-mentioned Stokes' hypothesis is equivalent to saying that the mean value of the viscosity normal stress vanishes. 2. In a Newtonian fluid, at any point, the components of bias stress tensor T are homogeneous linear functions of the components of a second-order tensor L , the gradient of the velocity, &/ax,, called deformation rate tensor. A velocity field determines not only the translation and rotation of fluid ele- ments viewed as rigid bodies, but also the relative motions between fluid particles. In other words, the velocity gradient determines the deformation rate, including dilata- tion and shearing. Obviously, the assumption of linearity between bias stress tensor and deformation rate tensor is a reasonable 3-D generalization of Newton's law in one space dimension, thus called generalized Newton law. Now shear stress is proportion- al to velocity gradient in a shear layer, whilst in the elastic-solid case stress is propor- tional to deformation, but not its rate. The deformation rate tensor L can be expressed as a sum of its symmetric part E and its asymmetric part D (an order-2 tensor called the rotation tensor, vorticity ten- sor or spin tensor) L = E + Q E = - ( L + LT) , D = - ( L - LT) 1 1 2 2 (1. 2. 32) (1. 2. 33) (1. 2. 34) (1. 2. 35) The relation between the symmetric term E and the velocity gradient tensor V Ti is written as 1 E = -[Vv 4- (VTi)T] 2 t. (1. 2. 36) 12 The diagonal elements of E denote the change rates of the relative extension or con- traction in the three coordinate directions respectively, and their sum (trace of matrix E) is just the divergence, V V , showing the inflation rate of the fluid element. A nondiagonal element e,, denotes the shear deformation rate of the angle made by a pair of coordinate axes (relative to the indices) divided by -2, or equivalently, a shear velocity in the x,-direction of two neighboring particles located at the x,-axis divided by -2. Among the components of E, a consistency condition must be satisfied, i. e. , the integrability condition that velocity can be obtained by integrating Eq. (1. 2. 36). In the 2-D case, it can be expressed as (1. 2. 37) As for the rotation tensor, Q, because of asymmetry there are only three inde- pendent components, which are equal to the components of vorticity V X V divided by 2, and can be used as a measure of the angular velocity of the fluid element rotat- ing as a rigid body. It should be noted that vorticity is a local property of the flow field, independent of the curvature of the streamlines. Indeed, the bias stress tensor T depends on the symmetric tensor E only, and not on the rotation tensor Q. For a Newtonian fluid a general relation holds ztj = t g j k l e k l (1. 2. 38) where T is an order-4 viscosity tensor whose elements are related to temperature on- ly, and not to u and E. As temperature is assumed here to be fixed, T is a constant tensor. In the 3-D case, T has 3'=81 elements. For an isotropic fluid, according to Eq. (1. 2. 1 6 ) , there are only two independent constants, denoted by A and ,u. By mathematical derivation, Eq. (1. 2. 38) can be expressed as z = 2 p ~ + a ( ~ - v ) z (1. 2. 39) so we have a,, = - pa,, t- 2w,, + k k k f f , ] (1. 2. 40) For incompressible fluids, as div V = V V = 0 (or e,, = 0 ) , and under the as- sumption that 3h+2p= 0 , (1. 2. 4 1 ) z = 2,uE - -(V V)Z = , u [VV + ( V V ) T ] - $(V V ) I c =- PI + 2hE (1. 2. 42) V . a = - V p + div(2,uE) (1. 2. 42a) ,u is called the dynamic viscosity. In the 1-D case, it is the proportional constant be- tween shear stress and shear deformation rate, and it can be interpreted as the mo- mentum of the viscous force exerted on a unit area. It is a macroscopic characteristic 2,u 3 13 of momentum exchanges stemming from molecular motions. Due to their common mechanism, there is a relation between viscosity ,u, heat conductivity k (featuring energy exchange) and mass diffusivity coefficient D (featuring material transporta- tion (1. 2. 43) The dimension of is MILT, while its unit is P = lg/(cm s ) = 1 dyne s/cm2 in the CGS system, or Pa s (or denoted by Pi) = 1 newton s/mz = 10 P in the SI u- nit system. The value ,u of water at 20 'C and 1 atmospheric pressure equals 0. 001 Pa s . When temperature varies only slightly, ,u can be viewed as a constant, and is correlated slightly with pressure. For isotropic fluids, a constant ,u can be applied to all coordinate directions in the 3-D case; otherwise, an order-2 viscosity tensor should be introduced. Moreover, we define kinetic viscosity by v=,u/p,with dimension L2/T and u- nit St (abbreviation of Stokes) = 1 cmz/s in the CGS system or m2/s( = l o 4 St) in the SI system. h is a parameter related to the inflation of volume. As stated before, the Stokes hypothesis states that the trace of the tensor z equals zero, so from Eq. (1. 2. 39), h and ,u must satisfy the condition 3h+2,u=O. The quantity I+2,u/3 is called the sec- ond viscosity, while p is the first viscosity. For a Stokes' fluid, the second viscosity is equal to zero, A= - 2p/3 (cf. Eq. ( 1. 2. 4 1 ) ) . Stokes' hypothesis is accurate e- nough in engineering applications, at least for gases and Newtonian fluids. For a Stokes' fluid, it can be proved that only the viscosity normal stress exists when inflation rate in the same direction is not equal to the mean value over those in the three coordinate directions, while viscosity shear stress is proportional to shear de- formation rate. Especially in the case of liquids, there are two sources of viscosity shear stress. One is the momentum exchange and viscosity dissipation stemming from random thermal motions of molecules, but its magnitude is relatively small. The oth- er is the continuous redistribution of molecules due to the velocity gradient, incurring a rotation of the net force field among contiguous molecules. k-,uC,,(or ,uCJ , D ~ P 1111. NAVIER-STOKES ( N S ) EQUATIONS 1. NS equation in vector form Introducing the constitutive equations for isotropic Newtonian Stokes fluid (i. e. , under the assumptions that ,u=const, 3I+2,u=0) , Eqs. (1. 2. 28) and (1. 2. 4 1 ) , into the equation of motion, Eq. (1. 2. 25 ) , we obtain the NS equation where V ( V V ) is a divergence of the gradient of vector V . In a rectangular coor- dinate system, the operator is called the Laplace operator (Laplacian) , also denoted by V *(or A ) , since it can be expanded by taking a formal operation like scalar prod- uct 14 (1. 2. 45) In this case, the components of V ( V V ) can be obtained by applying V to the correspondingcomponentsof V , i .e. , V (V V ) = V 2 V = V VV (SOV V is called a quasi-vector) . However, this relation does not hold in a more general coordi- nate system, when the operator should be expanded into (1. 2. 46) For an incompressible fluid, the last term on the right-hand side of Eq. (1. 2. 44) e- quals zero. If we have also an ideal fluid ( p = 0) , a reduced form of the NS equa- tions is called Euler equations. But even if p# 0, when all the derivatives &,/ax, are sufficiently small, the flow can still be approximated by the Euler equations. When DV/Dt= 0, the NS equations are reduced to a linear system, Stokes equations. The associated flow, Stokes flow (not to be confused with the Stokes fluid) , occurs when the inertial force is much smaller than the viscous force. We note in passing that historically only the equations of motion were called the NS equation. Now this term often means the complete system, including the equation of continuity (and additionally, the equation of energy for a general compressible fluid). There are important differences in structure and solution between the NS equa- tions for compressible and incompressible fluids. For incompressible fluids, only ve- locity appears in the equation of continuity, containing no density and pressure. This equation, as it is not in a complete form, is only partly coupled with the equation of motion and thus can be considered as a constraint on the velocity field. Pressure is in an unequal position to velocity, and cannot be determined by some thermodynamic condition (e. g. , equation of state). In numerical solutions, pressure and velocity are often obtained alternately in an iterative process. A difficulty lies in that the ve- locity obtained from the equation of motion under a given pressure generally cannot fulfill the equation of continuity. As for a compressible fluid, density appears in the equation of continuity, which is fully coupled with the equation of motion. If we take a fixed region in a flow, the pressure gradient at its boundary determines a velocity of the fluid leaving that re- gion, based on the equation of motion. That velocity incurs conversely a variation of density in the region based on the equation of continuity, and this has a feedback ef- fect on the pressure through the equation of state. In numerical solutions, pressure and velocity are often solved for simultaneously, and the 'elastic' constraint provided by the equation of continuity make the procedure easier due to inter-adjustability be- tween pressure and velocity. 2. NS equations in rectangular coordinate system for incompressible flow Later in this section we discuss incompressible flow only, but thereafter we shall study 2-D shallow-water flow from the viewpoint of compressible flow. The equation of continuity is (1. 2. 47) 15 Only the equation of motion in the x1 direction is listed below For the sake of brevity, it may be expressed by the use of a Cartesian tensor Eqs. (1. 2. 47) and (1. 2. 48) may be written concisely as (1. 2. 49) (1. 2. 5 0 ) (1. 2. 51) (1. 2. 52) We note in passing that space coordinates and velocity components will occasion- ally be denoted by (xl , x 2 , x3) and (u l , u 2 , us> , and sometimes by ( x , y , z ) and ( u , v , w). The latter are chiefly used in the 2-D case. 3. The NS equations expressed in terms of the stream function for 2-D incompre sible inviscid irrotational flows In the 2-D cases (plane or axis-symmetry) , it is possible to eliminate the pres- sure p from Eq. (1. 2. 52). Moreover, since for incompressible flows, p does not ap- pear in the equation of continuity, u and v can be combined into one unknown func- tion, so that only one equation is left. Defining the stream function @ by (1.2. 53) and introducing the above equation into the NS equations with p= 0 , we obtain (1. 2 .54) (1.2. 55) Taking derivatives of the above two equations with respect to y and x respective- ly , and then eliminating p , we have (1. 2. 56) When the flow is vortex-free, i. e. , the vorticity equals zero everywhere, the condition V X V = 0 can be expanded in the 2-D case, yielding 16 = o a u a v --- ag ax (1.2. 57) Introducing Eq. (1. 2. 57) into Eq. (1. 2. 53) , we get a Laplace equation satis- fied by the stream function v2*= 0 (1. 2. 58) Such an irrotational flow is called a potential flow. Substitute Eq. (1. 2. 58) into Eq. (1. 2. 561, and solve for the stream function under a given boundary condition. Then by taking a derivative in Eq. (1. 2. 53) we get a solenoidal velocity field also satisfying the NS equations. Lastly, substitute the velocity into the NS equations to obtain the pressure gradient. The whole problem has now been solved. An advantage of this approach is that the difficulty stemming from nonlinear convective terms can be avoided. In an extension to the 3-D case, a scalar velocity potential function can be in- troduced to replace the three velocity components. For a steady, irrotational and isentropic flow with pressure as the only external force, the NS equations can again be reduced to a single PDE in terms of the velocity potential. 4. Simplification of the equation of motion for incompressible flows For a general incompressible flow, the stream function cannot be used advan- tangeously, so it is necessary to adopt another approach in order to cancel out the e- quation of continuity and so eliminate pressure from the equation of motion. Firstly, we introduce a theorem : Any vector field w defined on a region can be uniquely decomposed into w=V+ V p , where vector field V satisfies: (i) div V= 0; (ii) V n= 0 , where n is an outward normal vector to the boundary of Sa. Physical- ly, it expresses the continuity requirement and the boundary condition that vector Ir is in its tangential direction. Based on the theorem, we define a linearorthogonal projection operator P , such that Pw = v , P ( V p ) = 0 (1. 2. 59) Then, we have v = P V , w = Pw + vp Applying P to the dimensionless NS equations (cf. Section 2. 31, we obtain (1. 2. 60) P (1. 2. 61) where Re=Reynolds number. In the above equation, p has been eliminated so that aV/& depends only on V . Let (1. 2. 62) w = - v . V V -f GAV Pw is just the time rate of change of velocity, which satisfies the above-mentioned two conditions. According to the definition of P , Pw=V, Eq. (1. 2. 61) is reduced to aV/at=V. Solve for 8 , introduce V into Eq. (1. 2. 62), and lastly, decompose w 1 17 as in Eq. (1. 2. 60) to get the pressure gradient, Vp=w-Pw=w-V. 1 . 3 TWO-DIMENSIONAL SYSTEM OF SHALLOW WATER EQUATIONS (2-D) SSWE I . NON-NEWTONIAN FLUID AND TURBULENCE In rheology a mathematical equation which expresses the relationship between shear stress and shear deformation rate is called the rheology equation. Its graphical representation has a name, a rheologic curve, which, in the case of Newtonian fluid, reduces to a straight line passing through the point of origin and invariable in time. All fluids that do not follow such a law are non-Newtonian. Among them, purely viscous fluids form the chief category, for which shear deformation rate at any point is a nonlinear function of shear stress at the same point only, and does not depend on other factors. Three simpler types are frequently encountered in this category, i. e. , Bingham fluid, pseudo-plastic fluid and expansive fluid. The rheologic curve for a Bingham fluid is a straight line (or a curve close to a straight line) which does not pass through the point of origin, that for a pseudo-plastic fluid is a upward-convex curve passing through the point of origin, while that for a expansive fluid is a down- ward-concave curve also passing through the point of origin. For a class of simple nonlinear viscous fluids, the condition that bias stress ten- sor z is proportional to velocity gradient tensor V V , is replaced by a general relation (1. 3. 1) Under isotropy condition, the stress tensor u has a general expression (1. 3. 2) For a compressible fluid, the coefficients a and /? are functions of temperature T and the two invariants of the symmetric deformation-rate tensor E. (The invariant of a tensor is defined as the quantity that does not change under orthogonal coordinate transformations. ) Here the two invariants are the determinant I E 1 and [(trEI2- trE2]/2. A fluid with Eq. (1. 3. 2) as the constitutive equation is called a Reiner- Rivlin fluid. Water flows with sediment concentration beyond some bound (such as wet clay and mud) can be treated as this kind of fluid. It has been found that the properties of some non-Newtonian fluids are more complicated than those admitted by Eq. ( 1. 3. 2). For example, a linear or nonlinear viscoelastic fluid has time-lagging , relaxation and wriggle characteristics, and these are the chief forms of energy dissipation. The stress depends not only on the present deformation, but also on the history. In other words, such fluids have their own ‘memory’. Besides, there is a class of viscoplastic fluids. When the stress exceeds a yield limit, the strain is proportional to the excess. A laminar flow in which molecular viscosity dominates can be treated as a New- tonian flow. A turbulent flow following some square-resistance law is similar to a non-Newtonian flow, but has a physically essential difference from it. Turbulence is a feature of the flow taking vortex size or mixing length as its length scale, while non-Newtonian property is a mathematical-physical description of the flow taking molecular mean free path as its length scale. u = u ( L , p ,TI u =- pl + aE + /?E2 18 I I . STATISTICAL TREATMENT OF TURBULENT FLBW- REYNOLDS EQUATIONS The fact that flow variables (e. g. , stress and velocity) change randomly in space and time is a basic feature of turbulent flow ; it can be treated as a sum of aver- age movement and irregular fluctuation. For example, velocity V at any instant may be decomposed into a time-averaged velocity 7 and a turbulent velocity I/' The mathematical expectation of Vt equals zero. We substitute Eq. (1. 3. 3) into the equation of motion, make a time-averaging, and simplify the resulting equations so that only time-averaged physical quantities appear. For incompressible flows, the e- quations can be expressed in Cartesian tensor form (Reynolds equations) as follows. V = B + V ' (1. 3. 3) 1. Equation of continuity - 0 ai, ax, _ - (1. 3. 4 ) which describes the mass conservation of the averaged flow. As for compressible flows, we have (1. 3. 4 a ) where the last term on the left-hand side represents mass transportation brought about by velocity pulsation around the averaged flow, often called dispersion. 2. Equation of motion Introduce a virtual stress tensor given by 6, 6, - Z,] = ,u - + - - p u',u; i 3x1 ax.) (1. 3. 5 ) (1. 3. 6 ) where z,, denotes the stress in the 2,-direction exerted on a plane perpendicular to the x,-axis. Hence (1. 3. 6 a ) It can be seen that the Reynolds equations are in the same form as Eqs. (1. 2. 23) and (1 . 2. 4 4 ) . Accordingly, the two terms on the right-hand side of Eq. (1. 3. 6 ) are called laminar viscosity stress z1 and turbulent viscosity stress zt , respectively. The latter is also termed Reynolds stress, which is a correction to laminar stress, and indeed is much larger than the former in almost the whole flow region, i. e. , except in a thin boundary layer. Strictly speaking, z1 is inherently not a stress, but comes from time-averaging the convective terms, leading to a better name of effective (or apparent) stress. It has a physical meaning of the influence of inertia which gives rise 19 to velocity pulsation , momentum transportation and energy dissipation. I I I . TURBULENT MODELS-THE CLDSURE PROBLEM OF REYNOLDS EQUATIONS The Reynolds equations are open since turbulent velocities appear. In order to describe a turbulent flow field, two types of information are required, i. e. , correla- tion both between random turbulent velocities in all coordinate directions at any point and between those at any two points. The former can be illustrated by an example. In a 2-D near-wall flow, denote velocity components tangential and normal to the wall by u and v , respectively. Normal Reynolds stress, -ur L u r , , is negative, whereas shear stress, -ur ,ur ~, is positive, so there exists a negative correlation between u and v , which determines various sizes of turbulent vortices and possible forms of energy transportation among them. However, the problem cannot be solved by a statistical approach. A most commonly-used approach is to establish a deterministic turbulence model so as to close the system of equations for the time-averaged flow and turbulent transportation. A turbulence model can be expressed as a system of algebraic or differential e- quations. According to the number of differential equations used, it can be classified into 0- , 1- , 2- and multi-equation models. There are two 0-equation models encoun- tered frequently : one introduces a turbulent viscosity analogous to molecular viscosi- t y , which is set to be constant everywhere in a flow field (or determined by using an empirical turbulent viscosity formula for inner- , overlay- and outer-layers respec- tively ) ; the other assumes that shear stress is proportional to the velocity gradient based on Prandtl's mixing length hypothesis. 1-equation models generally establish a transport equation in terms of turbulent kinetic energy k , which is related to turbulent shear stress z t , to describe the production, convection, diffusion and dissipation of k. The most famous %equation model is perhaps the k-E model, consisting of two trans- port equations in terms of k and E(dissipationrate of k ) , or two other similar equa- tions (e. g. , in terms of k and mixed length l ) . One of the most noticeable multi-e- quation models is a 3-equation model taking stress z, , k and E as unknowns, and an- other one is expressed in terms of turbulent shear stress and transport flux. Boussinesq' s approximation is the most classical technique for the establishment of an empirical relationship between the Reynolds stress and time-averaged velocity (a 0-equation model). By introducing a turbulent viscosity coefficient the Reynolds stress is written as (1. 3. 7) fiL and v1 are called dynamic and kinetic turbulent viscosity (or eddy viscosity, turbulent momentum exchange, turbulent diffusivity ) coefficients, cooresponding to /L and v respectively. As stated above, they do not originate from the viscosity prop- erty of the fluid but from the vortices produced by turbulence, and they depend on the characteristics of the flow field such as the velocity gradient. Turbulent and molecular viscosity are also quite different from each other in the scale and strength of fluid motion, as p<<pl, generally with a ratio about lo-*. Then the for- mula of turbulent stress Eq. (1. 3. 6 ) , (1. 3. 6a) become 20 (1. 3. 8 ) (1. 3. 8a) In a 2-D shallow-water flow computation, the simplest 0-equation model has been used in most cases, and even the Reynolds stress term is replaced by some em- pirical hydraulic resistance formula (cf. Section 1. 4 ) . When the turbulent structure varies gradually along a flow, a 2-D flow computation based on the Reynolds equa- tions and some turbulence model has been more or less successful. For example, at- tempts have been made to use a depth-averaged k-& model for calculating the distribu- tion of depth-averaged turbulent viscosity. It seems that at present the k-e model is the most promising one for calculating flow fields with small-scale circulation, espe- cially when there exists a separation between the fluid and the solid wall, or there are wakes behind an object. N . TURBULENT EDDY VISCOSI!lT In general, turbulent viscosity, v, , is an order-2 asymmetric tensor composed of nine components, whose values are subject to large variations with the flow behavior (near-wall turbulent flow or free turbulent flow-including jet, wake, etc. ) and ve- locity distribution, depending on the mechanism of turbulent energy production. In a flow field, the larger the shear deformation due to velocity gradient, the greater the value of v,. It is also affected by wind stress exerted on the water surface. Thus, vis- cosity should be a function of space and time, and it also depends on the space-time scale of flow (it increases with the scale of flow). In numerical solutions, however, because vortices are reproduced by using.a mesh, which imposes a limitation on their size, vt is often assumed to be a constant for the sake of simplification. In this case, the solution may be close to that for a viscous laminar flow on account of the similari- ty between the two expressions for z, so sometimes it is called a quasi-laminar flow. The value of v, may be chosen empirically according to the state of flow or through laboratory experiment, but an estimate based on observations can only be applied to dynamically similar flows. Of course, this choice would be influenced by personal judgement, leading to a range as large as several orders of magnitude and yielding quite different solutions unavoidably. Substituting Eq. ( 1 . 3. 8) with k= 0 into Eq. (1. 3. 5 ) , and assuming that all the components of turbulent viscosity are identical, we have (1. 3. 9 ) When turbulent viscosity varies in a flow field, the order-2 viscosity terms in Eq. ( 1. 3. 9 ) , should be written in the form of - ( v 2) . For brevity of notation, we shall omit below all the time-averaging marks over the symbols of various flow variables. a a ; axJ ax, 21 Seeing that in shallow-water flows, vertical velocity and acceleration are much smaller in magnitude than horizontal ones, we distinguish horizontal viscosity V~ from vertical viscosity v,, only, so that 1 a; =- - - + ipBi + v,&, v, = va( i = 1 , 2 ) or v,"(Z = 3) (1. 3. 10) at P ai In applications their values adopted are in a range wide enough, reaching 1 - 1 O3 and 1 OP4 - 1 0-'m2/s, respectively. The commonly-used range of values of vth is 5- 100 m2/s, while that of vtV in a middle layer of seawater is 10-2-10-4, and at the sea surface 5 X 1 O-'- 5 X lo-', corresponding to wind speed of 5-22 m/s. In a flow of geophysical scale, the orders of magnitude for vth and vl. are 1 O2 and 1 0-2, respec- tively. The cause of such a great difference can be interpreted as follows. Horizontal momentum exchange is affected chiefly by the vortices generated by the geometry of the boundary of a water body, while vertical vortices leading to vertical exchange have three sources, i. e. , underwater topography, wind over the water surface, and turbulence due to a vertical gradient of horizontal velocity. As shallow-water is char- acterized by its relatively large horizontal length scale, the energy and sizes of hori- zontal vortices are also much larger than vertical ones, but with much lower wave frequencies. The former influences mainly the distribution of the velocity, and has only a slight effect on water level, while frequency is the chief mechanism of dissipa- tion. Besides the simplest assumption of constancy of v w , there are several empirical formulas in use. ( 1 ) Neglect the variation of v , ~ with depth, and assume that v, is a bilinear func- tion of depth and depth-averaged velocity, e. g. (1. 3.11) where C=an empirical coefficient. ( 2 ) Assume that va is a linear function of the horizontal gradient of horizontal velocity. (3) Based on the Prandtl-Karman mixing-length theory, vw is estimated by vlh = Ch +/m (1. 3. 12) where the Karman constant K-0. 4 , and y is a height above bottom. Empirical formulas for vrV are exemplified in the following. ( 1) In the studies of European offshore waters in recent years, ignore the varia- tion of vth with depth, and assume that it is approximately proportional to the square of the depth-averaged velocity k V," = - (UZ + u2) U (1. 3. 13) where the dimensionless constant k= 2 X 1 O-5. Another constant u, the frequency of longwave (similar to the Coriolis coefficient of geostrophic motion), may be taken as u = lo-' 1/s. If it is required to consider the variation of vtu with depth, the above value is applicable to the middle layer of the seawater. At the sea surface, we need to 22 consider the effects of wind speed and water waves, whose height and period are re- lated to wind speed w,(m/s) and wind run L (km). So vt,, can be estimated by an empirical formula, e. g. , vt,, = (0. 1695 X 1 0-4 )Lo . wi, 6. At the sea bottom, v,,, de- pends on velocity and characteristic roughness length. From the above a vertical pro- file of vt,, can be established. (2 ) Let v,,, be a function of the square of the vertical gradient of local horizontal velocity, e. g. (1 . 3. 1 4 ) where I= mixing length showing the turbulent scale, determined empirically. ( 3 ) For a large-scale shallow-water flow, vt,, can simply be taken as a linear or parabolic function of depth. In a computation for the Aegean Sea, the following for- mula was used ( 1 . 3. 15) where u , =frictional velocity a t water surface, u* = m p ; t,” =shear stress a t water surface; 3, = calibration coefficient of the same order as 1 , h = 0 ( 1 ) ; z = height above mean sea level; d=height between sea bottom and mean sea level. If the fluid density varies in the vertical direction, the right-hand side of the above e- quation should be multiplied by a certain empirical function of the Richardson number (cf. Section 2. 2). A deeper study should utilize some turbulence model, to get the spatial distribu- tionand time-variation of vth and vtrl. The procedure, however, has its limitations and difficulties. This is because there are both large and small vortices appearing in a tur- bulent flow. Large vortices are nonisotropic, take energy from the main flow and continue to stretch in the course of motion until their breakdown into many small vor- tices. Small vortices, on the other hand, are isotropic, and dissipate energy when eventually reaching the smallest size. Therefore, it is difficult to simulate all sizes of vortices with only one model. V . DERNATIO,V OF 2 - 0 S A W E Two approaches can be made in the derivation of the 2-D SSWE. 1. Integrate the 3-D system of equations (except the equation of motion in the x3-direction) over water depth (or make a depth-averaging). In the integration the following boundary conditions are used. At the top of the water body ( i ) kinetic boundary condition (ii)dynamic boundary condition 7 = (a,,, At the bottom of the water body p = p. (1. 3. 16) ( 1 . 3. 17 ) U ] = 112 = U3 = 0 For example, the equation of continuity becomes 23 (1.3. 18) Exchange the order of differentiation with integration and intrpduce the boundary conditions, yielding (1. 3. 19) where h = water depth, h = z - z,, ; z = water level ; zb = bottom elevation ; N , V = depth-averaged velocity in the x-and y-directions , namely hI1 , hV =discharge per unit width in the x-and y-directions , also denoted by qz and qy. In oceanography they are called mass transportations or total currents, which can be taken as dependent variables to get a system of total-current equations (cf. Section 1 . 5). The 2-D equations of motion can be derived similarly. In the following we shall adopt another direct approach, which has a more explicit physical meaning. 2. Consider the 3-D flow as the following compressible plane flow. On the x-y plane, a virtual fluid with unit height in the z-direction flows a t the original depth- averaged velocity. It has a virtual density ph , such that the mass in a rectangular flu- id column with both sides of unit length, remains unchanged. Pressure p exerted on a lateral face of the column equals the total pressure over the whole water depth pgh2/2 (under the hydrostatic pressure distribution hypothesis). Now we proceed to derive the equations of motion for the column (Fig. 1. 1 ). We have : Fig. 1. 1 2-D shallow-water flow model Rate of change of momentum of the fluid column 24 Pressure difference between the two lateral faces at a distance dx Difference in shear stresses on the top and bottom faces Difference in shear stresses on the two vertical side faces at a distance dy z s z - z b z Body force in the r-direction : and po =atmospheric pressure at the water surface ; z., = r-component of wind stress at the water surface ; zbz = r-component of bottom friction force ; z,, , z,, = depth-av- eraged bias stresses; FR=body force per unit mass. phF8, From Newton's second law, we have (1. 3. 2 1 ) Then by analogy with the derivation of the NS equations, for the last term on the right-hand side we have which, under the assumptions that ,k=const and that h=O, reduces to (1. 3. 2 2 ) X J av & ay and under the additional assumption that - f - = 0 , reduces further to P A [ [ . An- other approximate derivation is based on the Eq. (1. 3. 8) with k= 0. Since by analo- gy with the Boussinesq approximation we have (1. 3 .22b) (1. 3. 2 2 c ) 25 The depth-averaging symbol over p, will be omitted hereafter, but keeping it in mind is necessary, because turbulent viscosity often varies significantly over a verti- cal. Moreover, we shall often use the symbols u and v instead of U and V . Substitut- ing Eq. (1. 3. 22) into Eq. (1. 3. 21) gives the desired 2-D SSWE, many of whose forms will be given in Section 1. 5. Now v, means the depth-averaged horizontal turbulent viscosity, which is differ- ent to the local viscosity in the 3-D case. According to a proposal made by the Delft Hydraulics Laboratory (DHL) , its value is about 1-2 m2/s, or it can be estimated by k2 v, = c, - e (1. 3. 23) where C,= an empirical coefficient; k=turbulent kinetic energy per unit mass, k = L(Lz f f G2) ; and &=dissipation rate of k. k and E can be calculated by using 2 a k-e model. Furthermore, in view of the difference in turbulent viscosity in the normal and tangential directions, make a correction to the formula of normal stress based on the turbulence theory (cf. Eq. (1. 3. 8 ) ) , so that a unified turbulent viscosity v, can be used on the horizontal plane, yielding 7,s a u 2 P ax 3 = 2v, - - -k - (1. 3. 24) (1. 3. 25) a2u a2v It can be seen that t, and T* come from depth-integration of p, - and pL - , a22 ,322 giving pi - and p, - respectively when viscosity is a constant. 3:: av az I ib Therefore, the physical meanings of and r,, are the vertical turbulent stresses at the top and bottom, which are expressed by the vertical gradient of horizontal ve- locity and vertical turbulent viscosity. When order-2 derivatives and nonhomogeneous terms vanish, the 2-D SSWE has the same form as the Euler equations for compressible ideal fluids. Since order-2 derivative terms are often neglected, and since the mathematical manipulation of non-homogeneous terms is relatively simple, we often use the latter model in theoreti- cal studies and numerical algorithms. It is noted in passing that for a 1-D gradually varied shallow-water open flow over a highly curved channel bed, the trajectories of fluid particles are curvilinear in the directions tangential to the surface and bed, so that vertical acceleration is not negligible. Horizontal velocity can no longer be assumed to be a uniform distribution over a vertical ; moreover, the pressure term includes, besides hydrostatic pressure, the effects due to the curvatures of streamlines. The 1-D shallow-water equations should be modified into, e. g. , the Boussinesq or Dressler equations containing, be- sides bed slope, additional terms involved with bed curvature and its space-variation. 26 I . 1 PHYSICAL MEANINGS OF VARIOUS TERMS IN 2-D SSWE I . LOCAL ACCELERATION Local inertial terms such as &/at , etc. , represent the time rate of change of ve- locity at any fixed position, and are the only terms showing nonstationarity of a flow. I I . CONVECTIVE ACCELERATION Convective acceleration terms such as u&/az, etc. , represent the effect of the space-gradient of velocity being transported together with a flow. Applying a curl op- erator to the equation of motion for deriving a differential equation of vorticity , also shows that these are the very terms governing the production and transportation of vorticity. They can be further categorized into two groups. The first group, non- cross convective terms such as u&/&, bear information on the velocity gradient in the same direction as the velocity. The second group, cross convective terms such as vaU/&, provides information on the velocity gradient in the other coordinate direc- tion. From the mathematical viewpoint, it is just these terms that make the system quasi-linear , so that a numerical solution could come to suffer from nonlinear insta- bility (cf. Section 10. 3 ) . The larger the Reynolds number, the more important role they play in the system. This fact explains the difficulties occurring in computations with high-Reynolds-number flows. For example, the coefficient matrix of the linear system of equations obtained by discretization may have no diagonal dominance (cf. Chapters 5 - 7 ) , so that the iterative process would show instability and thus noncon- vergence. The sum of the above two terms is just a material derivative, e. g. h/&, which represents the total acceleration of fluid particles, called the inertial term. By neglecting the convective acceleration, the system becomes linear. This ap-proximation suits the case of low Reynolds numbers (low velocity or high viscosity). Experiments show that such an approximation is favorable to computational stability, and that satisfying a linear stability criterion would be sufficient. However, a mech- anism of vorticity production and transport is then lost, so that various types of circu- lations and vortices can no longer be simulated. Surface slope terms such as gaz/&, etc. , represent the action of gravity. For a water flow with a free surface, these terms are often the chief driven forces, so the associated waves in an open flow are called gravity waves. As stated above, they o- riginate from the assumption of hydrostatic pressure. In theoretical studies, we often 27 ah azb ax ax prefer to decompose it into pressure gradient and bottom slope, g(- + -) . The first part shows the pressure gradient due to variation of depth, while the second part shows the effect of underwater topography, which acts as an external force. In this form, h (not z ) is an unknown in all equations. But in practical computations, they are often combined into one term, &/ax, in order to minimize discretization errors, as the surface slope is usually much smaller than the bottom slope. It should be emphasized that it is just this term that results in one of the most sig- nificant differences between the SSWE and the Euler equations for compressible flows. The Euler equations are homogeneous, while in the SSWE, even though all other forces except gravity are not present, the bottom-slope term in general would appear in the momentum equation. Of course, the structure and properties of the so- lution to the SSWE in homogeneous form remain the same as for the Euler equations, however, nonhomogeneous terms may have an important effect on the quantization of the solution. A system containing all the above three terms, is called a dynamic wave model. Neglecting convective terms leads to a diffusive wave model. Neglecting pressure gra- dient leads further to a kinetic wave model. If the bottom slope and bottom friction are neglected, a gravity-wave model results. IV. ATMOSPIIERIC PREiWlJRE GRADIENT 1 ape Terms such as - - , etc. , represent the action of the atmospheric pressure field. They can be ignored in many cases, but in storm-surge forecasting they may be one of the chief factors that must be considered. An atmospheric pressure field can be given based on observation or on forecast. In the absence of data, we may use a field associated with an ideal typhoon (hurri- cane) or cyclone model exemplified below. P a x 1. Fujita formula (1952) The atmospheric pressure at a distance T from the center of a typhoon can be es- timated by (1. 4. I ) where p , =environmental atmospheric pressure not affected by the typhoon ; p o = pressure at the center of the typhoon ; and R=radius associated with maximum wind speed, i. e. , the extent of typhoon, the value of which varies in the course of devel- opment of the typhoon, and also varies during the day (being larger in the day than at night). It can be taken directly from observed data of the atmospheric pressure field, or estimated iteratively with the least-square method. Similar formulas suit symmetric atmospheric pressure fields. p, (r> = p- - ( p - - Po>/ 41 + ( T / R > * 2. Takahashi formula p G ( T ) = p m - ( p - - p o ) / ( l + T / R ) (1. 4.2) 3. Meyer formula 28 4 . Telesnianski formula (1. 4. 3) (1. 4. 4 ) The variation of atmospheric pressure causes a rise or fall of water levels. In the equilibrium state, a water level has 1 cm of difference for every 1 mb variation of atmospheric pressure. (mb is abbr. of milli-bar. In the SI unit system, 1 bar= l o 5 Pa(N/mZ)=O. 986923 atm= 10. 1972 mH,O. ) When a water level at the bound- ary of an open sea cannot be given, it may be estimated according to the difference between the actual pressure pa and standard pressure p a = 1 atm= 1012 mb. Generally speaking, the larger the water depth, the stronger the influence of the atmospheric pressure field. An analysis made by Timmerman leads to the conclusion : When wind speed w, is higher than 15 m/s and water depth h is greater than 200 m , the effect of the atmospheric pressure field exceeds that of wind stress exerted on the water surface. In contrast, a case study shows that, when h<50 m and the pressure variation is small, its influence can be neglected. In addition , in the deep sea, an at- mospheric pressure field may create a flow extending down to the sea bottom, a fea- ture which is different to that caused by surface wind stress. V . SURFACE WIND STRESS Surface wind stress terms such as zo,/ph, etc. , represent the drag force pro- duced by wind over the water surface. Turbulence appears at the air-water interface due to instability, so that water moves up and down forming a rough surface, which, similarly to sand waves on a river bed, affects the value of the drag force. Its action is inversely proportional to water depth, so it is important for shallow-water flow. z. can be estimated by a semi-theoretical formula, based on a similarity hypothe- sis proposed by Karman for turbulent flow fields. At a point at a given distance z from a rough interface, turbulent properties are independent of molecular viscosity, and similar turbulent structures exist in all cases only differing on the space-time scales, and depending on the order- 1 and order-2 derivatives of the average flow with respect to z. From this a logarithmic law for average turbulent velocity has been de- rived, showing that the turbulent friction force at the interface z,, is proportional to the squared turbulent velocity u (now wind speed) (1. 4. 5) where ?~=Karman constant: p,=air density, at 0 ‘C and 1 atm p,=l . 293X lop3 g/cm3; k=roughness length of water surface. In oceanography, when the wind force is medium or strong, it is usual to take k = 0. 6 cm regardless of wind speed. Substi- tute z= 10 or 15 m into Eq. (1. 4 . 5) , giving z a = POC, I zo. 120. (1. 4. 6) 29 where w.=wind speed at height z with unit m/s, ( lm/s=3. 6 km/h, 1 knot= 1. 852km/h). According to the extended Beaufort wind force table, the wind speed at a height of 10m is related to the wind scale as follows: wind scale 0 1 2 3 4 5 windspeed 0-0. 2 0. 3-1.5 1.6-3. 3 3. 4-5. 4 5. 5-7.9 8. 0-13.8 wind scale 6 7 8 9 10 11 wind speed 10. 8-13. 8 13.9-17. 1 17. 2-20. 7 20. 8-24. 4 24. 5-28. 4 28. 5-32. 6 wind scale 12 13 14 15 16 17 wind speed 32.7-36.9 37. 0-41. 4 41.5-46. 1 46.2-50.9 51. 0-56. 0 56. 1-61. 2 C,=a dimensionless drag coefficient (i. e. , resistance coefficient at height z ) . Under the condition that the atmosphere is in a neutrally state of stability, C,, can be calcu- lated by using one of the following empirical formulas: 1. 1 X l o - 3 ; when w,=20 m/s, Cn=2. 6 X polation is used C, = (0. 9 + 0. 0 8 2 ~ ~ ) X Cn = (0. 63 f 0. 0 6 6 ~ ~ ) X l o p 3 1. Wilson formula (1960) Take z=10 m. When w,=2. 8 m/s , Ca= In the above range, linear inter- (1. 4 . 7) (1. 4. 8 ) (1. 4. 9) 2. 10s (Institute of Oceanographic Sciences, UK) formula 3. Garrat formula (1971) 4. Heaps formula (1965) Take z= 15 m. c,, = (0 .75 + 0 . 0 6 7 ~ ~ ) x 10-3 cD = 0.565 x 10-3 (w. < 5 CD =- 0. 12 + 0. 13720. C,, = 2. 513 X lop3 (w. > 19. 22 m/s> (wo = 5 - 19.22 m / s ) (1. 4. 10) In a storm-surge-tide computation for the North Sea, take z=10 m. When wa<15 m/s, Cp= 1. 7 X l o T 3 ; when w.>20 m/s, Cp= 2. 5 X 1 0-3 ; when w, is in the above range, linear interpolation is used. Wind field data may be taken from observation or forecasts. In addition, we may calculate geostrophic winds based on the atmospheric pressure field, and then estimate the sur- face wind speed by using some empirical relation between geostrophic wind and sur- face wind. To estimate the wind direction, it is necessaryto take into consideration a bias angle between geostrophic wind and real wind ( e. g. , about 18" counter-clock- wise in the northern hemisphere). As regards typhoons, in the literature it is some- times assumed that the wind velocity is directed toward its center and makes an angle with the isobar lines, which is taken as 30" for moderate-latitude zone in the nothern hemisphere. The wind field thus obtained is called gradient wind, which is added to the velocity of the moving center. In a numerical simulation of storm-surge tides for the North Sea made in the Netherlands, first of all, the actual atmospheric pressure field on the water surface every 3 hours and a forecast thereof every 6 hours are giv- en. Each pressure field is approximated by a n incomplete Fourier series, i. e. , one consisting of a finite number of terms. At any other instant, a field can be obtained by linear or quadratic interpolation on these Fourier coefficients. Calculation of the 5. Rijkswaterstaat formula 30 wind velocity at all grid points based on the pressure field data can be carried out by using - - - _ dw2 - d t p. ax ap + lw, - cw, (1. 4. 1 1 ) (1 .4 . 12) where I= f +bsin fi and c= bcos 6; b=ratio between friction force and wind speed; fi =angle made by negative wind velocity vector ( - w.) with friction force, measured counter-clockwise; f=Coriolis coefficient. c and 1 depend on the stability degree of the atmosphere, which is in turn related to the difference of air temperature and wa- ter surface temperature To-T,. To get 1 and c , read off T,,-T,, from synoptic charts and look up meteorological diagrams and tables. The above two formulas can be reduced to Let then in Eqs. (1. 4. 13) aG aG A , = AG, f BG, - C -f - D ‘I at at aG, aG A2 =- BG, f AG, t D - - C - at at aG aG A3 = 1 f C 2 f D 2 ax ax aG, aG A d = l - D - + C 2 ay aG X, a9 ay A5 = C L f D- aG aG A6 = - D _f f C -2 ax ax ( 1 . 4. 13) (1. 4. 1 4 ) We may also use the approximate formulas 31 0 . 7 a p 0 . 7 ap = f ax' f ay w = - - w = - - (1. 4. 1 5 ) where 0. 7 is the approximate value of friction coefficient. model. The following is one of such (in the CGS system) In the absence of data, we may use some wind field formula as an ideal typhoon w,, = c,l.',exp 7 - cz(rs inp + ycosp) [ 5 x T T 0 ] where c1 and c2 =empirical coefficients in the range of ( 4 / 7 - 6 / 7 ) and (0. 6 - 0. 8) , respectively, depending on the features of the typhoon ; V,, Vy =x- , y-compo- nents of the velocity of the typhoon center located at the origin ; cp= angle made by gradient with isopiestic lines; f=Coriolis coefficient. Moreover, ( 1 . 4. 1 7 ) In India, an idealized cyclonic wind field was described by the following formula (1. 1. 18) where V is a horizontal wind speed. In a simulation of the cyclone attacking the Andhra Pradesh Coast in Nov. 1977, took Vo=70m/s, R=80 km, c=240km. I ' I . BOTTOM FRICTION Bottom friction terms such as rbz/h, etc. , have a nonlinear effect of retarding the flow. When the strength of turbulence becomes stronger, the effect of molecular viscosity becomes relatively smaller, while viscous boundary layer near a solid wall becomes thinner, and may even appear not to exist. At any instant establish the 2-D Reynolds equation with p= 0 in the direction normal to the wall and make a time-in- tegration , yielding Zh = - p 21'20' ( 1 . 4 . 19) where u' =the pulsatile velocity along the wall, and w' = the pulsatile velocity per- pendicular to the wall (directed to the water body). This equation shows that the bot- tom friction is equal to the bottom turbulent stress. In 1942, Prandtl made the assumption that bottom friction resistance is a known constant ; in other words, there is a layer with constant stress. He also assumed that turbulent viscosity, v,, varies linearly with depth; i. e. , vL=xz, where ~ = K a r m a n - h constant and z is a height above bottom. Integrating the equation v, - = rh 9 we az know that the velocity profile over a vertical follows a logarithmic law. If it is possi- 32 ble to measure velocities at two heights zl and z 2 , we can estimate zb by the formula (1. 4. 20) However, we usually estimate q by using an empirical or semi-empirical formula, since the bottom turbulent stress is not well understood, in addition, the vertical dis- tribution of horizontal velocity cannot readily be obtained. 1. Hydraulics approach. Recall that in the 1-D system of Saint-Venant equations the term s b / p h can be expressed as qS,, where S, denotes the slope of hydraulic friction. Assume that the frictional force in a 2-D unsteady open flow can be estimated by referring to the for- mulas for I-D uniform flows in open channels, e. g. , by the Chezy formula (1. 4 .21) where C=the Chezy coefficient ( G / s ) ; R =the hydraulic radius ( m ) ; and u= cross-sectional average velocity. C / 6 expresses the ratio of u and the frictional velocity u * , sou" = a . When there is no wind stress, it is also possible to es- tablish a relation between C and the vertical turbulent viscosity (1. 4. 22) C is determined by experience and its value for natural rivers is about 20-70 (from floodplain to deep channel). Perhaps, the Manning formula below has been the most frequently used of the square-of-velocity resistance laws of turbulent flow ; it has the form T,,, = 0. 07 . / g u h / C (1. 4 . 23) where n is the Manning roughness coefficient often viewed as a constant. However, for alluvial rivers, especially those with a sand bed, the situation is more complicated in cases where the frictional resistance is closely related to flow behavior near the bed, e. g. , the evolution of sand waves may double the roughness. An approach is to establish a functional relation for roughness, e. g. , the Qian-Mai composite resistance formula, n = K,; / ' /A , where K , is the size of roughness, which may be set to D65 taken from the graduation curve of bed silt. A is related to those flow variables that govern the development of sand waves, and approaches a constant value with dimin- ishing sand waves; it varies for different rivers. Conversely, when roughness is fixed, the power in the formula has to be changed. Based on studies of channels with moving beds in India, Lacey suggested replacing the Manning formula by u = 16R2J3Sj/3 , while Malhotra proposed another version, u = 18R5/3S:/3 . Both formu- las have their own conditions of applicability, and a general form can be written as u = CR'S!. Substitute Eq. (1. 4. 23) into Eq. (1. 4. 21 ) , yielding S, = n2u I u I (1. 4. 24) Here we write u2 as u I u 1 so as to express the direction of S, , so that it can be used for 33 to-and-fro flows. It is easily seen that the formula can be generalized to the 2-D case only approximately. In 1-D flows we do not differentiate between bottom and lateral friction, while in 2-D flows we often take a unit-width channel ( R = h ) and only consider bottom friction. Moreover, there may exist some interactions between flows in the x- and y-directions. Keep these in mind and note that the projections of S, in the x- and y-directions S, and S , satisfy s,= (1. 4 . 25) then we obtain Zb, = pghS,, = pgU d m / c 2 = P U d w = RbU (1. 4. 26) zbY = pghS,, = pgi' A / w / C 2 = pv A/- = Rbv (1. 4. 27) where the dimensionless bottom-friction coefficient T = g/C2. The bottom-friction co- efficient Rb = pr d-. When we use the Manning formula, Rb = gn2 d V / h 1 I 3 (1. 4. 28) There is an empirical relation between Rb and bottom vertical turbulent viscosity (1. 4. 29) The order of magnitude of Rb/p is 10-4-10-3m/s. Based on observations and investigations into offshore waters and estuaries, the mean value of T is in the range (2. 5-3. 3 ) X l o p 3 , corresponding to C=62. 6-54. 5. Indeed, T is not a constant. For example, it may
Compartilhar