6. Temporal Discretization and Algebraic Solution Methods

Introduction and Context

The governing equations of fluid mechanics are continuous in space and time, whereas a computational fluid dynamics (CFD) program must represent them on a finite mesh and advance them through a finite sequence of computational operations. Spatial discretization converts differential equations into algebraic relations between neighboring control volumes or nodes. Temporal discretization performs the corresponding conversion in time: derivatives with respect to time are approximated using values at one or more discrete time levels.

This conversion is not merely a programming detail. The choice of time-integration method determines stability, accuracy, computational cost, memory requirements, and the extent to which a simulation can capture transient physical phenomena. Once the temporal derivative has been discretized, the resulting equations may be explicit, requiring only known quantities, or implicit, requiring the solution of a system of algebraic equations at the new time level. The latter leads directly to the second major concern of this chapter: efficient solution of sparse linear and nonlinear systems.

A successful unsteady CFD calculation therefore involves two related decisions. First, the numerical method must advance the solution in time without unacceptable truncation error or instability. Second, the algebraic equations generated at each step must be solved accurately and economically. These issues are especially important for viscous flows, compressible flows, reacting flows, multiphase systems, and flows containing widely separated time scales. In such problems, stiffness may make an apparently accurate explicit method impractical, while an implicit method may be limited by the cost or robustness of its linear solver.

Learning Objectives

After studying this chapter, the reader should be able to:

• distinguish explicit, implicit, single-step, and multistep time-integration methods;
• derive and interpret first- and second-order temporal approximations;
• explain CFL restrictions, stiffness, numerical damping, and temporal convergence;
• formulate the algebraic systems produced by implicit CFD methods;
• compare Jacobi, Gauss–Seidel, Krylov, Newton, multigrid, and preconditioned solvers;
• select convergence and stopping criteria using residuals and engineering quantities;
• diagnose common failures in transient and steady CFD calculations;
• evaluate practical trade-offs among accuracy, stability, memory, and computational cost.

From the Semi-Discrete Equations to a Time-Marching Problem

After spatial discretization, a broad class of conservation equations can be written in semi-discrete form as

\[

\mathbf{M}\frac{d\mathbf{U}}{dt}=\mathbf{R}(\mathbf{U},t),

\]

where \(\mathbf{U}\) is the vector of discrete unknowns, \(\mathbf{M}\) is a mass matrix or a diagonal matrix of control-volume measures, and \(\mathbf{R}\) is the spatial residual. In a finite-volume formulation, the residual is obtained from the net fluxes through cell faces and from source terms. For a steady solution, the desired state satisfies \(\mathbf{R}(\mathbf{U})=0\). For an unsteady solution, the residual determines the rate of change.

Dividing by the mass matrix, when appropriate, gives

\[

\frac{d\mathbf{U}}{dt}=\mathbf{F}(\mathbf{U},t).

\]

This form reveals the relationship between CFD time integration and ordinary differential equation (ODE) integration. Nevertheless, CFD systems are usually much larger, sparse, nonlinear, and more strongly coupled than the ODE systems encountered in elementary numerical analysis.

A time step \(\Delta t\) introduces discrete levels \(t^n\), \(t^{n+1}=t^n+\Delta t\), and so forth. The temporal truncation error depends on how the derivative is approximated. Stability, by contrast, concerns the behavior of numerical errors as the calculation progresses. A method can be formally high-order accurate but unusable if its stability region is incompatible with the physical and mesh time scales.

Explicit Time Integration

The forward Euler method approximates the derivative at the known time level:

\[

\frac{\mathbf{U}^{n+1}-\mathbf{U}^n}{\Delta t}=\mathbf{F}(\mathbf{U}^n,t^n),

\]

or

\[

\mathbf{U}^{n+1}=\mathbf{U}^n+\Delta t\,\mathbf{F}(\mathbf{U}^n,t^n).

\]

