- Email: [email protected]
A Three-Dimensional Finite Element Model for the Study of Steady and Non-Steady Natural Flows
359
A THREE-DIMENSIONAL FINITE ELEMENT MODEL FOR THE STUDY OF STEADY AND NON-STEADY NATURAL FLOWS J.-L. ROBERT and Y. OUELLET Civil Engineering department, Research Center on Computer application in Engineering, Lava1 University, Sainte-Foy, Quebec, CANADA, G1K 7P4 ABSTRACT This paper will present the basic concepts of a tridimensional finite element model for the simulation of natural flows. In a first step, the fundamentals of the formulation will be described with an emphasis on initial assumptions, mobile layer concept, type of element, discrete variables of the system. Then, special attention will be taken to expose the original technique developed to include, on the entire depth of water, the effect of the bottom friction and some justifying examples will be shown. Afterwards, a brief discussion on the strategy of resolution for such non linear system will point the specific importance of friction and inertial convective terms. This will be followed by the presentation of some examples showing the good general behaviour of the model, particularly its capability to well simulate the vertical circulation of a flow. 1 INTRODUCTION
After working since many years with finite difference and more recently with finite element depth average two dimensional models, some specific environmental problems in coastal and estuarine zones, involving the effect of density and wind on the hydrodynamic circulation, raise the need of a three dimensional model. The finite element approach was retained for mainly two reasons. The first one was its ability to discretize with variable density of computational points an irregular domain, as it can be found in the nature. This is an important advantage in a way that it allows us to choose the position of artificial boundary conditions far enough to minimize their influence and without using a prohibitive number of computational points. The second reason was linked with the possible future development of the model, specially the possibility of its coupling with existing 2-D F.E. hydrodynamic model. The major goal we kept during the development of the model was its capability to simulate the largest possible range of free surface natural flows, large or small scale, steady or non-steady, tide or wind driven. The domain of application itself has a direct influence on the formulation, in fact natural flows have a small verticalhorizontal aspect ratio. Then, the discretization technique, as it will be shown later, has been chosen differently depending on the vertical and the horizontal directions.
360 2 FORMULATION OF THE 3-D MODEL 2 . 1 Governing Equations
The physical behaviour of a fluid may be expressed by a set of equations based on the forces equilibrium and the principle of mass conservation. In a three dimension (x, y, z) domain the momentum equations are:
where: u, v, w : the velocity components; p : the pressure; p : the density; X, Y, Z : the external forces (Coriolis, gravity) T~~... :
the Reynolds stresses
The Reynolds stresses express the viscous similarity of the effect of the turbulence. In the model, they are included using the Boussinesq assumption using an eddy viscosity coefficient A:
T
x x = 2A 11
T
xixj
= A
(4)
(5)
Note must be taken that in the vertical direction, those terms will receive special treatments. The mass conservation in an incompressible fluid gives the continuity equation:
361 2.2 Vertical Integration
Since the model is specially designed to simulate natural domains, some characteristics of them can be observed and then stated:
- small vertical dimension in respect of horizontal ones; -
vertical acceleration of the flow negligible; small vertical velocity component.
Those statements eliminate from the range of applications high slope super critical flows, closed conduits and turbomachinery flows. With this restriction the assumption of hydrostatic pressure can be taken and the acceleration and stresses term of the third momentum equation vanish keeping only the pressure and the gravity terms. In this form the set of governing equations (1). ( 2 ) , (3) and (4) is integrated in the vertical z direction from the bottom position -h to the water level n. But, in order to get some information about the vertical structure of the flow, the integration process is split into several terms:
by the introduction of intermediate integration limits (fig. 1) which have their position located by: Lk =
n-(e) k
k = 1, 2 ,
... b
The presence of 11 in the equation (8) means that the position of those intermediate level must move proportionnally with the water level, consequently the relative thickness of a layer is constant. An uniform distribution of those level has been chosen here; however, it is possible to change those relative position by any other relationship including the water level to get another distribution. It is important to note that the intermediate integration surfaces do not have any physical meaning and allow the mass exchange between them. This is specifically prescribed by the intermediate boundary condition (Robert 1983):
n = k+l
a
Y
362
Fig. 1 Definition of the vertical integration boundaries. where Un and Vn are mean velocity components in a layer n. At the extremity of the vertical domain boundary conditions must be defined using a streamline condition: atl
atl
at
w(-h)
=
-u(-h)
ah - v(-h) ax
ah
The introduction of those conditions in the int gration pro ess using th Leibniz rule, does not eliminate all boundary velocity terms at the intermediate level, the assumption that those terms are included in the internal stresses term is made (Robert 1983). Then the model is defined as:
- Continuity equation
363
- Movement equations +
+ U k aauk x+Vkay auk -fV
e [>+
+
( TxLk-1 -
a
T
xLk
U kavk T + Vk avk + fUk
k = 1
))]
+
g
=
ahkTxyk aY
0
-- -
Y
The vertical velocity component is computed from equation ( 9 ) . The formulation of the horizontal stresses will be discussed in section 3. 2 . 3 Finite Element Formulation
This short overview of the finite element formulation will point out the specific choices made to obtain the discrete form of the model. The integral form of equations ( 1 2 ) , (13) and (14) is obtained by a Galerkin weighted residual method.
As one can see from the model equations, two dimensional elements can be used by associating at each node the water level
rl
and the two horizontal
components of the mean velocity in each layer Uk and Vk. Two types of interpolation functions were then considered. Based on the work of Cochet 1979 on 2-D hydrodynamic model, a quadratic function for the velocities and a linear function for the water level were tested. The result shown that in a multilayer configuration some instabilities appeared fig. 2. Then a full quadratic function was used and improved considerably the performance of the model. Two full quadratic elements have been developed, a six nodes triangular and a nine nodes quadrilateral elements. Special attention must be kept in the discretization of movement equations (13) and (14). Since the interpolation functions applie only on discrete parameter associated with one layer and are zero for any other parameters, the assembly process is equivalent to the sommation present in those equations. Then the resulting set of algebraic equation will have as much equations as discrete variables. The resulting non linear system in non steady form can be expressed by:
[Ml{il + CK ($)I {$I
=
0
The solution method will be presented In section 4.
364
.
-
I
I
Quadratic element Quadratic-linear element
5
Fig. 2
Computed velocities in the center of a schematic estuary presented at Fig. 5.
3 HORIZONTAL STRESSES 3.1 Formulation
One of the most important factor in the general behaviour of a natural mass of water is the horizontal shear stress. The predominant one is the effect of the bottom roughness over the entire depth of water and in some cases, the effect of the wind at the free surface may be determinant. In such case, this specific stress component may be expressed by:
365 Tx L = CtPa wx IWI 1 0 i where Ct is a drag coefficient, pa the air density and W the component of wind xi velocity W in the xi direction. For the bottom friction term, ChLzy formulation has been chosen:
TxiLb
3 vx i Ivl
=
Here, in the case of a two dimensional model the Ch6zy coefficient C , takes in account the overall effect of the bottom roughness on the flow over all the thickness of the water layer. However, when this depth is discretized, this term will affect only the bottom layer and the coefficient C must be corrected in way such as the instantaneous discharge must stay the same whatever the number of layer. To ensure the transfer of the bottom roughness to the upper layer, an internal shear stress component must appear at the interface of each layer. This may be expressed by discrete form of the Prandtl internal shear stress formulation; in the x direction we have at the upper interface of a layer k : =
TXLk-1
hk
Auk 4 AU; + AVk2
where Auk and AV are the differences of the flow velocity between the layer k k and k-1 and II is the mixing length. Many authors (Wang 1977, Kokayashi 1980, Leenderstee 1975) use a similar formulation, but the II coefficient is considered as an experimental adjustable coefficient. 3.2 Mixing Length Parameter
As we do not really know what is the value of the mixing length, expect near the bottom where Prandtl assume that it may be proportional to the distance from the bottom, we must find an usable formulation of the mixing length up to the free surface. Montgomery 1943 suggests, from some experiments in particular flow that:
x
is the Von-Karman coefficient which is defined by x = 0.4. As shown on fig.
3 , this distribution presents a plane symetry on the z/h = 0,5 plane. This
366
involves the assumption that the distribution of the horizontal stress is similar near the free surface and near the bottom, which is not really satisfactory. To remedy such situation, Yalin 1977 suggests a distribution of the mixing length based on the assumption of a logarithmic velocity profile:
A s the Montgomery 1943 distribution, this one agrees with the Prandtl assump-
tion near the bottom:
a
=
x z
(21)
Fig. 3 Mixing length vertical distribution. The major advantages of the Yalin 1977 distribution are: a) the mixing length decreases less quickly near the free surface than near the bottom, assuming that the turbulent fluctuation in the neighbourhood of the surface are less damped than by the presence of a rigid boundary at the bottom. b) the horizontal shear stress becomes zero without the necessity of assuming a zero vertical velocity gradient at the surface. This formulation was retained and adapted to the present model; in a discrete form, for the lower interface of a layer k, equation ( 2 0 ) becomes: k'
=
x
hk (b
-
k)
(22)
367
and the corrected bottom friction factor can be written:
where Cr is the Ch6zy reference coefficient for a two dimensional model. In order to verify the good overall behaviour of this approach, a test was run on the model to simulate a steady state uniform turbulent flow with a five level configuration. The velocity profile obtained is presented at the figure 4 .
t
L
y'm
K
..I
-*-
Fig. 4
.
-
-.u, m/s
Velocity profile for a uniform steady flow.
4 SOLUTION METHOD
In steady state configuration, the model of the flow is dominated by a quadratic friction term. The Newton Raphson procedure must be chosen rather than a substitution scheme to solve the non linear system. In fact, it is obvious that the later will not converge because the quadratic term. The incremental form of the method is:
368
where the tangent matrix is:
The non-steady flow is based on an implicit Euler scheme with an iterative Newton-Raphson stabilization procedure:
These techniques need much CPU time, and some possibilities have been set to reduce it. The tangent matrix can be, with a little loss of accuracy, recomputed at the beginning of each time steps rather than each iteration. In some case, computing time reduction may reach a factor three. 5 APPLICATION Many numerical tests have been performed with this model. The simulation of
the tridimensional structure of the flow in a 180'
curved flume in steady
state has been compared with a scale model in laboratory. The dimensions used in the mathematical model were the same than for the physical model to prevent any scale effect on the results. The test bring us the confirmation of the rightness of the formulation chosen here. The results of this comparison can be found in Robert 1983.
5.1 Schematic estuary In the testing process of the model, one goal was to verify that the vertical behaviour of the simulation agrees with the results of a two dimensional vertical models (width averaged) particularly with the Sunderman 1976 formulation based on cubic interpolation function in the vertical direction. The test case can be described as a schematic estuary closed at one end and submitted at the other one to a sinusoidal tidal cycle of 1 meter of amplitude. The flume has a length of 92 km and an width of 4.6 km. The depth is 50 meters except for a central portion where the bottom raise up to a depth of 20 meters. The finite element grid is composed of 10 quadrilateral elements and three layers were specified. The results of some time steps of a complete cycle (fig. 5) show a good agreement with those presented by Sundermann 1976 ones despite the fact that the vertical discretisation of the present model is far less accurate. The time interval of computation was. in this case, 1200 sec. for a full implicit scheme. This case was an occasion to test the validity of the mixing length distribution. The same test was run with a constant interfacial friction factor
369
Fig. 5 Vertical flow structure in a schematic estuary.
and shown some irregularities of the vertical velocity profile as an erroneous position of the maximum velocity during the ebb flow. That situation was cured by the use of the distribution presented by eq. 20. 5 . 2 Wind induced circulation
The second example focuses on wind driven circulation. Some preliminary numerical tests agreed with the Leenderstee 1975 samples in the same condition. As for preceeding example the simulation was performed as a propagation problem; a constant 40 km/h s-W wind blowing over a 10 km x 5 km x 1Gm rectangular bassin. The result, shown in fig. 6 , gives the circulation over the five layers, the water surface elevation and a transversal cross section in the middle of the bassin after six minutes of real time. At this time, the steady state is not totally achieved and the free surface is still moving. The corresponding computing time for this non-linear non-steady 2 2 2 1 degrees of freedom was about 2h40 on a VAX 850 for a full implicit scheme. Nevertheless, some techniques as the computation of the global matrix only at the beginning of each time step can reduce the computing time to one third and the use of a vectorial computer may be usefull. 6 . CONCLUSION
One of the most important concluding remark of this work, is that the automatic adjustment of the interlayer friction factor by the use of a well chosen mixing length distribution gives excellent vertical behaviour of the model. This vertical pattern of the flow is then relatively accurate and the formulation is well adapted to the natural flows. The use of the model is made easy by the fact that the finite element grid is only two dimensions and consequently a coupling with a two dimensional model is provided. This two and three dimensions models should be used to predict the vertical structure of the flow only where it is necessary.
371
-0,6 mrn
d
-
\
.
,
L
I
r
-
Fig. 6
r
c
Y
c
e
c
-
L-
c
+ Wind driven c i r c u l a t i o n i n a rectangular bassin.
0,51
312
7. REFERENCES Cochet. J.F., 1979. Modglisation d'gcoulements stationnaires et non stationnaires par la mlthode des 6lLments finis. ThPse de Docteur-Inghieur, Universitg de Technologie de CompiPgne, France, 142 p. Kobayashi, M., Nakota, K. and Kawahara, M., 1980. A Three Dimensional MultiLeveled Finite Element Model for Density Current Analysis. In: Third Int. Conf. on Finite Element in Flows Problems, Vol. 11. Banff, Canada, pp. 80-92. Leenderstee, J.J., Liu, S.K., Alexander, R.C. and Nelson, A.B., 1975. Threedimension Model for estuaries and Coastal Seas, Vol 11, Aspects of Computation, R-1764-OWRT, The Rand Corporation, California, USA, 165 p. Montgomery, R.B., 1943. Generalization for Cylinder of Prandtl's Linear Assumption for Mixing Length. Annuals of the New York Academy of Sciences, XLIV: 89-103. Robert, J.L., 1983. ModElisation tridimensionnelle des Ccoulements a surface libre, permanents et non permanents, par la mgthode des Blgments finis. ThSse de doctorat, Universitg Laval, Qudbec, Canada, 233 p. Sunderman, J., 1976. Computation of Barotropic Tides by the Finite Element Method. In: Finite Element in Water Resources, First Int. Conf. Gray, Pinder and Brebbia ed., Princeton University, Pentech Press, London, U.K.. pp. 4.31-4.67. Wang, H.P., 1977. Multi Leveled Finite Element Hydrodynamic Model of Block Island Sound. In: Finite Element in Water Resources, First Int. Conf., Gray Pinder and Brebbia ed., Princeton University, Pentech Press, London, U.K. Yalin, M.S., 1977. Mechanics of Sediment Transport, Second ed., Pergamon Press, Oxford, U.K., 259 p.
A THREE-DIMENSIONAL FINITE ELEMENT MODEL FOR THE STUDY OF STEADY AND NON-STEADY NATURAL FLOWS J.-L. ROBERT and Y. OUELLET Civil Engineering department, Research Center on Computer application in Engineering, Lava1 University, Sainte-Foy, Quebec, CANADA, G1K 7P4 ABSTRACT This paper will present the basic concepts of a tridimensional finite element model for the simulation of natural flows. In a first step, the fundamentals of the formulation will be described with an emphasis on initial assumptions, mobile layer concept, type of element, discrete variables of the system. Then, special attention will be taken to expose the original technique developed to include, on the entire depth of water, the effect of the bottom friction and some justifying examples will be shown. Afterwards, a brief discussion on the strategy of resolution for such non linear system will point the specific importance of friction and inertial convective terms. This will be followed by the presentation of some examples showing the good general behaviour of the model, particularly its capability to well simulate the vertical circulation of a flow. 1 INTRODUCTION
After working since many years with finite difference and more recently with finite element depth average two dimensional models, some specific environmental problems in coastal and estuarine zones, involving the effect of density and wind on the hydrodynamic circulation, raise the need of a three dimensional model. The finite element approach was retained for mainly two reasons. The first one was its ability to discretize with variable density of computational points an irregular domain, as it can be found in the nature. This is an important advantage in a way that it allows us to choose the position of artificial boundary conditions far enough to minimize their influence and without using a prohibitive number of computational points. The second reason was linked with the possible future development of the model, specially the possibility of its coupling with existing 2-D F.E. hydrodynamic model. The major goal we kept during the development of the model was its capability to simulate the largest possible range of free surface natural flows, large or small scale, steady or non-steady, tide or wind driven. The domain of application itself has a direct influence on the formulation, in fact natural flows have a small verticalhorizontal aspect ratio. Then, the discretization technique, as it will be shown later, has been chosen differently depending on the vertical and the horizontal directions.
360 2 FORMULATION OF THE 3-D MODEL 2 . 1 Governing Equations
The physical behaviour of a fluid may be expressed by a set of equations based on the forces equilibrium and the principle of mass conservation. In a three dimension (x, y, z) domain the momentum equations are:
where: u, v, w : the velocity components; p : the pressure; p : the density; X, Y, Z : the external forces (Coriolis, gravity) T~~... :
the Reynolds stresses
The Reynolds stresses express the viscous similarity of the effect of the turbulence. In the model, they are included using the Boussinesq assumption using an eddy viscosity coefficient A:
T
x x = 2A 11
T
xixj
= A
(4)
(5)
Note must be taken that in the vertical direction, those terms will receive special treatments. The mass conservation in an incompressible fluid gives the continuity equation:
361 2.2 Vertical Integration
Since the model is specially designed to simulate natural domains, some characteristics of them can be observed and then stated:
- small vertical dimension in respect of horizontal ones; -
vertical acceleration of the flow negligible; small vertical velocity component.
Those statements eliminate from the range of applications high slope super critical flows, closed conduits and turbomachinery flows. With this restriction the assumption of hydrostatic pressure can be taken and the acceleration and stresses term of the third momentum equation vanish keeping only the pressure and the gravity terms. In this form the set of governing equations (1). ( 2 ) , (3) and (4) is integrated in the vertical z direction from the bottom position -h to the water level n. But, in order to get some information about the vertical structure of the flow, the integration process is split into several terms:
by the introduction of intermediate integration limits (fig. 1) which have their position located by: Lk =
n-(e) k
k = 1, 2 ,
... b
The presence of 11 in the equation (8) means that the position of those intermediate level must move proportionnally with the water level, consequently the relative thickness of a layer is constant. An uniform distribution of those level has been chosen here; however, it is possible to change those relative position by any other relationship including the water level to get another distribution. It is important to note that the intermediate integration surfaces do not have any physical meaning and allow the mass exchange between them. This is specifically prescribed by the intermediate boundary condition (Robert 1983):
n = k+l
a
Y
362
Fig. 1 Definition of the vertical integration boundaries. where Un and Vn are mean velocity components in a layer n. At the extremity of the vertical domain boundary conditions must be defined using a streamline condition: atl
atl
at
w(-h)
=
-u(-h)
ah - v(-h) ax
ah
The introduction of those conditions in the int gration pro ess using th Leibniz rule, does not eliminate all boundary velocity terms at the intermediate level, the assumption that those terms are included in the internal stresses term is made (Robert 1983). Then the model is defined as:
- Continuity equation
363
- Movement equations +
+ U k aauk x+Vkay auk -fV
e [>+
+
( TxLk-1 -
a
T
xLk
U kavk T + Vk avk + fUk
k = 1
))]
+
g
=
ahkTxyk aY
0
-- -
Y
The vertical velocity component is computed from equation ( 9 ) . The formulation of the horizontal stresses will be discussed in section 3. 2 . 3 Finite Element Formulation
This short overview of the finite element formulation will point out the specific choices made to obtain the discrete form of the model. The integral form of equations ( 1 2 ) , (13) and (14) is obtained by a Galerkin weighted residual method.
As one can see from the model equations, two dimensional elements can be used by associating at each node the water level
rl
and the two horizontal
components of the mean velocity in each layer Uk and Vk. Two types of interpolation functions were then considered. Based on the work of Cochet 1979 on 2-D hydrodynamic model, a quadratic function for the velocities and a linear function for the water level were tested. The result shown that in a multilayer configuration some instabilities appeared fig. 2. Then a full quadratic function was used and improved considerably the performance of the model. Two full quadratic elements have been developed, a six nodes triangular and a nine nodes quadrilateral elements. Special attention must be kept in the discretization of movement equations (13) and (14). Since the interpolation functions applie only on discrete parameter associated with one layer and are zero for any other parameters, the assembly process is equivalent to the sommation present in those equations. Then the resulting set of algebraic equation will have as much equations as discrete variables. The resulting non linear system in non steady form can be expressed by:
[Ml{il + CK ($)I {$I
=
0
The solution method will be presented In section 4.
364
.
-
I
I
Quadratic element Quadratic-linear element
5
Fig. 2
Computed velocities in the center of a schematic estuary presented at Fig. 5.
3 HORIZONTAL STRESSES 3.1 Formulation
One of the most important factor in the general behaviour of a natural mass of water is the horizontal shear stress. The predominant one is the effect of the bottom roughness over the entire depth of water and in some cases, the effect of the wind at the free surface may be determinant. In such case, this specific stress component may be expressed by:
365 Tx L = CtPa wx IWI 1 0 i where Ct is a drag coefficient, pa the air density and W the component of wind xi velocity W in the xi direction. For the bottom friction term, ChLzy formulation has been chosen:
TxiLb
3 vx i Ivl
=
Here, in the case of a two dimensional model the Ch6zy coefficient C , takes in account the overall effect of the bottom roughness on the flow over all the thickness of the water layer. However, when this depth is discretized, this term will affect only the bottom layer and the coefficient C must be corrected in way such as the instantaneous discharge must stay the same whatever the number of layer. To ensure the transfer of the bottom roughness to the upper layer, an internal shear stress component must appear at the interface of each layer. This may be expressed by discrete form of the Prandtl internal shear stress formulation; in the x direction we have at the upper interface of a layer k : =
TXLk-1
hk
Auk 4 AU; + AVk2
where Auk and AV are the differences of the flow velocity between the layer k k and k-1 and II is the mixing length. Many authors (Wang 1977, Kokayashi 1980, Leenderstee 1975) use a similar formulation, but the II coefficient is considered as an experimental adjustable coefficient. 3.2 Mixing Length Parameter
As we do not really know what is the value of the mixing length, expect near the bottom where Prandtl assume that it may be proportional to the distance from the bottom, we must find an usable formulation of the mixing length up to the free surface. Montgomery 1943 suggests, from some experiments in particular flow that:
x
is the Von-Karman coefficient which is defined by x = 0.4. As shown on fig.
3 , this distribution presents a plane symetry on the z/h = 0,5 plane. This
366
involves the assumption that the distribution of the horizontal stress is similar near the free surface and near the bottom, which is not really satisfactory. To remedy such situation, Yalin 1977 suggests a distribution of the mixing length based on the assumption of a logarithmic velocity profile:
A s the Montgomery 1943 distribution, this one agrees with the Prandtl assump-
tion near the bottom:
a
=
x z
(21)
Fig. 3 Mixing length vertical distribution. The major advantages of the Yalin 1977 distribution are: a) the mixing length decreases less quickly near the free surface than near the bottom, assuming that the turbulent fluctuation in the neighbourhood of the surface are less damped than by the presence of a rigid boundary at the bottom. b) the horizontal shear stress becomes zero without the necessity of assuming a zero vertical velocity gradient at the surface. This formulation was retained and adapted to the present model; in a discrete form, for the lower interface of a layer k, equation ( 2 0 ) becomes: k'
=
x
hk (b
-
k)
(22)
367
and the corrected bottom friction factor can be written:
where Cr is the Ch6zy reference coefficient for a two dimensional model. In order to verify the good overall behaviour of this approach, a test was run on the model to simulate a steady state uniform turbulent flow with a five level configuration. The velocity profile obtained is presented at the figure 4 .
t
L
y'm
K
..I
-*-
Fig. 4
.
-
-.u, m/s
Velocity profile for a uniform steady flow.
4 SOLUTION METHOD
In steady state configuration, the model of the flow is dominated by a quadratic friction term. The Newton Raphson procedure must be chosen rather than a substitution scheme to solve the non linear system. In fact, it is obvious that the later will not converge because the quadratic term. The incremental form of the method is:
368
where the tangent matrix is:
The non-steady flow is based on an implicit Euler scheme with an iterative Newton-Raphson stabilization procedure:
These techniques need much CPU time, and some possibilities have been set to reduce it. The tangent matrix can be, with a little loss of accuracy, recomputed at the beginning of each time steps rather than each iteration. In some case, computing time reduction may reach a factor three. 5 APPLICATION Many numerical tests have been performed with this model. The simulation of
the tridimensional structure of the flow in a 180'
curved flume in steady
state has been compared with a scale model in laboratory. The dimensions used in the mathematical model were the same than for the physical model to prevent any scale effect on the results. The test bring us the confirmation of the rightness of the formulation chosen here. The results of this comparison can be found in Robert 1983.
5.1 Schematic estuary In the testing process of the model, one goal was to verify that the vertical behaviour of the simulation agrees with the results of a two dimensional vertical models (width averaged) particularly with the Sunderman 1976 formulation based on cubic interpolation function in the vertical direction. The test case can be described as a schematic estuary closed at one end and submitted at the other one to a sinusoidal tidal cycle of 1 meter of amplitude. The flume has a length of 92 km and an width of 4.6 km. The depth is 50 meters except for a central portion where the bottom raise up to a depth of 20 meters. The finite element grid is composed of 10 quadrilateral elements and three layers were specified. The results of some time steps of a complete cycle (fig. 5) show a good agreement with those presented by Sundermann 1976 ones despite the fact that the vertical discretisation of the present model is far less accurate. The time interval of computation was. in this case, 1200 sec. for a full implicit scheme. This case was an occasion to test the validity of the mixing length distribution. The same test was run with a constant interfacial friction factor
369
Fig. 5 Vertical flow structure in a schematic estuary.
and shown some irregularities of the vertical velocity profile as an erroneous position of the maximum velocity during the ebb flow. That situation was cured by the use of the distribution presented by eq. 20. 5 . 2 Wind induced circulation
The second example focuses on wind driven circulation. Some preliminary numerical tests agreed with the Leenderstee 1975 samples in the same condition. As for preceeding example the simulation was performed as a propagation problem; a constant 40 km/h s-W wind blowing over a 10 km x 5 km x 1Gm rectangular bassin. The result, shown in fig. 6 , gives the circulation over the five layers, the water surface elevation and a transversal cross section in the middle of the bassin after six minutes of real time. At this time, the steady state is not totally achieved and the free surface is still moving. The corresponding computing time for this non-linear non-steady 2 2 2 1 degrees of freedom was about 2h40 on a VAX 850 for a full implicit scheme. Nevertheless, some techniques as the computation of the global matrix only at the beginning of each time step can reduce the computing time to one third and the use of a vectorial computer may be usefull. 6 . CONCLUSION
One of the most important concluding remark of this work, is that the automatic adjustment of the interlayer friction factor by the use of a well chosen mixing length distribution gives excellent vertical behaviour of the model. This vertical pattern of the flow is then relatively accurate and the formulation is well adapted to the natural flows. The use of the model is made easy by the fact that the finite element grid is only two dimensions and consequently a coupling with a two dimensional model is provided. This two and three dimensions models should be used to predict the vertical structure of the flow only where it is necessary.
371
-0,6 mrn
d
-
\
.
,
L
I
r
-
Fig. 6
r
c
Y
c
e
c
-
L-
c
+ Wind driven c i r c u l a t i o n i n a rectangular bassin.
0,51
312
7. REFERENCES Cochet. J.F., 1979. Modglisation d'gcoulements stationnaires et non stationnaires par la mlthode des 6lLments finis. ThPse de Docteur-Inghieur, Universitg de Technologie de CompiPgne, France, 142 p. Kobayashi, M., Nakota, K. and Kawahara, M., 1980. A Three Dimensional MultiLeveled Finite Element Model for Density Current Analysis. In: Third Int. Conf. on Finite Element in Flows Problems, Vol. 11. Banff, Canada, pp. 80-92. Leenderstee, J.J., Liu, S.K., Alexander, R.C. and Nelson, A.B., 1975. Threedimension Model for estuaries and Coastal Seas, Vol 11, Aspects of Computation, R-1764-OWRT, The Rand Corporation, California, USA, 165 p. Montgomery, R.B., 1943. Generalization for Cylinder of Prandtl's Linear Assumption for Mixing Length. Annuals of the New York Academy of Sciences, XLIV: 89-103. Robert, J.L., 1983. ModElisation tridimensionnelle des Ccoulements a surface libre, permanents et non permanents, par la mgthode des Blgments finis. ThSse de doctorat, Universitg Laval, Qudbec, Canada, 233 p. Sunderman, J., 1976. Computation of Barotropic Tides by the Finite Element Method. In: Finite Element in Water Resources, First Int. Conf. Gray, Pinder and Brebbia ed., Princeton University, Pentech Press, London, U.K.. pp. 4.31-4.67. Wang, H.P., 1977. Multi Leveled Finite Element Hydrodynamic Model of Block Island Sound. In: Finite Element in Water Resources, First Int. Conf., Gray Pinder and Brebbia ed., Princeton University, Pentech Press, London, U.K. Yalin, M.S., 1977. Mechanics of Sediment Transport, Second ed., Pergamon Press, Oxford, U.K., 259 p.