Background Information

Program Features

Free and Moving Boundary Capabilites

Goma is a general purpose program designed for the solution of both steady and transient, two and three-dimensional problems involving heat, mass, and momentum (solid and fluid) transport. A unique feature is the treatment of all boundaries and interfaces as free (position unknown) or moving (position unknown or prescribed, but variable). If the material domain of interest is a solid, a Lagrangian formulation (i.e., the computational mesh follows the motion of material) of the momentum equations naturally leads to mass conservation and a natural parameterization of the boundaries and interfaces as material surfaces. If the material domain of interest is a fluid, then an Arbitrary-Lagrangian-Eulerian (ALE) formulation allows the boundaries to respond to constraint equations, hereafter referred to as distinguishing conditions. These conditions are responsible for determining the location of all boundaries and interfaces, providing the necessary mathematical closure of the system of equations governing the free boundary problem. Distinguishing conditions available to the user fall into several classes, as described below.

Since publication of the Goma 2.0 manual in 1998 (and more recently the Goma 4.0 manual in 2002), the ALE formulation has been extended to solid-material regions (viz. the TALE algorithm, Schunk, 2000) and purely Eulerian front tracking schemes based on the method of level-sets have been incorporated for free surfaces with large deformations; moreover, both algorithms have been implemented in a completely-coupled way. Of course Eulerian schemes are inherently transient and less accurate in capturing interfacial physics, even though they are more robust and even optimal for a certain class of problems. It is fair to say that of all the available mechanics codes, Goma provides the greatest breadth of free and moving boundary tracking formulations and options.

With regard to the ALE algorithms, the fully-implicit, pseudo-solid, unstructured mesh deformation algorithm sets Goma apart from other finite element programs. All surfaces, internal and external, together with other geometric features such as corners and junction points, are permitted to move as part of the algorithm. The movement of boundaries, interfaces, and geometric features is dictated by a weighted residual statement of the distinguishing conditions, whether based on specific physical constraints or arbitrary conditions described by the analyst. The internal mesh deforms as if it were embedded in a deforming elastic solid continuum; with the mechanics of the solid governed by either infinitesimal (linear) or finite (nonlinear) deformation theory. Through Newton’s method, the deformation of the mesh is determined simultaneously with all of the other physics of the problem.

The key connection between the mesh deformation and the physics of interest is accomplishedthrough a library of distinguishing conditions. Currently, those conditions include (a) kinematic (material surface of a fluid), (b) isotherm (phase transition temperature, such as melting), (c) isoconcentration and (d) geometric (either smooth plane curves or fixed point specifications). As part of the required input for Goma, the analyst specifies the associations between the particular distinguishing conditions and corresponding sets of material points of the initial pseudo-solid used to embody the mesh. Chapter 4 describes this process in more detail. Essentially, the algorithm causes smooth boundaries of the pseudo-solid to slide tangentially in a “frictionless” fashion. Further details of this algorithm and the corresponding equations can be found in several references (e.g., Sackinger, Schunk, and Rao, 1995).

Coordinate Systems and Frames of Reference

Coordinate systems accessible through this version of Goma include two-dimensional and threedimensional Cartesian coordinates, cylindrical coordinates for axisymmetric problems, spherical coordinates, and a swirling option for two-dimensional axisymmetric problems with a (swirling) velocity component in the third dimension. A limited framework has been built within Goma to use arbitrary orthogonal curvilinear coordinate systems, but this has not yet been extensively tested. As for frame of reference, all conservation equations are cast in an inertial frame (viz. nonaccelerating) but with extensions to allow for arbitrary frame velocities that may or may not be related to the material motion. Hereafter, when we refer to the frame/mesh motion type to be of the Eulerian variety, we mean the mesh is fixed with respect to all material motion, which basically means it is fixed in the laboratory frame. For now, we allow this frame of reference for fluid systems and are researching ways to allow this frame for solid systems. The ALE frame of reference, as mentioned above, allows for independent mesh motion in the interior of the domain, but seeks to maintain a material frame of reference on the boundary. This means that the mesh will move to accommodate material boundary motion. Currently, the ALE frame is allowed for all classes of materials (cf. Schunk, 2000). Finally, a pure Lagrangian frame of reference implies that our mesh moves with the material. This formulation is quite common in solid mechanics and is one advocated here for truly solid regions.

Problem Physics and Thermophysical Properties