Because the right-hand side contains only known quantities, the method is explicit. Each cell can generally be updated without solving a global matrix system, although boundary conditions and coupled physics may still require communication or local operations.

The principal attraction of explicit methods is low cost per time step. They are relatively simple to parallelize and can be highly effective for wave propagation, compressible flows, and problems in which the allowable stability time step is comparable to the physical time scale of interest. Their central limitation is conditional stability. The time step must be sufficiently small relative to the fastest signal or diffusion process represented on the mesh.

For the one-dimensional advection equation

\[

\frac{\partial u}{\partial t}+a\frac{\partial u}{\partial x}=0,

\]

an explicit upwind discretization produces a Courant number

\[

C=\frac{|a|\Delta t}{\Delta x}.

\]

A typical stability requirement is \(C\leq 1\), although the exact limit depends on the spatial and temporal discretizations. In multiple dimensions, the relevant restriction commonly involves the sum of directional Courant numbers. For compressible flow, the wave speed includes the local fluid velocity and acoustic speed, so a useful estimate is based on \(|u|+c\), where \(c\) is the speed of sound.

Diffusion creates a different restriction. For the one-dimensional diffusion equation, an explicit central approximation often requires a condition of the form

\[

\frac{\nu\Delta t}{\Delta x^2}\leq C_d,

\]

where \(\nu\) is diffusivity and \(C_d\) depends on the discretization and dimensionality. Since the restriction scales with the square of the smallest cell size, highly refined meshes can make explicit methods prohibitively expensive for viscous flows.

Higher-order explicit Runge–Kutta methods improve temporal accuracy without abandoning the explicit structure. A two-stage method, for example, evaluates the residual at intermediate states before forming the final update. Strong-stability-preserving Runge–Kutta schemes are designed to retain useful monotonicity properties under appropriate time-step restrictions. Their stability regions are generally larger than that of forward Euler, but they remain conditionally stable.

Implicit and Semi-Implicit Integration

Backward Euler evaluates the residual at the new time level:

\[

\frac{\mathbf{U}^{n+1}-\mathbf{U}^{n}}{\Delta t}=\mathbf{F}(\mathbf{U}^{n+1},t^{n+1}).

\]

The method is first-order accurate in time and, for many linear dissipative systems, unconditionally stable in the linear stability sense. Unconditional stability does not mean that any time step is accurate. A very large step may strongly damp transient behavior and produce a solution that is stable but physically misleading.

The second-order backward differentiation formula (BDF2) is

\[

\frac{3\mathbf{U}^{n+1}-4\mathbf{U}^{n}+\mathbf{U}^{n-1}}{2\Delta t}

=\mathbf{F}(\mathbf{U}^{n+1},t^{n+1}).

\]

It uses two previous solution levels and is particularly common in stiff CFD applications. It is more accurate than backward Euler for smooth solutions, but it requires a starting procedure and retains an implicit nonlinear system at every step.

The trapezoidal or Crank–Nicolson method averages the residual between old and new levels:

\[

\frac{\mathbf{U}^{n+1}-\mathbf{U}^{n}}{\Delta t}

=\frac{1}{2}\left[\mathbf{F}(\mathbf{U}^{n+1})+\mathbf{F}(\mathbf{U}^{n})\right].

\]

It is second-order accurate and has relatively little numerical damping. This property is valuable for oscillatory phenomena but can also permit nonphysical oscillations when the solution contains sharp gradients or when the time step is large. A weighted family of methods, often called the theta method, uses

\[

\mathbf{F}^{n+\theta}=\theta\mathbf{F}(\mathbf{U}^{n+1})+(1-\theta)\mathbf{F}(\mathbf{U}^{n}),

\]

with backward Euler corresponding to \(\theta=1\) and Crank–Nicolson to \(\theta=1/2\). Choosing \(\theta>1/2\) introduces additional damping and can improve robustness.

