Solution Procedures

Overview of Numerical Scheme

Icepak will solve the governing integral equations for mass and momentum, and (when appropriate) for energy and other scalars such as turbulence. A control-volume-based technique is used that consists of:

• Division of the domain into discrete control volumes using a computational grid.
• Integration of the governing equations on the individual control volumes to construct algebraic equations for the discrete dependent variables (unknowns) such as velocities, pressure, temperature, and conserved scalars.
• Linearization of the discretized equations and solution of the resultant linear equation system to yield updated values of the dependent variables.

The governing equations are solved sequentially (that is, segregated from one another). Because the governing equations are non-linear (and coupled), several iterations of the solution loop must be performed before a converged solution is obtained. Each iteration consists of the steps illustrated in Overview of Solution Method and outlined below:

1. Fluid properties are updated, based on the current solution. (If the calculation has just begun, the fluid properties will be updated based on the initialized solution.)
2. The u, v and w momentum equations are each solved in turn using current values for pressure and face mass fluxes, in order to update the velocity field.
3. Because the velocities obtained in Step 2 may not satisfy the continuity equation locally, a "Poisson-type" equation for the pressure correction is derived from the continuity equation and the linearized momentum equations. This pressure correction equation is then solved to obtain the necessary corrections to the pressure and velocity fields and the face mass fluxes such that continuity is satisfied.
4. Where appropriate, equations for scalars such as turbulence, energy, and radiation are solved using the previously updated values of the other variables.
5. A check for convergence of the equation set is made.

These steps are continued until the convergence criteria are met.

Overview of Solution Method

Linearization

The discrete, non-linear governing equations are linearized to produce a system of equations for the dependent variables in every computational cell. The resultant linear system is then solved to yield an updated flow-field solution.

The manner in which the governing equations are linearized takes an "implicit" form with respect to the dependent variable (or set of variables) of interest. For a given variable, the unknown value in each cell is computed using a relation that includes both existing and unknown values from neighboring cells. Therefore each unknown will appear in more than one equation in the system, and these equations must be solved simultaneously to give the unknown quantities.

This will result in a system of linear equations with one equation for each cell in the domain. Because there is only one equation per cell, this is sometimes called a "scalar" system of equations. A point implicit (Gauss-Seidel) linear equation solver is used in conjunction with an algebraic multigrid (AMG) method to solve the resultant scalar system of equations for the dependent variable in each cell. For example, the x-momentum equation is linearized to produce a system of equations in which u velocity is the unknown. Simultaneous solution of this equation system (using the scalar AMG solver) yields an updated u-velocity field.

In summary, Icepak solves for a single variable field (for example, p) by considering all cells at the same time. It then solves for the next variable field by again considering all cells at the same time, and so on.

Spatial Discretization

Icepak uses a control-volume-based technique to convert the governing equations to algebraic equations that can be solved numerically. This control volume technique consists of integrating the governing equations about each control volume, yielding discrete equations that conserve each quantity on a control-volume basis.

Discretization of the governing equations can be illustrated most easily by considering the steady-state conservation equation for transport of a scalar quantity φ. This is demonstrated by the following equation written in integral form for an arbitrary control volume v as follows:

Equation 106

where

ρ = density

= velocity vector (= in 2D)

= surface area vector

Γφ = diffusion coefficient for φ

φ = gradient of in 2D)

S_φ = source of φ per unit volume

Equation 106 is applied to each control volume, or cell, in the computational domain. The two-dimensional, triangular cell shown in Control Volume Used to Illustrate Discretization of a Scalar Transport Equation is an example of such a control volume. Discretization Equation 106 on a given cell yields

Equation 107

where

N faces = number of faces enclosing cell

φ_f = value of φ convected through face f

= mass flux through the face

= area of face f(A) (= in 2D))

(φ)_n = φ normal to face f

V = cell volume

The equations solved by Icepak take the same general form as the one given above and apply readily to multi-dimensional, unstructured meshes composed of arbitrary polyhedra.

Control Volume Used to Illustrate Discretization of a Scalar Transport Equation

Icepak stores discrete values of the scalar φ at the cell centers (c0 and c1 in Control Volume Used to Illustrate Discretization of a Scalar Transport Equation). However, face values φ_f are required for the convection terms in Equation 107 and must be interpolated from the cell center values. This is accomplished using an upwind scheme.

