(
1
)
Dahong Chen1,
Thomas Ackermann, Christian Jokiel and J¨¹rgen Köngeter
Institute of Hydraulic Engineering and Water Resources Management,
Aachen University of Technology, Kreuzherrenstr., 52056 Aachen, Germany
1Current address: Department of River Engineering, Wuhan University, Wuhan 430072, China
Abstract: Free surface flows are very common flow phenomena in the hydraulic engineering. In this paper, a newly developed model to simulate such flows is presented. To solve the free surface problems, the model applied the arbitrary Lagrangian-Eulerian (ALE) kinematic description of the fluid domain, in which the movement of the mesh nodes is independent on the fluid motion. In the model, a new algorithm is developed for the mesh rezoning. It has been found that the ALE formulation is very efficient and allows a more accurate treatment of free surfaces than the purely Eulerian method. Compared to the Lagragian method it allows the handling of greater distortions in the fluid motion.
Keywords: Navier-Stokes equations, incompressible viscous flow, projection methods, finite element, ALE method, free surface
Problems with free surfaces, such as flows in rivers and seas, gain more and more importance in hydraulic engineering. In comparison to problems with fixed boundaries, for such problems not only the usual variables, i.e. velocity and pressure etc., but also the location of the free surface has to be determined, making calculations more difficult.
To determine the location of the free surface, some methods have been studied in the last decades. First treatment of the free surface problem is reported by Harlow & Welch1 who used the marker and cell (MAC) method. This simple method is widely applied. Second method, the volume-of-fluid (VOF) method, can be seen as an enhancement of the MAC method. It uses an approach similar to the MAC method, but employs only one variable, i.e. the VOF function, to determine the status of each cell. It was first introduced by Hirt & Nichols2 and overcomes important disadvantages of the MAC, e.g. it requires less memory and a shorter computing time for the determination of the free surface. However, the numerical calculation of the VOF function must be performed in a special way because of its discontinuity. Another popular, but more simple method to determine the location of the free surface is the height function method. The only disadvantage is that problems with multiple free surfaces such as surface breaking and merging can not be treated. Nevertheless, it is applicable for most problems.
Basically, the MAC and the VOF method are connected with the governing equations in the Eulerian formulation. The calculation is based on a stationary mesh. In contrast to the Eulerian formulation, another fundamental formulation is the so-called Lagrangian formulation in which the calculation is performed on a mesh that moves with the fluid. The height function method is often associated with the Lagrangian formulation. One of the disadvantages of the Eulerian formulation is the difficult treatment of the free surface boundary. Using the Lagrangian formulation the free surface can be treated more easily, but unacceptable degrees of mesh distortion may occur and potential entanglements may lead to an unsuccessful calculation.
The arbitrary Eulerian-Lagrangian (ALE) method, which is first presented by Hirt et al3., combines the advantages of the two methods mentioned above and avoids their impairments. In the past years, this method is mainly used to solve large-deformation problems in the area of solid mechanics and fluid-structure problems. Researchers have been attracted by the advantages of this method and some studies can be found in hydraulic engineering.4-6
This paper presents the development of a new model to simulate two-dimensional incompressible viscous fluid flows with a free surface by using the ALE method.
The equations of conservation of momentum and mass for incompressible Newtonian fluids in the ALE formulation are given as follows:
where uk is the fluid velocity in the Eulerian co-ordinate system, wj is the mesh velocity, and fk is the volume force. P denotes the kinematic pressure (pressure p divided by density r). z denotes the referential co-ordinate system, and x the Eulerian co-ordinate system.
In general, there are four different types of boundaries, i.e. inlet, outlet, wall and free surface, to be considered the equations above. The boundary conditions can be specified as follows: (a) give the velocity for the inlet; (b) apply the ¡°no-slip¡° for the wall condition and (c) give the stress tensor sij as
(3)
for the outlet and free surface.
The movement of the free surface is described by the kinematic boundary condition as follows:
where
S is the free surface,
is the velocity on the free
surface and
is the referential velocity for
the free surface. If a height function is applied to equation (4) the
location of the free surface is determined by
where h is the height function.
The fluid domain is discretized by using a mesh of quadrilateral isoparametric elements with four nodes. The basic functions for the velocity are linear and C0 continuous while the pressure is assumed to be elementwise constant and C-1 continuous. If the Galerkin finite element method is applied to (1) and (2) a weak formulation is obtained:
where u is a vector of the velocities at the nodes, P is a vector of pressures in elements, M is the consistent mass matrix (CM matrix), D is the diffusion matrix, V is the advection matrix, and fv is a force vector that includes the BCs on both outlet and free surface. The matrixes are given as follows:
(8)
(9)
(10)
(11)
(12)
where ji and yj are the finite element basis functions for node i and for element j, respectively. To discretize the free surfaec, linear isoparametric elements are used. The basic functions for the height function are linear and C0 continuous. According to (6) and (7) the formulation for the height function is given by
with
(14)
(15)
where h is a vector of height functions at the nodes on the free surface, Msf and Vsf are the CM matrix and advective matrix, respectively. fi represents the finite element basis function for node i.
The velocity and pressure in the spatial discretized momentum and continue equations (6) and (7) are decoupled using the projection methods (projection 1 and projection 2 method) and for the time integration a semi-implicit scheme is used.7 The equation (13) is integrated in time by using an implicit scheme.
To solve the resulting systems of algebraic equations the diagonally scaled preconditioned conjugate gradient (CG) method is applied for the velocity, the direct Cholesky solver for the pressure, and the Gaussian elimination method for the height function. Furthermore, in each time step the velocity field is calculated first. Then it is corrected with respect to the conservation of mass and the pressure is calculated if necessary. Finally, the height function is calculated.
The mesh rezoning is a significant process for a
simulation by using the ALE method. A new mesh will be constructed after this
processing. Basically, the new mesh can be determined by giving either the mesh
velocity
or the fluid velocity
(in the referential co-ordinate
system) because there is a relation between these velocities:
The Lagrange-Euler matrix method8 and the mixed formulation9 are based on the fluid velocity in the referential co-ordinate system. The mesh velocity or mesh displacement is then obtained by solving equation (16). But there is one distinct difference between these two methods. Using the Lagrange-Euler matrix method the mesh velocity or the mesh displacement is calculated within the whole fluid domain, whereas the mesh velocity is only determined on free surface using the mixed formulation. Donea et al.10 determine the mesh velocity at one node by using an averaged mesh velocity of all neighbouring nodes in the last time step. To improve the mesh quality, the gained mesh velocity must be modified by using a correction term.
In the work presented here, a simple and efficient method to perform the mesh rezoning is developed. After the height function is obtained by solving equation (13) the position of the interior nodes is calculated. Using this algorithm unacceptable mesh distortion and potential entanglement will be avoided.11
Fig.1 Schematic view and boundary conditions
For purposes of convenience non-dimensional variables are used in the following sections if not stated otherwise. Therefore, the velocity is ui/U, time is (U/H)¡Át, etc., where U stands for a referential velocity and H for a referential length. In addition, two non-dimensional parameters, the Reynolds number Re = U¡ÁH/n and the Froude number Fr = U2/g¡ÁH, are used. Physically, Re represents the ratio of inertia to friction, and Fr represents the ratio of inertia to gravity.
Fig.1 shows the schematic overview of the geometry. The fluid is bounded by three walls. The wall at the left side moves upwards with a constant speed U. The others are fixed. Due to friction the fluid near the left wall moves up, too. Therefore, the free surface changes. For this calculation the Reynolds number is set to 1.0, and the Froude number is Fr = 0.01. The flow domain is covered by a 41x41 mesh.
(a)
(b)
Fig. 2 Comparison of the location of the free
surface at time
(a) t = 0.51 and
(b) t = 2.21
As initial condition, the fluid is in a steady state. The velocity at the left side is increased instantaneously to the final state U. Fig. 2 represents the free surface at time t = 0.51 and t = 2.21. The numerical results by Frederiksen & Watts12 are plotted in the same diagram. The comparison of the results reveals good agreement.
The following description deals with a more complicated free surface problem, open channel with a sill at the bottom. It does not only have a complex geometry but also have a more practical relevance to hydraulic engineering. Its geometry and boundary conditions are given in Fig. 3. The inlet is at the left side of the channel, the outlet at the right side. The length of the channel is 30H, where H is the water depth at the inlet. The slope of the channel bed is 10%. The sill is placed in 7.75H £ x £ 8.75H. It has a height of 0.25H and its width is H. For the numerical calculation the ¡°no-slip¡° boundary condition at the bottom of the channel is used. For the free surface the stress tensor is set to zero. At the outlet the stress tensor in normal direction is set equal to the static pressure and in tangential direction it is set to zero. The velocity distribution at the inlet is given, too (Fig. 3). To discretize the domain, the mesh has 3960 elements with 4100 nodes.
Fig. 3 Schematic view and boundary conditions
To generate an appropriate velocity profile for the inlet and to determine a suitable Froude number Fr for a given Reynolds number Re, the obstacle is neglected in the first calculation. This is necessary to minimize the influence of the velocity distribution at the inlet and the Froude number Fr for the final calculation with the obstacle which causes a change of the free surface. The calculated velocity profile and the corresponding Re and Fr ensure that the open channel has a constant water depth and a regular velocity profile for the flow configuration without sill.
Fig. 4
depicts that the greater change of the free surface will be close to the
obstacle. To study the influence of the Reynolds number Re on the velocity and
the free surface, two additional calculations for Re = 100 and
Re = 500 are performed. In Fig. 5 the
free surfaces for these Reynolds numbers are plotted. The calculations
demonstrate that the main change of the water elevation is located near the
obstacle, but the influence caused by the obstacle reaches far downstream. From Fig. 5, the
location of the free surface is affected by the Reynolds number Re and Froude
number Fr obviously.
Fig. 4 Location of the free surface in the open channel for Re = 1000 and Fr = 25.- - - line = simulation without obstacle
Fig. 5 Effects of Reynolds number Re and Froude number Fr on the location of the free surface
A newly developed model based on the ALE formulation for the partial flow domain is presented. It used a new rezoning algorithm for controlling the change of the mesh during the calculation. Comparisons with data from literature indicated that the method is suitable to represent the free surface flows. The ALE method has been found particularly useful for such flows.
References
[1] F. H. Harlow and J. E. Welch, 'Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface', The Physics of Fluids, 8, 2182-2189, 1965.
[2] C. W. Hirt and B. D. Nichols, 'Volume of fluid (VOF) method for the dynamics of free boundaries', J. Comput. Phys., 39, 201-225, 1981.
[3] C. W. Hirt, A. A. Amsden and J. L. Cook, 'An arbitrary Lagrangian-Eulerian computing method for all flow speeds', J. Comput. Phys., 14, 227-253, 1974.
[4] B. Ramaswamy, 'Numerical simulation of unsteady viscous free surface flow', J. Comput. Phys., 90, 396-430, 1990.
[5] A. Soulaimani, M. Fortin, G. Dhatt and Y. Ouellet, 'Finite element simulation of two- and three-dimensional free surface flows', Comput. Methods Appl. Mech. Eng., 86, 265-296, 1991.
[6] P. Szabo and O. Hassager, 'Simulation of free surfaces in 3-D with the arbitrary Lagrangian-Eulerian methods', Int. J. Num. Meth. Eng., 38, 717-734, 1995.
[7] P. M. Gresho, 'On the theory of semi-implicit projection methods for viscous incompressible flow and its implementation via a finite element method that also introduces a nearly consistent mass matrix. Part 1: theory', Int. J. Num. Meth. Fluids, 11, 587-620, 1990.
[8] T. J. R. Hughes, W. K. Liu and T. K. Zimmermann, 'Lagrangian-Eulerian Finite Element Formulation for Incompressible Viscous Flows', Comput. Methods Appl. Mech. Eng., 29, 329-349, 1981.
[9] A. Huerta and W. K. Liu, 'Viscous fluid structure interaction', J. Pressure Vessel Technol., 110, 15-21, 1988.
[10] J. Donea, S. Giuliani and J. P. Halleux, 'An arbitrary Lagrangian-Eulerian finite element method for transient dynamics fluid-structure interactions', Comput. Methods Appl. Mech. Eng., 33, 689-723, 1982.
[11] D. Chen, Numerische Simulation von Strömungsvorgängen mit der "Arbitrary Lagrangian-Eulerian Method" (ALE-Methode); Mitteilungen des Instituts f¨¹r Wasserbau und Wasserwirtschaft, RWTH Aachen; Heft 110, 1997.
[12] C. S. Frederiksen and A. M. Watts, 'Finite-Element method for time-dependent incompressible free surface flow', J. Comput. Phys., 39, 282-304, 1981.