In an implicit method, the new solution appears inside fluxes, material properties, source terms, and possibly turbulence or reaction models. Consequently, the update is not a direct formula but a nonlinear algebraic problem. Linearization produces a system such as

\[

\mathbf{A}\,\delta\mathbf{U}=\mathbf{b},

\]

where \(\delta\mathbf{U}\) is a correction to the current approximation. The matrix \(\mathbf{A}\) is sparse because a control volume interacts directly with only a limited number of neighboring volumes, but it may be nonsymmetric, poorly conditioned, and strongly coupled among variables.

Multistep Schemes, Accuracy, and Startup

A single-step method uses only the current state to construct the next state, whereas a multistep method uses several earlier states. Multistep methods can achieve higher order with relatively few residual evaluations per step. Their disadvantages include storage of previous solutions, sensitivity to variable time steps, and the need for a special startup method.

The order of a method describes how its local truncation error decreases as \(\Delta t\) becomes small. If the global temporal error behaves as

\[

E_t\approx K(\Delta t)^p,

\]

then \(p\) is the observed temporal order under suitable smoothness and stability conditions. A second-order method should approximately reduce the error by a factor of four when the time step is halved, provided spatial error, iterative error, and nonlinear error are already smaller than the temporal error.

Variable time stepping is useful when the flow evolves at different rates. A small step may be employed during rapid transients and a larger step during slow evolution. However, variable-step multistep schemes require modified coefficients and careful control of step-size ratios. Abrupt changes in \(\Delta t\) can reduce accuracy or affect stability.

CFL Restrictions, Stiffness, and the Meaning of Stability

The Courant–Friedrichs–Lewy (CFL) concept expresses a compatibility requirement between the numerical time step and the propagation of information. In simple advection, a numerical domain of dependence must contain the physical domain of dependence. In practical CFD, local CFL numbers may be estimated as

\[

\mathrm{CFL}_i\sim \frac{\lambda_i\Delta t}{\Delta x_i},

\]

where \(\lambda_i\) is an appropriate characteristic speed. For finite-volume methods, the estimate may be formed from the sum of face mass fluxes divided by the cell volume.

A system is stiff when some modes decay or evolve much faster than the modes of primary interest. Viscous diffusion on a fine mesh, chemical reactions, turbulence-model source terms, low-Mach pressure adjustment, and strongly coupled multiphysics can all introduce stiffness. Explicit methods must resolve the fastest stable mode even if the engineer is interested only in a slow mode. Implicit methods reduce this restriction by treating rapidly varying terms at the new time level, but they replace a small-step requirement with a more difficult algebraic solve.

Stability must be distinguished from accuracy and convergence. Stability concerns bounded propagation of numerical errors. Accuracy concerns approximation of the governing equations. Iterative convergence concerns how closely the algebraic equations are solved. A computation may be stable but inaccurate, formally accurate but insufficiently converged, or algebraically converged to a discretized model that is physically inappropriate.

Nonlinear Algebraic Systems and Linearization

Let the fully discrete equations at a new time level be written as

\[

\mathbf{G}(\mathbf{U}^{n+1})=\mathbf{0}.

\]

A common approach is to use Newton's method. Given an iterate \(\mathbf{U}^{(k)}\), linearization gives

\[

\mathbf{J}(\mathbf{U}^{(k)})\delta\mathbf{U}^{(k)}

=-\mathbf{G}(\mathbf{U}^{(k)}),

\]

followed by

\[

\mathbf{U}^{(k+1)}=\mathbf{U}^{(k)}+\alpha\delta\mathbf{U}^{(k)},

\]

where \(\mathbf{J}\) is the Jacobian and \(\alpha\) is a damping factor. Full Newton methods can converge rapidly near the solution, but forming and storing exact Jacobians may be expensive. CFD codes therefore often use approximate Jacobians, defect-correction methods, Picard or fixed-point linearization, or pseudo-time iterations.