Upwinding means that the face value φ_f is derived from quantities in the cell upstream, or "upwind," relative to the direction of the normal velocity v_n in Equation 107. Icepak allows you to choose from two upwind schemes: first-order upwind, and second-order upwind. These schemes are described below.

The diffusion terms in Equation 107 are central-differenced and are always second-order accurate.

First-Order Upwind Scheme

When first-order accuracy is desired, quantities at cell faces are determined by assuming that the cell-center values of any field variable represent a cell-average value and hold throughout the entire cell; the face quantities are identical to the cell quantities. Thus when first-order upwinding is selected, the face value φ_f is set equal to the cell-center value of φ in the upstream cell.

Second-Order Upwind Scheme

When second-order accuracy is desired, quantities at cell faces are computed using a multidimensional linear reconstruction approach [T. J., Barth. The design and application of upwind schemes on unstructured meshes. Technical Report AIAA-89-0366. AIAA 27th Aerospace Sciences Meeting. Reno, Nevada:] . In this approach, higher-order accuracy is achieved at cell faces through a Taylor series expansion of the cell-centered solution about the cell centroid. Thus when second-order upwinding is selected, the face value φ_f is computed using the following expression:

Equation 108

where φ and ∇φ are the cell-centered value and its gradient in the upstream cell, and Δs is the displacement vector from the upstream cell centroid to the face centroid. This formulation requires the determination of the gradient ∇φin each cell. This gradient is computed using the divergence theorem, which in discrete form is written as

Equation 109

Here the face values are computed by averaging φ from the two cells adjacent to the face. Finally, the gradient ∇φ is limited so that no new maxima or minima are introduced.

Linearized Form of the Discrete Equation

The discretized scalar transport equation (Equation 107) contains the unknown scalar variable φ at the cell center as well as the unknown values in surrounding neighbor cells. This equation will, in general, be nonlinear with respect to these variables. A linearized form of Equation 107 can be written as

Equation 110

where the subscript nb refers to neighbor cells, and a_p and a_nb are the linearized coefficients for φ and φ_b.

The number of neighbors for each cell depends on the grid topology, but will typically equal the number of faces enclosing the cell (boundary cells being the exception).

Similar equations can be written for each cell in the grid. This results in a set of algebraic equations with a sparse coefficient matrix. For scalar equations, Ansys Icepak solves this linear system using a point implicit (Gauss-Seidel) linear equation solver in conjunction with an algebraic multigrid (AMG) method algebraic multigrid (AMG) solvermultigrid which is described in Restriction, Prolongation, and Coarse-Level Operators.

Under-Relaxation

Because of the nonlinearity stability of the equation set being solved by Icepak, it is necessary to control the change of φ. This is typically achieved by under-relaxation, which reduces the change of φ produced during each iteration. In a simple form, the new value of the variable φ within a cell depends upon the old value, φ_old , the computed change in φ, Δφ, and the under-relaxation factor, α, as follows:

Equation 111

Discretization of the Momentum and Continuity Equations

In this section, special practices related to the discretization of the momentum and continuity equations and their solution are addressed. These practices are most easily described by considering the steady-state continuity and momentum equations in integral form:

Equation 112

Equation 113

where is the identity matrix, is the stress tensor, and is the force vector.

Discretization of the Momentum Equation

The discretization scheme described earlier in this section for a scalar transport equation is also used to discretize the momentum equations. For example, the x-momentum equation can be obtained by setting φ = u:

Equation 114

If the pressure field and face mass fluxes were known, Equation 114 could be solved in the manner outlined earlier in this section, and a velocity field obtained. However, the pressure field and face mass fluxes are not known a priori and must be obtained as a part of the solution. There are important issues with respect to the storage of pressure and the discretization of the pressure gradient term; these are addressed later in this section.

Icepak uses co-located scheme, whereby pressure and velocity are both stored at cell centers. However, requires the value of the pressure at the face between cells c0 and c1, shown in Control Volume Used to Illustrate Discretization of a Scalar Transport Equation. Therefore, an interpolation scheme is required to compute the face values of pressure from the cell values.

Pressure Interpolation Schemes