This brief section summarizes the physics capabilities in Goma and the thermophysical properties and constitutive equations available to the user. The rest of the manual is designed to greatly expand on all material parameter options, boundary condition options, and equation options; perusing Chapter 4 and Chapter 5 is recommended to extract more detail.

The class of problems treated by Goma are those described by any one or a combination of the incompressible form of the momentum conservation equation for generalized Newtonian fluids, the momentum conservation and differential stress constitutive equations for viscoelastic fluids, saturated and unsaturated flow equations cast for rigid or deformable porous media, the energy conservation equation, the equations of quasi-static equilibrium of an elastic solid, and any number of additional or auxiliary species convection-diffusion-reaction equations. Goma has been tested with the following types of fluid mechanics, solid mechanics, and heat transfer problems: (a) mixed convection with mesh parameterization of an isotherm, (b) melting, with a parameterization of the liquidus and solidus isotherms, (c) coating and related flows (slide coating, curtain coating, etc.), (d) polymer processing (viscoelastic) flows (e.g. fountain flow, planar and axisymmetric extrusion, simple mold filling, contraction flow), (e) neutral or charged species transport in multicomponent concentrated systems, (f) partially saturated flow in poroelastic systems, (g) suspension flows, (h) drying and shrinking of gelled polymer films (with creep and elastic recovery), and (i) microfluidic systems with fluid-structure interaction (e.g. MEMS device performance).

Thermophysical properties in the bulk for all equations may be taken as constant or variable, with dependencies on any of the dependent and independent variables of the problem. General property variation models of this sort can be implemented with a user-defined subroutine capability. Moreover, a growing number of often-used standard models are supported within the core routines. These include a Carreau-Yasuda model for the generalized Newtonian viscosity and a Boussinesq source term for the fluid momentum equation that provides a means for simulating flows with thermal and solutal buoyancy. A plethora of other constitutive models and properties are available, including viscoelasticity, elastoviscoplasticity, nonFickian diffusivity, etc.

To enhance the capability for modeling problems in capillary hydrodynamics, e.g., coating flows, a boundary condition expressing the normal stress balance for two-dimensional Cartesian and axisymmetric problems has been implemented and tested. When capillary forces are activated, a pressure jump term (proportional to the mean curvature) is added to the normal component of the momentum flux balance at specified fluid material interfaces in a natural fashion. At three-phase boundaries (points in two dimensions) a contact angle condition and a surface tangent force condition may be applied. The former is used in place of a specified position on the mesh motion equations and is best used to set static and dynamic contact angles, and the latter is an additional endpoint force which is added to the momentum balance, necessitated because the curvature term is integrated by parts. The current version of Goma also includes the ability to model tangential shear forces along capillary surfaces, i.e., those originating from surface tension gradients caused, for example, by variations in temperature or species concentration. To access this capability requires a constitutive equation for the surface tension. A powerful low-level capability has been implemented which allows the user to select which degree of freedom, or variable, is associated with a particular boundary condition. Such a capability is useful at dynamic contact lines, where it is often desirable to replace the liquid-phase momentum equations with auxiliary constraint conditions.

Generalized interphase boundary conditions that allow for discontinuous field variables are supported through a multiple degree-of-freedom capability. The prime targets for this capability include flowing vapor-liquid equilibrium problems for which there are concentration and velocity jumps between phases due to change in density and solute partitioning through the phase diagram and multiphase/multicomponent corrosion problems. A series of boundary conditions which allow for the application of ideal and non-ideal vapor/liquid equilibrium (e.g. Raoult’s law and Flory-Huggins theory), latent heat release/adsorption, and discontinuous velocity components due to evaporation/condensation have been implemented. In the future this capability can be extended to thermal contact resistance, which often involves a temperature jump at an interface.

Recently the solid mechanics module of Goma, which was originally installed as a part of the pseudo-solid ALE mesh motion algorithm, has been exploited to solve problems in transport in deformable porous media and other outstanding problems of elastohydrodynamics. For modeling flow in non-deformable porous media, the Brinkman terms in the fluid momentum equations (cf. Gartling, et. al., 1996) may be activated. Since Goma 2.0, generalized Darcy transport equations for multiphase components (solid, liquid, gas) have been added and can be used for simulations of deformable poroelastic media. For incompressible but deformable solids, a pressure term was added to the solid momentum balance (e.g. rubber). In continuous shrinking or swelling solids, the dilation is proportional to changes in solvent concentration. In deformable porous media, the solid deformation is coupled to the pressure in the fluid-filled interstices of the porous matrix. Several boundary conditions exist to apply normal tractions (i.e. compressive, tensile, or shear boundary forces) to solid surfaces. To effectively simulate coupled fluid/solid interaction problems, boundary conditions which balance the surface tractions exerted by the liquid and solid phases at the common interface have been incorporated as have been the appropriate interface impregnation/expulsion conditions at boundaries between porous and continuous media.