Under-relaxation replaces a full correction by a fraction of it. It can stabilize nonlinear iterations when the equations are strongly coupled, but excessive under-relaxation increases computational cost and may conceal a poor discretization or boundary-condition problem. A robust solver strategy often begins with conservative relaxation and gradually becomes more aggressive as the solution improves.

Sparse Linear Solvers

The linear system \(\mathbf{A}\mathbf{x}=\mathbf{b}\) generated by a CFD discretization is sparse. Direct factorization is reliable for small systems, but fill-in during factorization and memory requirements become prohibitive for large three-dimensional meshes. Iterative methods are therefore central to industrial CFD.

Stationary iterations include Jacobi and Gauss–Seidel methods. In Jacobi iteration, each component is updated using values from the previous iteration. Gauss–Seidel uses newly updated values as soon as they are available. Successive over-relaxation introduces a factor selected to accelerate convergence in favorable cases. These methods are simple and useful as smoothers or building blocks, but they usually converge slowly for large, ill-conditioned systems.

Krylov-subspace methods construct successive approximations in spaces generated by residuals and matrix-vector products. Conjugate Gradient (CG) is appropriate for symmetric positive-definite matrices. Generalized Minimal Residual (GMRES) is applicable to nonsymmetric systems and minimizes the residual over an expanding Krylov space, but memory and orthogonalization costs grow with iteration count. Restarted GMRES limits memory at the possible expense of convergence. BiCGSTAB is another option for nonsymmetric systems and may require less storage, although its convergence can be less predictable.

The choice of solver must match matrix properties. Using CG for a strongly nonsymmetric convection-dominated system is mathematically inappropriate. Conversely, using a high-cost general nonsymmetric solver for a well-conditioned symmetric problem wastes resources. The pressure-correction system in incompressible flow, momentum equations, energy equations, and coupled block systems may each require different treatments.

Preconditioning: Making the Linear System Easier

A preconditioner transforms the system into one with more favorable spectral properties. If \(\mathbf{M}\) approximates \(\mathbf{A}\) and can be solved economically, a left-preconditioned system is

\[

\mathbf{M}^{-1}\mathbf{A}\mathbf{x}=\mathbf{M}^{-1}\mathbf{b}.

\]

The goal is not to compute \(\mathbf{A}^{-1}\) exactly, but to reduce the number of expensive Krylov iterations.

Jacobi or diagonal scaling is inexpensive and parallel-friendly, but often weak. Incomplete LU (ILU) factorization retains selected entries of an approximate factorization and can be effective for structured sparse systems. Its strength depends on the level of fill and drop tolerances, while parallel implementation may be more difficult. Algebraic multigrid (AMG) constructs coarse representations directly from matrix information and is particularly effective for elliptic pressure and diffusion operators. Domain-decomposition methods partition the mesh and solve local subproblems, making them attractive for parallel computing.

Preconditioning is a major performance decision rather than a minor solver option. A method that requires twice as many floating-point operations per iteration may still be preferable if it reduces the iteration count by an order of magnitude. The preconditioner must also be updated when the matrix changes substantially; reusing an outdated preconditioner can cause stagnation.

Multigrid and the Removal of Error at Multiple Scales

Iterative smoothers rapidly reduce oscillatory, short-wavelength error but are inefficient against smooth, long-wavelength error. Multigrid addresses both by transferring the problem between fine and coarse levels. A typical cycle consists of:

1. smoothing the fine-grid equations;

2. computing the residual or defect;

3. restricting the residual to a coarser level;

4. solving or approximately solving the coarse problem;

5. prolongating the correction to the fine level;

6. applying post-smoothing.

Geometric multigrid uses a known hierarchy of meshes. Algebraic multigrid constructs a hierarchy from the sparse matrix and is therefore suitable for unstructured grids. A V-cycle visits each level in a relatively economical pattern, whereas a W-cycle performs more coarse-level work and may be more robust for difficult systems.

