|
FEFLOW Basic information
Physical concepts
Computational concepts
Other issues
PHYSICAL CONCEPTS
1. Code name:
FEFLOW - A finite-element simulation package for 2D and 3D density-dependent flow, mass and heat transport processes in
groundwater.
2. Advective-dispersive or sharp-interface:
Solving the advection-dispersion balance equations for contaminant mass (salinity) and heat in
two and three dimensions; also thermohaline flow.
3. Saturated or unsaturated aquifer:
Allowing the analysis of both saturated (groundwater) and unsaturated aquifers. For unsaturated
conditions, it solves the nonlinear Richards equation for optional parametric models such as van Genuchten-Mualem, Brooks-Corey and Haverkamp. Alternatively, for 3D regional unconfined aquifers, a multiple free
surface (perched water table) approach using moving meshes has been developed.
4. Multi or single phase:
Based on a single phase approach (salinity is considered as a miscible component).
5. Adsorption or conservative:
Adsorption effects can be modeled via Henry, Freundlich or Langmuir isotherms. For the conservation
equations of mass and heat, both divergence and convective formulations are available which lead to different boundary condition formulations in the form of total or dispersive boundary fluxes, respectively.
6. Transient or steady-flow conditions:
The flow can be either transient or steady-state.
7. Transient or steady-transport conditions:
The transport equations can by solved for either transient or steady-state conditions.
8. Isothermal or non-isothermal conditions:
The transport process can be considered as isothermal or non-isothermal processes; however,
this is optional. It can be considered that viscosity and fluid density are affected by salinity, and temperature can be considered. Boussinesq and extended Boussinesq approximations are able to be selected. The
latter is appropriate for strong density coupling (large density contrasts). Furthermore, a variable density expansion for the temperature effects is capable of solving large temperature ranges and to tackle the 4
degrees Centigrade anomaly (cold groundwater, also under variable saturation). Salinity and heat can be simultaneously coupled as given for thermohaline flow (double diffusion convection DDC).
9. Multi or single component (species):
The conservation equation is solved for a single component governed by advection, diffusion
(Fick's law), mechanical dispersion (Bear-Scheidegger dispersion model), retardation (Henry, Freundlich or Langmuir isotherms), zero and first-order reaction terms.
10. Chemical reactions or non-reactive:
Within the context of single species, a first order nonequilibrium reaction, in the form of a
decay rate, and a zero order reaction term can be specified. Non-reactive processes represent a special case within the approach.
COMPUTATIONAL CONCEPTS
1. Dimensionality:
FEFLOW is capable of solving 2D and 3D problems. Among 2D problems, vertical plane, horizontal plane and axisymmetric
domains can be schematized.
2. Numerical method:
A Galerkin-based finite-element method for unstructured meshes is used. For convection-dominant transport problems,
upwind techniques are available if needed such as streamline upwinding and shock capturing. For 3D hexahedral and pentahedral elements of a trilinear and triparabolic approximation, for 2D quadrilateral and
triangular elements of a bilinear and biparabolic accuracy are available. Different time marching schemes can be chosen: Automatic time stepping schemes based on the Gresho-Lee-Sani (GLS) predictor-corrector time
integrator for second order (Adams-Bashforth/trapezoid rule) and first order (forward Euler/backward Euler) accuracy as well as Crank-Nicolson or fully implicit schemes at fixed (user-defined) time increments.
Nonlinearities are solved either by the Newton or Picard iteration methods. The error-controlled predictor-corrector techniques have proven to be very powerful for nonlinear problems such as saltwater intrusion,
geothermal and thermohaline flows. A fully adaptive approach is provided for 2D based on mesh refinement and derefinement (AMR) techniques to enhance the reliability of the numerical simulation. An entropy-based
posterior error estimator controls the spatial discretization by using optimality criteria given by Zienkiewics & Zhu and Onate & Bugeda. Consistent and lumped mass matrices are optional (latter is
recommended for unsaturated media).
3. Boundary conditions:
Dirichlet (1st kind), Neumann (2nd kind) and Cauchy-type (3rd kind) boundary conditions can be specified for
flow, mass and heat. A so-called 4th kind of boundary condition exists for singular (pumping or injection) wells. These boundary conditions can be arbitrarily placed at nodal points of a 2D and 3D mesh. All boundary
conditions (from 1st to 4th kind) can be specified either as steady-state or transient. Furthermore, for each boundary condition, constraint formulations can be combined. They represent limitations of boundary
conditions and result from the requirement that boundary conditions should only be valid as long as minimum and maximum bounds are satisfied (e.g., a Dirichlet-type concentration boundary condition should only be
valid if the flow enters the boundary portion while a natural boundary condition has to be set up as soon as the flow releases the boundary). The constraints can also be time-dependent (for instance, if modeling a
spatially and/or temporally moving boundary condition). For groundwater (saturated aquifer), free-surface and seepage face boundary conditions are provided. For Cauchy-type boundary conditions, the accompanying
transfer (leakage) coefficient differs between infiltration and exfiltration conditions which depend on the direction of flow through the boundary. The leakage coefficients can also be time-dependent. In the
convective form of the transport equations, the Neumann and Cauchy boundary fluxes are mimicked by dispersive (derivative) fluxes while the divergence form of the transport equations allow the description of a total
(dispersive plus) convective) flux rate. The latter represents a totally conserved formulation; however, it needs specific techniques and additional effort for boundary conditions at outflowing portions. This is
practically solved via a postprocessing balance analysis to split up the convective part in the boundary condition in the form of an additional 'reaction' term.
4. Specified internal conditions:
There are no spatial and temporal limitations in prescribing boundary conditions and their
corresponding constraints of a mesh. The boundary conditions are node-related while material parameters are element-related. Both nodal and elemental quantities can be transient if desired.
5. Stresses:
Rainfall or evaporation (areal groundwater recharge) are normally modeled by sink/source formulations. They may also be
time-dependent. Pumping or injection conditions of singular wells (occurred at selected nodal points) are described by a boundary condition of the 4th kind. Multi-layer pumping or injection wells are also
accommodated in 3D.
6. Parameter input:
All parameters are handled on an element level, i.e., they can differ from element to element, and it is possible to
consider those as steady-state or transient quantities (e.g., conductivity can be changed with time). Accordingly, parameter heterogeneity, zoning or constant settings are arbitrarily possible. The input and
assignment of the parameter (including such distributed parameters like capillary pressure or relative conductivity curves with all related fitting coefficients) is done via a graphical problem attribute editor.
This attribution encompasses different graphical assignment methods; also interpolation techniques such as Akima, kriging or inverse distance or joining techniques for topological data related to the measured point
data as a primary input. Otherwise, there is a programming interface where the user can implement his own rules for parameter assignment and generation (e.g., random fields, alternative kriging techniques). The
programming interface is more general. It also allows manipulation of parameters during the simulation run which is appropriate in coupling a parameter identification routine. Anisotropy of the hydraulic
conductivity is supported for both 2D and 3D. In 3D, a specific technique is used to compute the overall rotation matrix from the stratigraphic layer structure by shape-derived principal directions of the 3D
prismatic finite elements.
7. The kind of equation solvers employed:
The solution of sparse equation systems can be solved either by iterative or direct techniques.
As default and especially for large systems (normally in 3D with more than hundred thousands of unknowns), the following methods are available: For symmetric equations, the preconditioned conjugate gradient (PCG) is
provided. Standard preconditioner such as the incomplete factorization (IF) technique and alternatively a modified incomplete factorization (MIF) technique based on the Gustafsson algorithm is used. Different
alternatives are available for the CG-like solution of the nonsymmetric transport equations: a restarted ORTHOMIN (orthogonalization-minimization) method, a restarted GMRES (generalized minimal residual) technique
and Lanczos-type methods such as CGS (conjugate gradient square), BiCGSTAB (bi-conjugate gradient stable) and BiCGSTABP (postconditioned bi-conjugate gradient stable). For preconditioning, the incomplete Crout
decomposition scheme is applied. For small (2D) or ill-posed problems, direct Gaussian elimination techniques for symmetric and nonsymmetric matrices are available. Here, the Reverse Cuthill-McKee (RCM) nodal
reordering technique is used to minimize fill-in entries of the matrices. The memory monitoring of the program is fully dynamic, accordingly, there are no software limits.
8. Nonlinearities:
Nonlinearities occur for density and viscosity coupling phenomena, unsaturated problems, free surface conditions,
boundary constrained conditions and/or nonlinear adsorption. In dependence on the used temporal discretization technique ((predictor-corrector versus marching with fixed predefined time steps), the Newton or Picard
iterative techniques are employed. The full Newton method is preferred for transport equations as soon as adaptive techniques are applied. This is normally the first choice. For the Richards equation (unsaturated
conditions), the Newton method is dropped to maintain the symmetric properties of the resulting flow equation systems. For fixed time stepping or steady-state approaches, the Picard iteration method is used in
combination with relaxation of primary variables.
Derivative fields, as occurred in computing the Darcy fluxes from potential results, are computed by either local or global smoothing techniques to
get continuous velocity distributions of a higher order accuracy. Such continuous velocity representations have shown to be superior to discontinuous approximations, especially for buoyancy-affected transport
problems.
OTHER ISSUES
1. Applicable graphical interfaces if they exist:
FEFLOW is an interactive, fully graphics-based commercial software system that contains
graphical editors and mesh generators for the geometric design and discretization as well as editors for the problem data attribution. It also supports a wide class of subsurface flow, mass and heat transport
problems, provides an open data interface for importing, exporting (GIS interface) and programming (interface manager IFM) as well as encompassing many graphical tools in the postprocessing analysis of the results.
There is a layer configurator to generate, refine and reposition stratigraphic and layer subdivisions as well as to inherit data. Tools for refinement and derefinement of meshes exist. The GIS-based data interface
is supported for maps and databases in the form of the ARC/INFO ASCII exchange format, ESRI Shape file format, Image (TIFF) format, AutoCAD (DXF) format and the HPGL data exchange format. Plotting tools are
integrated. Besides the 2D plotting equipment (isolines, diagrams, etc.), 3D visualization techniques are incorporated such as volume and surface displays, rotation, translation, shading, 3D cursor, arbitrary cross
sections, fence displays, isosurfaces, map integration, 2D and 3D pathline computations (by Runge-Kutta or Pollock method), flow vector patterns, isochrones and others.
2. Platforms and operation systems on which code will run:
FEFLOW is available for UNIX and PC.
The supported UNIX platforms are:
- SGI (IRIX 5.3 and above)
- SUN (Solaris 2.4 and above)
- DEC (DIGITAL UNIX 4.0 and above)
For PC:
- Microsoft Windows and Microsoft NT 4.0 and higher for INTEL PC's
The code is written in C and C++ and currently encompasses more than one million code statements.
3. Method of verification (concerning saltwater intrusion):
The development of FEFLOW began in 1979. Since that time, FEFLOW has
been continuously developed, modified, extended and improved for different practical tasks and theoretical studies. FEFLOW was extensively tested, compared to analytical and numerical solutions, benchmarked with
experimental data and successfully applied to a series of practical studies and problems as referenced in the user's manual. With regard to saltwater intrusion and density-dependent problems, the provings and
benchmarks have been considered:
- Comparison in cellular free convection problems (Elder problem, Benard problem) for solute, thermal and thermohaline processes in 2D and
3D
- Salt dome problem (HYDROCOIN case 5 level 1)
- Henry problem (saltwater encroachment)
- Saltwater upconing below pumping wells (comparison with sharp-interface model results and experimental data) and development of interception
techniques (so-called 'saltwater prevention wells')
- Flushing of saltwater horizons (comparison with analytical solutions).
4. Perceived strengths and weaknesses of the model:
FEFLOW is very attractive in practice and research due to the advanced interactive
graphical working environment, the integration of powerful and modern numerical techniques together with tools highly useful in the 'daily' handling of data and computational results, the network-based
implementation and the open data interface concept. The features of the package even allow very complex and large problems to be solved where other programs fail. Otherwise, due to the interactive merits and the
easy-to-use program handling, it is also appropriate for teaching purposes to study physical phenomena and effects of numerical approaches. Practical modeling needs a number of reruns to find out optimal design
variables for remediation and interception purposes as well as to study sensitivities of parameters and boundary conditions. Such a feedback procedure is highly supported by the package and makes engineering work
cost-effective and shortens the time required for a project (the user concentrates on his actual problem and does not waste time with burdensome and error-sensitive data works). FEFLOW provides an online-help
system. It is continuously being further developed and improved. Service and support of the package are worldwide.
FEFLOW is no longer a small program. It requires more resources and computational power from the hardware to profit from the implemented features of
the software. Nevertheless, the design of the package follows the standard capabilities of graphics workstations and PC's available today and is not necessarily dependent on high-end or super computing technology.
|