The default pressure interpolation scheme in Icepak is the standard scheme. This scheme interpolates the pressure values at the faces using momentum equation coefficients . This procedure works well as long as the pressure variation between cell centers is smooth. When there are jumps or large gradients in the momentum source terms between control volumes, the pressure profile has a high gradient at the cell face, and cannot be interpolated using this scheme. If this scheme is used, the discrepancy shows up in overshoots/undershoots of cell velocity.

Flows for which the standard pressure interpolation scheme will have trouble include flows with large body forces, such as in strongly swirling flows and in high-Rayleigh-number natural convection. In such cases, it is necessary to pack the mesh in regions of high gradient to resolve the pressure variation adequately.

Another source of error is that Icepak assumes that the normal pressure gradient at the wall is zero. This is valid for boundary layers, but not in the presence of body forces or curvature. Again, the failure to correctly account for the wall pressure gradient is manifested in velocity vectors pointing in/out of walls.

The other scheme available in Ansys Icepak is the body-force-weighted scheme. This scheme computes the face pressure by assuming that the normal acceleration of the fluid resulting from the pressure gradient and body forces is continuous across each face. This works well if the body forces are known a priori in the momentum equations (for example, buoyancy and axisymmetric swirl calculations). This scheme is good for high-Rayleigh-number natural convection flows.

Discretization of the Continuity Equation

Equation 112 may be integrated over the control volume in Control Volume Used to Illustrate Discretization of a Scalar Transport Equation to yield the following discrete equation

Equation 115

where J_f is the mass flux through face f, ρv_n.

As described in Overview of Numerical Scheme, the momentum and continuity equations are solved sequentially. In this sequential procedure, the continuity equation is used as an equation for pressure. However, pressure does not appear explicitly in Equation 114 for incompressible flows, since density is not directly related to pressure. The SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithm  is used for introducing pressure into the continuity equation. This procedure is outlined below.

To proceed further, it is necessary to relate the face values of velocity v_n to the stored values of velocity at the cell centers. Linear interpolation of cell-centered velocities to the face results in unphysical checker-boarding of pressure. Icepak uses a procedure similar to that outlined by Rhie and Chow to prevent checkerboarding. The face value of velocity v_n is not averaged linearly; instead, momentum-weighted averaging, using weighting factors based on the a_P coefficient from , is performed. Using this procedure, the face flow rate J_f may be written as

Equation 116

where P_c0 and P_c1 are the pressures within the two cells on either side of the face, and Ĵ_f contains the influence of velocities in these cells (see Control Volume Used to Illustrate Discretization of a Scalar Transport Equation). The term d_f is a function of ā_P , the average of the momentum equation a_P coefficients for the cells on either side of face f.

Pressure-Velocity Coupling with SIMPLE

Pressure-velocity coupling is achieved by using Equation 114 to derive an equation for pressure from the discrete continuity equation (Equation 114). Icepak uses the SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) pressure-velocity coupling algorithm. The SIMPLE algorithm uses a relationship between velocity and pressure corrections to enforce mass conservation and to obtain the pressure field.

If the momentum equation is solved with a guessed pressure field p^∗ , the resulting face flux J_f ^∗ computed from Equation 114

Equation 117

does not satisfy the continuity equation. Consequently, a correction J_f ^′ is added to the face flow rate J_f ^∗ so that the corrected face flow rate J_f

Equation 118

satisfies the continuity equation. The SIMPLE algorithm postulates that J_f ^′ be written as

Equation 119

where is the cell pressure correction.

The SIMPLE algorithm substitutes the flux correction equations (Equation 118 and Equation 119) into the discrete continuity equation (Equation 114) to obtain a discrete equation for the pressure correction in the cell:

Equation 120

where the source term b is the net flow rate into the cell:

Equation 121

The pressure-correction equation (Equation 120) may be solved using the algebraic multigrid (AMG) method described in Restriction, Prolongation, and Coarse-Level Operators. Once a solution is obtained, the cell pressure and the face flow rate are corrected using

Equation 122

Equation 123

Here α _p is the under-relaxation factor for pressure (see Equation 113 and related description for information about under-relaxation). The corrected face flow rate J_f satisfies the discrete continuity equation identically during each iteration.

Coupled Pressure-Velocity Formulation