A complete rewrite of the species transport equations has been undertaken since the release of Goma 2.0 that allows for generalized phase/species formulations on multimaterial problems. Accommodating an arbitrary number of species, each of which can exist in an arbitrary number of phases, was the goal of this development in order to model corrosion and charged species transport.

Of course there are many more material property models and constitutive equations, specialized boundary conditions, and more esoteric differential equations that can be solved for just about any mechanics problem. Many of these capabilities are not cited in this manual because they were under development at the time of publication. Interested readers should inquire about the status of the following capabilities: generalized solid-model geometry features, wetting and spreading models for Eulerian front tracking schemes, Eulerian/Eulerian fluid-structural interaction capability, multiphase porous energy equation, Generalized surface and volume user-defined Lagrange multiplier constraints, and much more.

Advanced Capabilities

Several developments in Goma that enable advanced engineering analysis of complex systems have been completed since the last major release. These developments include a complete, generalized capability of automated parameter continuation (zeroth-order, first-order, arclength, multiparameter, user-defined parameter continuation, etc.) using the LOCA library (Salinger, et. al., 2002), linear stability analysis of any dynamic system using normal modes, and augmenting condition capability. It is recommended that the user consult a separate manual (Gates et. al., 2000; contact authors for a more recent version) for a complete user description of these features. The input record sections required to activate these features are not covered in this document.

Numerical Methods

With over 150 different boundary conditions for 70 plus differential equation types, Goma’s algorithms are very extensive for any brief discussion. In this section we simply point out the foundation algorithms. A developer’s manual, advanced capabilities manual, and tutorial memos can be consulted for more details (see Goma Document List in the Appendix for the citations.).