Multigrid is especially powerful for Poisson-type pressure equations and diffusion-dominated operators. Strong convection, anisotropy, discontinuous coefficients, and poor mesh quality can reduce its effectiveness. A successful implementation depends on compatible restriction and prolongation operators, appropriate smoothers, and a coarse-grid problem that captures the important global modes.

Residuals, Convergence, and Stopping Criteria

For a linear system, the residual is

\[

\mathbf{r}^{(k)}=\mathbf{b}-\mathbf{A}\mathbf{x}^{(k)}.

\]

For nonlinear equations, it is the defect \(\mathbf{r}^{(k)}=-\mathbf{G}(\mathbf{U}^{(k)})\), subject to sign convention. A normalized residual may be reported as

\[

\hat r^{(k)}=\frac{\|\mathbf{r}^{(k)}\|}{\|\mathbf{r}^{(0)}\|},

\]

although the normalization must be stated because different choices can produce very different plotted values.

Residual reduction is necessary but not sufficient evidence of a trustworthy CFD result. A small algebraic residual may coexist with significant discretization error, an incorrect boundary condition, or an under-resolved physical feature. Conversely, a residual may stagnate at a nonzero level while integral quantities such as drag, mass flow, or heat transfer have become insensitive to further iterations.

Good monitoring combines several measures:

• normalized equation residuals;
• global mass and energy imbalance;
• changes in forces, moments, pressure drop, and heat-transfer rates;
• maximum changes in primary variables;
• boundedness and physical admissibility;
• sensitivity to time-step and mesh refinement.

For an unsteady calculation, one should distinguish inner iterations within a time step from changes between physical time levels. The inner system should be solved sufficiently accurately that its error does not dominate the temporal truncation error. Driving every time step to an extremely small residual may be wasteful when the time step itself is large or when the solution is intentionally first-order accurate.

Practical Method Selection

Method selection should begin with the physical question. If the objective is wave-resolved transient prediction, the time step must resolve the relevant frequencies and an explicit method may be attractive. If the objective is a slow thermal evolution coupled to fast diffusion, an implicit or semi-implicit method is usually more efficient. For steady solutions, pseudo-time stepping can provide a controlled route to the steady residual equations; local or adaptive pseudo-time steps can accelerate convergence, but the pseudo-time solution is not automatically a physical transient.

The following trade-offs are typical:

• Explicit methods have simple updates and strong parallel scalability but restrictive stability limits.
• Fully implicit methods permit larger steps and handle stiffness better but require nonlinear and linear solves.
• Crank–Nicolson preserves transient detail well but may produce oscillations without suitable damping.
• Backward Euler is robust and strongly damping, making it useful for startup and difficult transients.
• BDF2 offers second-order accuracy with good stiff-problem behavior but requires history and careful startup.
• Newton–Krylov methods can be powerful for strongly nonlinear systems, provided Jacobian actions and preconditioning are effective.
• Multigrid can provide near-optimal scaling for suitable elliptic problems but is sensitive to operator and implementation quality.

A practical workflow is to begin with a stable method and conservative time step, verify boundary conditions and conservation, monitor both residuals and engineering outputs, and then increase temporal order or step size only after the solution behaves consistently.

Illustrative Case Study: Transient Cooling in a Narrow Channel

Background

Consider an illustrative numerical study of forced convection through a narrow channel containing a heated wall. The fluid enters at a prescribed temperature and velocity, while the wall receives a step change in heat flux. The channel is represented by a nonuniform mesh refined near the wall, where the thermal boundary layer develops rapidly.

Situation and Problem

An explicit scheme initially appears attractive because each control volume can be updated independently. However, the near-wall cells are much smaller than the core cells. The diffusive stability restriction based on the smallest cell forces a time step far below that suggested by the bulk flow residence time. The simulation becomes dominated by the cost of advancing many tiny steps.

Analysis