Selecting Coupled pressure-velocity formulation from the Solver options group box indicates that you are using the pressure-based coupled algorithm which is described in this section.

The pressure-based solver allows you to solve your flow problem in either a segregated or coupled manner. Using the coupled approach offers some advantages over the non-coupled or segregated approach. The coupled scheme obtains a robust and efficient single phase implementation for steady-state flows, with superior performance compared to the segregated SIMPLE solution scheme. For transient flows, using the coupled algorithm is necessary when the quality of the mesh is poor, or if large time steps are used.

The pressure-based segregated algorithm solves the momentum equation and pressure correction equations separately. This semi-implicit solution method results in slow convergence.

The coupled algorithm solves the momentum and pressure-based continuity equations together. The full implicit coupling is achieved through an implicit discretization of pressure gradient terms in the momentum equations, and an implicit discretization of the face mass flux, including the Rhie-Chow pressure dissipation terms.

In the momentum equations (Equation 114), the pressure gradient for component k is of the form

Equation 124

Where is the coefficient derived from the Gauss divergence theorem and coefficients of the pressure interpolation schemes. Finally, for any th cell, the discretized form of the momentum equation for component is defined as

Equation 125

In the continuity equation, the balance of fluxes is replaced using the flux expression in Equation 114, resulting in the discretized form

Equation 126

As a result, the overall system of equations (Equation 125 and Equation 126), after being transformed to the -form, is presented as

Equation 127

where the influence of a cell i on a cell j has the form

Equation 128

and the unknown and residual vectors have the form

Equation 129

Equation 130

Note Equation 127 that is solved using the coupled AMG.

Pseudo Transient Under-Relaxation

The pseudo transient under-relaxation method is a form of implicit under-relaxation. Here, the under-relaxation is controlled through the pseudo time step size. The pseudo time step size, which is automatically calculated by the solver, can be the same or different for different equations solved.

Equation 131

where is the pseudo time step.

Time Discretization

In Icepak the time-dependent equations must be discretized in both space and time. The spatial discretization for the time-dependent equations is identical to the steady-state case (see Spatial Discretization). Temporal discretization involves the integration of every term in the differential equations over a time step Δt. The integration of the transient terms is straightforward, as shown below.

A generic expression for the time evolution of a variable φ is given by

Equation 132

where the function F incorporates any spatial discretization. If the time derivative is discretized using backward differences, the first-order accurate temporal discretization is given by

Equation 133

where:

φ = a scalar quantity

n + 1 = value at the next time level, t + Δt

n = value at the current time level, t

Once the time derivative has been discretized, a choice remains for evaluating F(φ): in particular, which time level values of φ should be used in evaluating F?

One method (the method used in Icepak) is to evaluate F(φ) at the future time level:

Equation 134

This is referred to as implicit integration since φⁿ⁺¹ in a given cell is related to φⁿ⁺¹ in neighboring cells through F(φⁿ⁺¹):

Equation 135

until φⁱ stops changing (that is, converges). At that point, φⁿ⁺¹ is set to φⁱ .

The advantage of the fully implicit scheme is that it is unconditionally stable with respect to time step size.

Multigrid Method

Approach

Icepak uses a multigrid scheme to accelerate the convergence of the solver by computing corrections on a series of coarse grid levels. The use of this multigrid scheme can greatly reduce the number of iterations and the CPU time required to obtain a converged solution, particularly when your model contains a large number of control volumes.

The Need for Multigrid

Implicit solution of the linearized equations on unstructured meshes is complicated by the fact that there is no equivalent of the line-iterative methods that are commonly used on structured grids. Because direct matrix inversion is out of the question for realistic problems and whole-field solvers that rely on conjugate-gradient (CG) methods have robustness problems associated with them, the methods of choice are point implicit solvers like Gauss-Seidel method. Although the Gauss-Seidel scheme rapidly removes local (high-frequency) errors in the solution, global (low-frequency) errors are reduced at a rate inversely related to the grid size. Thus, for a large number of nodes, the solver stalls and the residual reduction rate becomes prohibitively low.

