dr. Erik A. Toorman
Hydraulics Laboratory, Katholieke Universiteit Leuven
INTRODUCTION
This text presents a description of the software package FENST- 2D, developed by the author at the Hydraulics Laboratory of the Katholieke Universiteit Leuven (e.g. Toorman, 1992).
FENST is an acronym, standing for Finite Elements for Navier- Stokes, Sediment Transport and Turbulence. It solves the equations in a 2-dimensional plane, hence the extension 2D to the name.
The code is based on a tutorial code from the University of Swansea, developed by Taylor & Hughes (1981). It has been modified and extended in many ways.
For sediment transport modelling, not a multi-phase, but a continuous phase (or mixture) approach is employed. This implies that the hydrodynamics is solved for the mixture of water and suspended sediment.
EQUATIONS
FENST has a modular structure. Each module corresponds to a set of equations, which matrix is solved with the same solver.
Module 1: Hydrodynamics
1.1. Continuity of incompressible fluid
fraction
1.2. Generalized Navier-Stokes equations
Generalized refers to the possibility of working with a non-
constant viscosity, i.e. in the case of turbulent or non-
Newtonian laminar flow. In the latter case, the coordinate system
independent shear rate intensity, which is the second invariant
of the shear (or strain, or deformation) rate tensor, is required.
The Boussinesq approximation is not employed to account for
momentum transfer due to density changes.
Module 2: Turbulence closure
FENST uses the two-equation k-epsilon model. The standard high-Reynolds,
various low-Reynolds (Lam-Bremhorst, new DNS-data based, ...)
and a hybrid two-layer formulation are implemented.
At present, isotropic turbulence is assumed.
Equations solved:
2.1. Conservation of turbulent kinetic energy k
2.2. Conservation of energy dissipation rate.
Module 3: Free surface movement
Mass conservation for the free surface is expressed by the kinematic condition. This module is used together with the Mixed Euler-Lagrange method to update the deforming mesh (Toorman, 1993).
Module 4: Sediment transport
Mass conservation of the sediment is solved. Erosion and deposition is accounted for by imposing the net sediment flux at the boundaries.
Module 5: Thixotropy
Time-dependent rheological behaviour of a visco-plastic fluid (e.g. to allow the study of laminar flow of fluid mud) is implemented according to the model proposed by Toorman (1997).
NUMERICAL STABILIZATION
Advection-dominated problems may show spurious oscillations as soon as the Peclet number (the element Reynolds number) exceeds the value 1. Therefore, standard and modified upwind stabilization techniques have been implemented, i.e. the standard Petrov-Galerkin method (Heinrich, 1981) and various versions of the Streamline-Upwind/PG method (Codina, 1993). As the kinematic condition for the free surface is a pure advection equation, upwinding is necessary for this equation under all circumstances.
The high degree of non-linearity of the k-epsilon equations makes it difficult to obtain monotonous convergence of these equations. Stabilization is obtained using two methods: by adding self-eliminating artificial diffusion (SEAD) and applying a pseudo-time stepping procedure (Rocabado et al., 1998).
BOUNDARY CONDITIONS
Many different types of boundaries can be distinguished, all of which can be handled by FENST. Either natural (or Neumann) type boundary conditions (NBC) or essential (or Dirichlet) type boundary conditions (EBC) are required.
Normally all pressure and stress conditions are of type NBC and all velocity conditions are of type EBC. In practice, an imposed pressure or stress distribution is set equal to the total traction (i.e. surface stress) on the corresponding boundary.
The need for a pressure EBC is extremely rare. The exceptional cases are those where no pressure boundary condition is needed at all. A typical case as such is the well known driven cavity flow (where all boundary conditions are imposed velocities). In order to have a reference pressure, the pressure may be imposed only in a wall corner node where both velocity are zero. In any other case a pressure EBC replaces the continuity equation in that node, which may lead to an erroneous solution (or no solution at all).
1. Solid wall boundaries
Hydrodynamics: In general no-slip conditions are used, i.e. all velocity components zero (EBCs). In some cases, a slip condition is useful: only the velocity normal to the boundary is set to zero. A wall friction stress NBC is used instead of imposing the tangential velocity. This is particularly useful in the case of free surface flow, where the fluid should be able to move up and down along the bounding vertical or sloping walls. In the latter particular case, wall friction is often neglected (perfect slip condition).
Turbulence: For turbulent flows solved with the high-Reynolds model, besides the traditional boundary conditions (which skips the visocus sublayer and employs wall functions), an unconventional method for the wall boundary conditions for k and epsilon has been implemented. All the equations are solved until the wall. Generally, at the wall u, v and k are set to 0 (EBC). The law of the wall is used to impose k and epsilon (EBC) at the first node away from each wall node, except when the wall node is a corner node. Epsilon at the wall is calculated. This yields excellent results for the few tested cases (pipe- and channel flow and backward-facing step).
For the low-Reynolds models, k and the normal gradient of epsilon are set zero at the wall.
Sediment transport: flux NBC, computed using the net flux of deposition (settling velocity x concentration) and erosion (e.g. Partheniades' formula), or in some exceptional cases an imposed concentration (e.g. at the bottom in the case of steady state open channel flow solved time-independent).
2. Inlet/outlet boundaries
Hydrodynamics: Different types of boundary conditions can occur, depending on the testcase. Inlet boundary conditions can be an imposed pressure distribution (NBC) or velocity profile (EBC).
Sediment transport: known concentration or flux profile.
Turbulence: imposed k and epsilon profile. In practice, these are often not known and must be approximated. A new method, which estimates initial values based on the velocity field of a preliminary laminar flow computation, has been implemented.
3. Open boundaries
This are boundaries with unknown conditions. By lack of better proposals, most often zero gradient NBCs are employed for all unknowns, even in cases where this is not true (e.g. backward facing step problem with short domain).
4. Symmetry plane
The mesh should be defined such that a symmetry plane is parallel to one of the coordinate axes. The boundary conditions are: zero velocity EBC perpendicular to the symmetry plane; zero gradient NBCs for all other unknowns.
5. Free surface
Hydrodynamics: a traction NBC is employed. For a flat surface, the traction is set equal to the surficial stress (e.g. wind stress). For a free surface without surficial stress, the traction is set equal to zero.
Turbulence: no good solution method exists. A opular method is to apply a zero-gradient NBC for k and compute epsilon assuming a reduced turbulent length scale, which in the case of openchannel flow is taken as a fraction of the water depth. In general this is not useful as the water depth may vary. It is preferred to apply some sort of wall-equivalent method by using the surficial shear stress from the computed velocity profile.
Sediment transport: the flux through the surface is imposed. For most cases this flux is zero, but not in cases such as mud dumping.
CLOSURE EQUATIONS
Closure equations can be modified, as we have developed the code
ourselves. Closure equations are used for the determination of:
- the settling rate of particles,
- erosion law (e.g. Partheniades' formula),
- turbulent Schmidt number,
- turbulence damping function (e.g. Munk-Anderson),
- rheological equation.
FURTHER DESCRIPTION
Element type: FENST uses 9-noded quadrilateral elements with mixed interpolation functions (Q2/P1 Lagragean): quadratic for u, v, k, epsilon and concentration, linear for p. These elements are second-order accurate in space.
Coordinate systems: FENST can use cartesian and cylindrical coordinates in an Eulerian frame. Mixed Euler-Lagrange formulations are possible to allow deformations of the mesh in special cases, such as moving interfaces, e.g. a free surface (Toorman, 1993).
Time scheme: The equations are solved implicitly and fully coupled. A first-order implicit finite difference discretization in time is implemented.
Numerical integration: The integration over the domain is done numerically using a 9-point Gauss rule.
Linearisation of the k-epsilon equations: The high degree of non-linearity leaves a great number of possibilities to linearize these equations. Based on intercomparison of the convergence behaviour of the different forms using turbulent open-channel flow as test-case the best scheme has been withheld.
Iterative procedure: Because of the non-linear nature of most of the equations and the coupling of the different sets of equations, an interative procedure is necessary. For sediment- laden turbulent flow, it solves the hydrodynamic, turbulence and sediment transport equations subsequently, before going to the next iteration step. Within each iteration step, the turbulence equations are solved iteratively until a certain minimal residual is reached, before continuing. This is necessary to overcome the intial bad convergence due to the difficulty to define good starting values for k and epsilon.
Solver: The sparse matrix is solved with a frontal solver using Gauss-Seidel elimination, which is a direct solver. This solver is relatively slow and restricts the size of the problem. A GMRES based iterative solver is under investigation, in cooperation with the Computer Science Department (Sedaghat et al., 1998).