The semi-discrete energy equation contains convective and diffusive contributions. The local diffusion number in the refined near-wall region is much larger than in the core for the same \(\Delta t\). Consequently, stability is controlled by the smallest cell, not by the average mesh size. The physical thermal response of the channel is relatively slow, so resolving the extremely fast numerical diffusion mode with an explicit method is inefficient.

Approach and Solution

The solver is changed to backward Euler during the initial rapid adjustment, followed by BDF2 after the thermal field becomes smoother. At every time step, the nonlinear energy equation is linearized and solved using restarted GMRES with an ILU preconditioner. A stronger AMG preconditioner is used when the mesh is refined further. The time step is increased gradually, but it is limited by a temporal-accuracy study rather than by stability alone. The calculation monitors wall heat flux, outlet temperature, global energy balance, and normalized residuals.

Outcome

The implicit method permits time steps substantially larger than the explicit diffusion limit while maintaining bounded temperatures and acceptable energy conservation. A comparison using successively halved time steps shows that the reported wall heat-transfer coefficient changes consistently with second-order behavior after the BDF2 phase begins. The startup phase remains first-order, but its influence decreases as the total simulated time increases.

Learning Points

The case illustrates that unconditional linear stability does not eliminate the need for time-step refinement. It also demonstrates the interaction between mesh design, stiffness, temporal discretization, and linear solver performance. A poor preconditioner can erase the advantage of an implicit method, while an unnecessarily tight inner residual tolerance can waste computation without improving the engineering result.

Common Failure Modes and Professional Practice

Several recurring mistakes deserve attention. Choosing a time step solely from a formal stability limit ignores accuracy and physical time scales. Reporting a low residual without checking conservation can conceal a defective implementation. Applying excessive under-relaxation may make a calculation appear stable while producing painfully slow progress. Reusing a preconditioner after the matrix has changed substantially can cause Krylov stagnation. Comparing two temporal schemes without controlling spatial and iterative errors does not establish temporal accuracy.

Poor mesh quality also affects algebraic convergence. Highly skewed or stretched cells can produce ill-conditioned matrices, while abrupt changes in coefficient values can challenge AMG coarsening. Boundary conditions must be included consistently in the residual and matrix; an algebraically converged solution with inconsistent boundary treatment is not a valid CFD solution.

Professional practice requires documenting the time-integration method, time-step strategy, CFL values, nonlinear and linear tolerances, preconditioner, residual normalization, and convergence monitors. Reproducible studies should distinguish physical time, pseudo-time, inner iteration count, and wall-clock cost. For safety-critical or design-critical computations, conservation checks and independent refinement studies are essential.

Technology and Emerging Directions

Modern CFD solvers increasingly combine matrix-free Krylov methods, block preconditioners, GPU-oriented sparse kernels, and algebraic multigrid. Matrix-free approaches avoid explicitly storing the full Jacobian and instead evaluate Jacobian-vector products, reducing memory use for large systems. Automatic differentiation can assist in generating consistent derivatives for Newton-type methods, although derivative cost and memory remain important considerations.

Adaptive time stepping is becoming more closely coupled with error estimation, local truncation-error indicators, and output-based accuracy requirements. In multiphysics simulations, partitioned and monolithic approaches are being compared according to stability, coupling strength, and parallel efficiency. Machine-learning techniques may assist parameter selection or predict preconditioner performance, but they do not remove the need for conservation, stability analysis, verification, and physical validation. Any automated solver decision must remain auditable and bounded by numerical safeguards.

Best-Practice Framework

A dependable temporal and algebraic solution strategy follows a sequence of checks:

1. Identify the fastest relevant physical and numerical scales.

2. Estimate convective and diffusive CFL restrictions.

3. Select explicit, implicit, or semi-implicit integration according to stiffness and accuracy needs.

4. Start with a robust scheme and verify boundary-condition implementation.

5. Match the linear solver to matrix symmetry, coupling, and conditioning.

6. Use preconditioning and multigrid where the operator supports them.