Multigrid techniques allow global error to be addressed by using a sequence of successively coarser meshes. This method is based upon the principle that global (low-frequency) error existing on a fine mesh can be represented on a coarse mesh where it again becomes accessible as local (high-frequency) error: because there are fewer coarse cells overall, the global corrections can be communicated more quickly between adjacent cells. Because computations can be performed at exponentially decaying expense in both CPU time and memory storage on coarser meshes, there is the potential for very efficient elimination of global error. The fine-grid relaxation scheme or smoother, in this case either the point-implicit Gauss-Seidel or the explicit multi-stage scheme, is not required to be particularly effective at reducing global error and can be tuned for efficient reduction of local error.

The Basic Concept in Multigrid

Consider the set of discretized linear (or linearized) equations given by

Equation 136

where φ_e is the exact solution. Before the solution has converged there will be a defect d associated with the approximate solution φ:

Equation 137

We seek a correction ψ to φ such that the exact solution is given by

Equation 138

Substituting Equation 138 into Equation 136 gives

Equation 139

Now using Equation 137 into Equation 139 we obtain

Equation 140

which is an equation for the correction in terms of the original fine level operator A and the defect d. Assuming the local (high-frequency) errors have been sufficiently damped by the relaxation scheme on the fine level, the correction ψ will be smooth and therefore more effectively solved on the next coarser level.

Restriction and Prolongation

Solving for corrections on the coarse level requires transferring the defect down from the fine level (restriction), computing corrections, and then transferring the corrections back up from the coarse level (prolongation). We can write the equations for coarse level corrections ψ^H as

Equation 141

where A^H is the coarse level operator and R the restriction operator responsible for transferring the fine level defect down to the coarse level. Solution of Equation 141 is followed by an update of the fine level solution given by

Equation 142

where P is the prolongation operator used to transfer the coarse level corrections up to the fine level.

Unstructured Multigrid

The primary difficulty with using multigrid on unstructured grids is the creation and use of the coarse grid hierarchy. On a structured grid, the coarse grids can be formed simply by removing every other grid line from the fine grids and the prolongation and restriction operators are simple to formulate (for example, injection and bilinear interpolation).

Multigrid Cycles

A multigrid cycle can be defined as a recursive procedure that is applied at each grid level as it moves through the grid hierarchy. Four types of multigrid cycles are available in Icepak: the V, W, F, and flexible (flex) cycles.

The V and W Cycles

V-Cycle Multigrid and W-Cycle Multigrid show the V and W multigrid cycles (defined below). In each figure, the multigrid cycle is represented by a square, and then expanded to show the individual steps that are performed within the cycle. You may want to follow along in the figures as you read the steps below.

V-Cycle Multigrid

W-Cycle Multigrid

For the V and W cycles, the traversal of the hierarchy is governed by three parameters, β₁ , β₂ , and β₃ , as follows:

1. β₁ " smoothings," (sometimes called pre-relaxation sweeps), are performed at the current grid level to reduce the high-frequency components of the error (local error).
   
   In V-Cycle Multigrid and W-Cycle Multigrid, this step is represented by a circle and marks the start of a multigrid cycle. The high-wave-number components of error should be reduced until the remaining error is expressible on the next coarser mesh without significant aliasing.
   
   If this is the coarsest grid level, then the multigrid cycle on this level is complete. (In V-Cycle Multigrid and W-Cycle Multigrid there are 3 coarse grid levels, so the square representing the multigrid cycle on level 3 is equivalent to a circle, as shown in the final diagram in each figure.)
   Note:

   In Icepak, β1 is zero (that is, pre-relaxation is not performed).

2. Next, the problem is "restricted" to the next coarser grid level using the appropriate restriction operator.
   
   In V-Cycle Multigrid and W-Cycle Multigrid, the restriction from a finer grid level to a coarser grid level is designated by a downward-sloping line.
3. The error on the coarse grid is reduced by performing β₂ multigrid cycles (represented in V-Cycle Multigrid and W-Cycle Multigrid and as squares). Commonly, for fixed multigrid strategies β₂ is either 1 or 2, corresponding to V-cycle and W-cycle multigrid, respectively.
4. Next, the cumulative correction computed on the coarse grid is "interpolated" back to the fine grid using the appropriate prolongation operator and added to the fine grid solution.
   
   In V-Cycle Multigrid and W-Cycle Multigrid, the prolongation is represented by an upward-sloping line.
   
   The high-frequency error now present at the fine grid level is due to the prolongation procedure used to transfer the correction.