Goma is based primarily on the Galerkin/finite element method. The element library currently includes (in two dimensions) 4- and 9-node isoparametric quadrilaterals (i.e., Q1 and Q2 interpolations) with available interpolations for linear discontinuous (P1) or piecewise constant (P0) variables, and (in three dimensions) 8-node isoparametric hexahedral elements and 27-node bricks, also available with piecewise constant interpolations. The overall solution algorithm centers around a fully-coupled Newton-Raphson iterative scheme for solving the nonlinear algebraic equations which results from the finite element discretization. That is, all active equations and boundary conditions are solved simultaneously in a single matrix system at the same time plane and during the same Newton iteration. The sparse matrix system is stored in a central element-level matrix data structure that is injected into one of three sparse matrix formats as dictated by the matrix solver chosen. The three formats are modified sparse row, MSR or compressed row format (Hutchinson, et. al., 1995, Schunk and Shadid, 1992), the variable block row, or VBR, format (see Heroux, 1992), or the frontal-solver element-level format (cf. Hood, 1976). If the matrix system is not too poorly conditioned, then iterative solvers of the generalized preconditioned conjugate gradient-type can be used to solve the system (see Tuminaro, et. al., 1999, Schunk and Shadid, 1992). A new matrix-services/solver-services library known as TRILINOS (http://www.cs.sandia.gov/Trilinos), has been installed to handle all iterative solver and preconditioner options. This package has greatly extended the robustness of iterative solvers to the class of problems that Goma solves. Virtually all problems and all finite element formulations are now solvable with these iterative schemes (see Schunk, et al., 2002). If all else fails, Goma deploys a suite of direct solvers that, even though not always efficient for large threedimensional problems, will always get a solution at the current Newton iteration. These solvers are known as Sparse 1.3 (lu), a classical LU decomposition (Gaussian elimination) method, and two frontal solvers, Umfpack (umf) and front; these are discussed in the next section.

The Galerkin least squares (GLS) method for pressure stabilization of Hughes and Franca (1987) has also been added to Goma. The GLS method adds the momentum residual, weighted by the gradient of the Galerkin weight function, to the standard Galerkin continuity equation, thus providing a diagonal term for the pressure. This is a first-order convergent and consistent method that enables the use of iterative solvers for incompressible equations over the entire range of Reynold’s numbers.

The overall differential-algebraic system of equations may be advanced in time with implicit time-integration techniques (simple backward Euler and Adams-Bashforth predictor, trapezoidal corrector algorithms for fluid systems, species transport and energy transport; and Newmark-Beta algorithms for solid dynamics). Time marching offers an alternative, albeit indirect, route to attaining solutions to steady equations, as well as providing the capability of simulating process transients directly. Automatic time step control based on current truncation error is also available.

Perhaps the most complicated part of the algorithm is the construction of the Jacobian sensitivity matrix. Because the mesh point positions are actually unknowns in a free or moving boundary problem, that matrix must include sensitivities of each weighted residual equation with respect to each of the mesh variable unknowns that can affect the value of the residual. Unfortunately, almost every term of the bulk equations and many boundary conditions contribute to this sensitivity. This occurs mainly through gradient operators and surface normal and tangent vectors (see Kistler and Scriven, 1983) and through dependencies on mesh position of the determinant of the elemental Jacobian transformation matrix that maps between a fixed unit element and any element in the computational domain. Great care has been taken to include analytical expressions for all of these mesh sensitivities. However, some of this task inevitably falls to the user when implementing user-defined boundary conditions, material property models, and constitutive equations, particularly when any of these quantities depends directly on spatial position or spatial gradients of other variables. In order to maintain the strong convergence properties of Newton’s method, these sensitivities must be specified in those user-defined routines. To aid in this task, a debugging option is available which computes a numerical finite-difference approximation of the global Jacobian matrix and compares it with its analytical counterpart. This tool enables users and developers to check the consistency of newly-created equations (whether bulk or boundary constraints) with their corresponding analytic Jacobian contributions.

Portability, Software Library Infrastructure, and Code Accessibility

Goma is written in the C programming language (specifically Kernighan and Ritchie, 1988, C with some ANSI extensions). It has been ported to a number of UNIX platforms including Solaris and Linux, with the Linux Enterprise-4 version being the most actively maintained. Most recent versions are aimed at Red-Hat RHEL5 and RHEL6 levels, almost exclusively. Many of the machine dependencies in the program have been isolated using C preprocessor directives. Some of the machine dependencies that occur in the I/O routines are insulated from the user by software libraries. Building Goma requires EXODUS II v2.02 (Schoof and Yarberry, 1994), SPARSE 1.3 (cf. Kundert and Sangiovanni-Vincentelli, 1988), NetCDF v2.3.2 (Rew, et. al., 1993) libraries, Umfpack direct solver libraries (Davis and Duff, 1997), and the TRILINOS 10.0 library (Tuminaro, et. al., 1999; http://software.sandia.gov/trilinos). The first of these is part of the SEACAS system at Sandia National Laboratories (Sjaardema, 1993); the latter two libraries are available publicly. Parallel processing is enabled by OPEN-MPI. The user should consult the build instructions for the most recent library revisitions. The most updated library needs are also made clear in the Goma makefile: Makefile. There are special versions of this makefile for building for the test suite (Makefile_guts) and debug mode (Makefile_debug). These are the most general makefiles that are deployed. Generally, pre- and post-processing is performed outside of Goma, although some post-processing of results is available within the program. This separation of the functionality permits the use of alternative solid-modeling and mesh-generation software and visualization packages of choice, insofar as they may be interfaced with the EXODUS II finite element data model.

Pre-processing options include mesh generation via CUBIT (http://cubit.sandia.gov), PATRAN (PDA, 1990), and SolidWorks (www.solidworks.com). The latter two require special plug-ins. These mesh generators currently support and will output a finite element database in the EXODUS II format.

Post-processing options include BLOT (see the SEACAS distribution, Gilkey and Glick, 1989), Paraview (www.paraview.org), and Ensight (www.mscsoftware.com.au/products/software/cei/ ensight).

Since Goma is built around the EXODUS II finite element data model, there are numerous options available for communication with other analysis codes that also exchange data via the same EXODUS II data model. Recent modifications to Goma permit not only the initialization of unknown values from an EXODUS II file, but also the ability to incorporate field variables into the analysis that are not unknowns. For example, the quasi-static and dynamic electromagnetic fields from codes such as ALEGRA can be used to compute electric fields and current fluxes on a specified finite element mesh that are input to Goma through the EXTERNAL FIELD data card.