7. Set inner tolerances in relation to temporal and spatial errors.

8. Monitor residuals together with conservation and engineering quantities.

9. Perform time-step, mesh, and solver-tolerance sensitivity studies.

10. Record all numerical settings needed for reproducibility.

Chapter Summary

Temporal discretization transforms semi-discrete CFD equations into a sequence of algebraic problems. Explicit methods are inexpensive per step but subject to CFL restrictions, particularly when fine meshes create severe diffusive limits. Implicit methods permit larger steps and are better suited to stiff systems, but they require nonlinear linearization and the solution of sparse algebraic systems. Backward Euler is robust and dissipative; Crank–Nicolson is second-order but less damping; BDF2 provides second-order accuracy for many stiff applications.

Stability, accuracy, and algebraic convergence are distinct requirements. Sparse iterative methods such as Gauss–Seidel, Krylov solvers, preconditioned GMRES, and multigrid are central to large CFD calculations. Preconditioning improves the conditioning experienced by the iterative method, while multigrid removes error components at multiple spatial scales. Reliable CFD practice combines residual monitoring with conservation checks, engineering outputs, and time-step and mesh refinement. The best method is selected by balancing physical fidelity, stiffness, computational cost, robustness, and reproducibility.

Key Terms

Temporal discretization: Approximation of time derivatives using values at discrete time levels.

Semi-discrete equation: An equation discretized in space but continuous in time.

Explicit method: A time-integration method whose new state is computed from known previous states.

Implicit method: A method in which the new state appears inside the equations that must be solved.

CFL number: A nondimensional measure relating time-step size to information propagation across a cell.

Diffusion number: A nondimensional measure relating a diffusive time step to the square of a spatial scale.

Stiffness: The presence of rapidly decaying or evolving modes that impose severe step restrictions on explicit methods.

Backward Euler: A first-order implicit time-integration scheme using the residual at the new time level.

Crank–Nicolson method: A second-order method based on averaging old- and new-time residuals.

BDF2: A second-order backward differentiation formula using two previous solution levels.

Multistep method: A time-integration method using several previous solution states.

Residual: The imbalance remaining when an approximate solution is substituted into the discrete equations.

Jacobian: The matrix of derivatives of nonlinear algebraic equations with respect to the unknowns.

Krylov-subspace method: An iterative linear solver based on spaces generated by matrix-vector products.

Preconditioner: An approximate inverse or transformation used to accelerate iterative convergence.

GMRES: A Krylov method that minimizes the residual for general nonsymmetric systems.

Multigrid: A strategy that removes error components on a hierarchy of spatial or algebraic scales.

Under-relaxation: Updating an iterate by only a fraction of the computed correction.

Pseudo-time stepping: Artificial time advancement used to obtain a steady solution.

Temporal order: The power governing how temporal discretization error decreases as the time step is refined.

Review Questions

Short Answer Questions

1. What is the difference between a semi-discrete and a fully discrete CFD equation?

2. Why does an explicit diffusion scheme become expensive on a highly refined mesh?

3. Define stiffness in the context of CFD time integration.

4. What information does a linear-system residual provide?

5. Why is preconditioning used with Krylov-subspace methods?

Long Answer Questions

1. Derive the forward Euler and backward Euler forms for a semi-discrete conservation equation and compare their stability and computational requirements.

2. Explain the differences among backward Euler, Crank–Nicolson, and BDF2 with respect to order, damping, startup, and suitability for stiff flows.

3. Describe how a nonlinear implicit CFD equation is linearized and solved using Newton or approximate Newton methods.

4. Explain the operating principles of multigrid and discuss why it is effective for pressure and diffusion equations.

5. Develop a systematic strategy for selecting time steps, linear solvers, preconditioners, and convergence criteria for an industrial transient simulation.

Analytical / Application Questions