5. In the final step, β₃ smoothings (post-relaxations) are performed to remove the high-frequency error introduced on the coarse grid by the β₂ multigrid cycles.
   
   In V-Cycle Multigrid and W-Cycle Multigrid, this relaxation relaxation scheme procedure is represented by a single triangle.

The Flexible Cycle

For the flexible cycle, the calculation and use of coarse grid corrections is controlled in the multigrid procedure by the logic illustrated in the following image. This logic ensures that coarser grid calculations are invoked when the rate of residual reduction multigrid solverresidual reduction rate on the current grid level is too slow. In addition, the multigrid controls dictate when the iterative solution of the correction on the current coarse grid level is sufficiently converged and should thus be applied to the solution on the next finer grid. These two decisions are controlled by the parameters α and β shown in the following image, as described in detail below. Note that the logic of the multigrid procedure is such that grid levels may be visited repeatedly during a single global iteration on an equation. For a set of 4 multigrid levels, referred to as 0, 1, 2, and 3, the flex-cycle multigrid procedure for solving a given transport equation might consist of visiting grid levels as 0-1-2-3-2-3-2-1-0-1-2-1-0, for example.

Logic Controlling the Flex Multigrid Cycle

The main difference between the flexible cycle and the V and W cycles is that the satisfaction of the residual reduction tolerance and termination criterion determine when and how often each level is visited in the flexible cycle, whereas in the V and W cycles the traversal pattern is explicitly defined.

The F Cycles

The multigrid F cycle is essentially a combination of the V and W cycles described in Multigrid Cycles.

Recall that the multigrid cycle is a recursive procedure. The procedure is expanded to the next coarsest grid level by performing a single multigrid cycle on the current level. Referring to V-Cycle Multigrid and W-Cycle Multigrid, this means replacing the square on the current level (representing a single cycle) with the procedure shown for the 0–1 level cycle (the second diagram in each figure). We see that a V cycle consists of:

and a W cycle:

An F cycle is formed by a W cycle followed by a V cycle:

As expected, the F cycle requires more computation than the V cycle, but less than the W cycle. However, its convergence properties turn out to be better than the V cycle and roughly equivalent to the W cycle. The F cycle is the default AMG cycle for the coupled equation set and for the scalar energy equation.

The Residual Reduction Rate Criteria

The multigrid procedure invokes calculations on the next coarser grid level when the error reduction rate on the current level is insufficient, as defined by

Equation 143

Here R_i is the absolute sum of residuals (defect) computed on the current grid level after the relaxation on this level. The above equation states that if the residual present in the iterative solution after i relaxations is greater than some fraction, β (between 0 and 1), of the residual present after the (i − 1)th relaxation, the next coarser grid level should be visited. Thus β is referred to as the residual reduction tolerance, and determines when to give up on the iterative solution at the current grid level and move to solving the correction equations on the next coarser grid. The value of β controls the frequency with which coarser grid levels are visited. The default value is 0.1. A larger value will result in less frequent visits, and a smaller value will result in more frequent visits.

The Termination Criteria

Provided that the residual reduction rate is sufficiently rapid, the correction equations will be converged on the current grid level and the result applied to the solution field on the next finer grid level.

The correction equations on the current grid level are considered sufficiently converged when the error in the correction solution is reduced to some fraction, α (between 0 and 1), of the original error on this grid level:

Equation 144

Here, R_i is the residual on the current grid level after the ith iteration on this level, and R₀ is the residual that was initially obtained on this grid level at the current global iteration. The parameter α , referred to as the termination criterion, has a default value of 0.1. Note that the above equation is also used to terminate calculations on the lowest (finest) grid level during the multigrid procedure. Thus, relaxations are continued on each grid level (including the finest grid level) until the criterion of this equation is obeyed (or until a maximum number of relaxations has been completed, in the case that the specified criterion is never achieved).

Restriction, Prolongation, and Coarse-Level Operators