1. A mesh contains near-wall cells one hundred times smaller than the core-cell width. Explain how this affects explicit convective and diffusive time-step restrictions and recommend a suitable strategy.

2. A GMRES calculation shows steadily decreasing residuals, but the predicted drag changes significantly between successive time steps and the global mass imbalance remains large. Diagnose the situation and propose corrective actions.

3. Compare two simulations of a thermal transient: one uses Crank–Nicolson with a large time step and the other uses backward Euler with a smaller time step. Explain which results should be examined before deciding that the first simulation is more accurate.

22. Multiple Choice Questions with Answers

Question 1

Which statement best distinguishes an explicit time-integration method from an implicit method?

A. An explicit method always has second-order accuracy.

B. An explicit method computes the new state using known previous states without solving a global new-time system.

C. An implicit method cannot be used for nonlinear equations.

D. An explicit method is unconditionally stable for diffusion.

Correct Answer: B. An explicit method computes the new state using known previous states without solving a global new-time system.

Question 2

For a fine mesh in a viscous flow, the explicit diffusive time-step restriction is primarily controlled by which quantity?

A. The largest cell volume

B. The mean flow residence time only

C. The smallest spatial length scale, generally through its square

D. The number of output files

Correct Answer: C. The smallest spatial length scale, generally through its square

Question 3

Which method is second-order accurate in time and uses two previous solution levels?

A. Forward Euler

B. Backward Euler

C. BDF2

D. Jacobi iteration

Correct Answer: C. BDF2

Question 4

A stable implicit calculation produces strongly damped temperature oscillations when the time step is increased. What is the most appropriate interpretation?

A. Stability guarantees temporal accuracy.

B. The calculation may be stable but temporally inaccurate because of numerical damping.

C. The residual must be exactly zero.

D. The mesh has automatically become second-order accurate.

Correct Answer: B. The calculation may be stable but temporally inaccurate because of numerical damping.

Question 5

Which Krylov solver is generally appropriate for a symmetric positive-definite matrix?

A. Conjugate Gradient

B. GMRES only

C. BiCGSTAB only

D. Newton's method

Correct Answer: A. Conjugate Gradient

Question 6

What is the principal purpose of a preconditioner?

A. To change the physical boundary conditions

B. To reduce the number of unknowns by deleting equations

C. To transform the linear system so that iterative convergence is improved

D. To increase the temporal truncation error

Correct Answer: C. To transform the linear system so that iterative convergence is improved

Question 7

Why does multigrid use coarse levels?

A. To remove smooth, long-wavelength error that fine-grid smoothing handles poorly

B. To eliminate the need for spatial discretization

C. To guarantee exact nonlinear convergence in one iteration

D. To replace all boundary conditions with approximate ones

Correct Answer: A. To remove smooth, long-wavelength error that fine-grid smoothing handles poorly

Question 8

Which observation provides the strongest evidence that a transient CFD result is not yet trustworthy, despite a decreasing algebraic residual?

A. The output file is large.

B. The engineering quantity of interest and conservation imbalance continue to change substantially.

C. The solver uses a sparse matrix.

D. The method is implicit.

Correct Answer: B. The engineering quantity of interest and conservation imbalance continue to change substantially.

Question 9

In Newton linearization of a nonlinear algebraic system, what does the Jacobian represent?

A. The matrix of derivatives of the equations with respect to the unknowns

B. The list of time-step values

C. The mesh connectivity alone

D. The normalized residual history

Correct Answer: A. The matrix of derivatives of the equations with respect to the unknowns

Question 10

A simulation contains fast chemical or diffusive modes but seeks a slow thermal response. Which strategy is generally most appropriate?

A. Use an explicit method with a time step based only on the slow response.

B. Use an implicit or semi-implicit method and verify accuracy through time-step refinement.

C. Ignore the fast modes because they do not affect stability.

D. Use no algebraic convergence criterion.

Correct Answer: B. Use an implicit or semi-implicit method and verify accuracy through time-step refinement.