The multigrid algorithm in Ansys Icepak is referred to as an "algebraic" multigrid (AMG) scheme because, as we shall see, the coarse level equations are generated without the use of any geometry or re-discretization on the coarse levels; a feature that makes AMG particularly attractive for use on unstructured meshes. The advantage is that no coarse grids have to be constructed or stored, and no fluxes or source terms need be evaluated on the coarse levels. This approach is in contrast with FAS (sometimes called "geometric") multigrid in which a hierarchy of meshes is required and the discretized equations are evaluated on every level. In theory, the advantage of FAS over AMG is that the former should perform better for non-linear problems since non-linearities in the system are carried down to the coarse levels through the re-discretization; when using AMG, once the system is linearized, non-linearities are not felt by the solver until the fine level operator is next updated.

AMG Restriction and Prolongation Operators

The restriction and prolongation operators used here are based on the additive correction (AC) strategy described for structured grids by Hutchinson and Raithby . Inter-level transfer is accomplished by piecewise constant interpolation and prolongation. The defect in any coarse level cell is given by the sum of those from the fine level cells it contains, while fine level corrections are obtained by injection of coarse level values. In this manner the prolongation operator is given by the transpose of the restriction operator

Equation 145

The restriction operator is defined by a coarsening or "grouping" of fine level cells into coarse level ones. In this process each fine level cell is grouped with one or more of its "strongest" neighbors, with a preference given to currently ungrouped neighbors. The algorithm attempts to collect cells into groups of fixed size, typically two or four, but any number can be specified. In the context of grouping, strongest refers to the neighbor j of the current cell i for which the coefficient A_ij is largest. For sets of coupled equations A_ij is a block matrix and the measure of its magnitude is simply taken to be the magnitude of its first element. In addition, the set of coupled equations for a given cell are treated together and not divided amongst different coarse cells. This results in the same coarsening for each equation in the system.

AMG Coarse Level Operator

he coarse level operator A^H is constructed using a Galerkin approach. Here we require that the defect associated with the corrected fine level solution must vanish when transferred back to the coarse level. Therefore we may write

Equation 146

Upon substituting Equation 137 and Equation 142 for d^new and φ^new we have

Equation 147

Now rearranging and using Equation 137 once again gives

Equation 148

Comparison of Equation 148 with Equation 137 leads to the following expression for the coarse level operator:

Equation 149

The construction of coarse level operators thus reduces to a summation of diagonal and corresponding off-diagonal blocks for all fine level cells within a group to form the diagonal block of that group’s coarse cell.

Solution Residuals

During the solution process you can monitor the convergence dynamically by checking residuals. At the end of each solver iteration, the residual sum for each of the conserved variables is computed and stored, thus recording the convergence history. This history is also saved in the data file. The residual sum is defined below.

On a computer with infinite precision, these residuals will go to zero as the solution converges. On an actual computer, the residuals decay to some small value ("round-off") and then stop changing ("level out"). For single precision computations (the default for workstations and most computers), residuals can drop as many as six orders of magnitude before hitting round-off. Double precision residuals can drop up to twelve orders of magnitude.

After discretization, the conservation equation for a general variable φ at a cell P can be written as

Equation 150

Here a_P is the center coefficient, a_nb are the influence coefficients for the neighboring cells, and b is the contribution of the constant part of the source term Sc in S = S_c + S_Pφ and of the boundary conditions. In Equation 150,

Equation 151

The residual R^φ computed by Icepak is the imbalance in Equation 150 summed over all the computational cells P. This is referred to as the "unscaled" residual. It may be written as

Equation 152

In general, it is difficult to judge convergence by examining the residuals defined by Equation 152 since no scaling of is employed. This is especially true in enclosed flows such as natural convection in a room where there is no inlet flow rate of φ with which to compare the residual. Icepak scales the residual using a scaling factor representative of the flow rate of φ through the domain. This "scaled" residual is defined as

Equation 153

For the momentum equations the denominator term a_Pφ_P is replaced by a_P , where v_P is the magnitude of the velocity at cell P.

The scaled residual is a more appropriate indicator of convergence, and is the residual displayed by Icepak.

For the continuity equation, the unscaled residual is defined as

Equation 154

R^{c} = \sum _{{\rm cells}\ P}{|{\rm rate\ of\ mass\ creation\ in\ cell\ P}|} \label {def-resid-eq5}\mylabel {def-resid-eq5}

The scaled residual for the continuity equation is defined as

Equation 155

The denominator is the largest absolute value of the continuity residual in the first five iterations.