Discretization Schemes – Upwinding


In what follows a problem which arises as a consequence of specific properties of solutions to difference equations shall be presented. In contrary to what is often thought, the phenomena which manifests itself as non-physical oscillation of the numerical solution from one grid point to the next does not arise simply from nonlinearities of Navier–Stokes equations (NSE).

“All About CFD…” – Index, Resources and Free Online Courses

The Cell-Péclet/Reynolds Problem

Discretized linear equations can exhibit the same kind of behavior as is seen in the NS equation solutions and to appreciate this it shall be presented (as it is usually presented) for the following 1-D steady-state convection-diffusion problem, a simplification homogeneous Burgers’ equation, one of the few nonlinear equations to possess an exact solution:

Burgers Equation

where U is a known constant.


Discretize the simplified equation with centered differences to obtain:


where h is the discretization step size. Intentionally multiplying by h/U leaves this in the form. Follow the following caption:


The starting point toward understanding the cell- Péclet problem stems from a property of elliptic and parabolic differential equations known as the maximum principle.

Starting from the basic fact from calculus that if a function f(x) satisfies f” > 0 on an interval [a, b], then it can only achieve its maximum on the boundary of that interval. For partial differential equations, the same idea allows to draw very useful conclusions from the properties of the solutions and the domain of a given problem and it shall allow a quick view on the properties of linear discretization of such entities.

We shall start by defining L as an elliptic linear partial differential operator (as the Laplacian operator for example). Formally for an elliptic operator with a bounded domain we shall write:

elliptic operator

and the boundary condition satisfies:

elliptic operator

Then the maximum principle states that the maximum and minimum values of u occur on the boundary.

Since it’s difference approximations are our concern, it shall be stated that a similar principle exists for difference approximations to elliptic equations of the following general form:

elliptic operator

To guarantee  such a maximum principle for discretizations two conditions must be satisfied:

  • Non-negativity:
  • diagonal dominance:

These requirements are not dependent upon the spatial dimension of the problem or on the specific formulation for differencing involved in the construction of the approximation. In particular, they are sufficient even on unstructured grids.
We shall observe the discretized version of the maximum principle in a local sense rather than the whole domain. Specifically, if we consider placing a boundary  around each grid, the maximum principle implies that the grid function value at the point (i, j) cannot exceed that of any of its nearest neighbors, and because there is a corresponding minimum principle, the grid-function value cannot be lower than that of any neighbor, meaning it cannot oscillate from point to point.


Remedies for the Cell-Péclet/Reynolds Problem

There are some several widely-used ways to treat the cell-Péclet/Reynolds problem. In what follows some popular approximations such as first and second order upwinding and QUICK shall be discussed.  Nearly all notable commercial codes provide means to employ these specific approximations including the standard centered difference approximation.

  • First-Order Upwind:
    This is probably the first approach to treat the cell-Péclet problem, originating from S. V. Patankar (father to the original SIMPLE algorithm).
    Patankar had taken an ad-hoc “physical reasoning” arguing that as it is clear from the representation of centered difference approximation to the velocity gradient both upstream and downstream flow behavior are sensed. Hence a least biased difference approximation, using only information carried in the flow direction in its simplest form is offered. After employing the simplification as in the above procedure we arrive to the following form:
    difference equation_2
    This difference equation satisfies non-negativity and hence the extermum conditions for difference equations.
    On the other hand, the specific approximation’s shortcoming may be viewed by performing a simple procedure.
    We shall take the non-dimensional for of the linearized steady-state Burger’s equation:
    difference equation_2
    and remember that upwind difference approximation we employ the following backward difference for the velocity gradient at grid i:
    difference equation_2expansion of u at grid point i-1 around grid i yields:
    difference equation_2
    which may then be replaced to the linearized steady-state Burger’s equation above to achieve the following form:
    difference equation_2
    we shall write the above equation as:
    difference equation_2
    this equation is termed the modified equation and it is actually it rather than the original partial differential equation which is solved to second order accuracy instead of the first order accuracy the discretization implies.
    In the above equation, Re* is an effective Péclet/Reynolds defined as:
    difference equation_2
    What is clear from the above analysis is that first-order upwind approximation effectively increases the coefficient of the diffusion (!), a phenomena often called numerical diffusion.
  • Second-Order Upwind:

    Here we follow (without the elaboration as taken for the first-order upwind) a reference from W. Shyy (1985) assisted by other pointed out references.
    The form of the linearized steady-state Burger’s equation after applying second-order one-sided differences:
    difference equation_2
    We note the switching applied in order to keep the difference into the flow direction.We observe and conclude the following important notions:


    • As is obvious, the discretization is second-order accurate. This means a five-point discrete representation, resulting in a more complicated matrix structure if implicit time stepping is employed in the solution algorithm.
    • As it turns out the leading terms of the modified equation are those of the original equation which means that the physics is not compromised by this approximation.
    • Because the second order upwind differencing connects two upstream nodes, it is necessary to follow special practices for the finite difference nodes adjacent to a boundary; otherwise reference will have to be made to a value outside the flow domain. There are two popular techniques to deal with that notion. M. Atias et al. (“Efficiency of Navier-Stokes Solver” – 1977) proposed switching back to first order upwind differencing for cells adjacent to a solid boundary.
      A different route proposed by W. Shyy and S. M. Correa (“A Systematic Comparison of Several Numerical Schemes for Complex Flow Calculations” – 1985) is to use the hybrid formulation near the boundaries. The hybrid formulation at boundaries is the same as the upwind formulation except that when the boundary cell -Péclet/Reynolds number is less than two, central differencing is used.
    • the roots of the corresponding difference equation are all non-negative, independent of the value of cell Re. Hence, formally, there is no cell-Péclet/Reynolds restriction associated with this approximation.
    • We further note that referring to the theoretical difference equation analysis above neither non-negativity of the difference equation coefficients nor diagonal dominance hold for this scheme, implying that it does not satisfy a maximum principle in the sense described earlier. This reflects the fact that those conditions pose satisfactory but not necessary conditions meaning that in order to not draw erroneous conclusions about the cell-Péclet problem more emphasis should be placed on finding the solutions to the difference equation.
    • Finally, as guidance for whether to make use of a second rather than a first order approximation we shall add that When the flow is aligned with the flow gradients (e.g., laminar flow in a rectangular duct modeled with a quadrilateral or hexahedral grid) the first-order upwind discretization may be acceptable. When the flow is not aligned with the flow gradients (i.e., when it crosses the grid lines obliquely), however, first-order convective discretization increases the numerical discretization error (numerical diffusion).This is quite an important topic of which there is a consensus among CFD practitioners, i.e. the impact of alignment of flow gradients on phenomena like numerical diffusion an dispersion.
      Not to delve to deep into the topic, looking at the pure convection equation:


      The above PDE’s modified equation (the PDE which the exact solution to the discretized equations satisfies) is:

      Meaning that the numerical behavior of the discretization scheme largely depends on the relative importance of dispersive and dissipative effects, specifically, a low order scheme along with a mesh which is not aligned with the flow’s gradients will tend to smear all over the place:

      Quad Vs. tri mesh: contours of axial velocity for the case of inviscid co-flow jet

      There is of course no benefit in this aspect when there is no dominant flow gradient direction:

      Quad Vs. tri mesh: contours of temperature for inviscid flow


      Poly mesh would not be as aligned with the flow gradients for very simple flows such as the above, but due to the multiplicity of faces it has 6 optimal direction (if it has 12 faces) in contrast with hex which have 3 optimal faces, it may be more accurate than hex mesh for recirculating flows even for those which seem optimal for hex as they have a cartesian domain such as the qubic lid-driven cavity:

      Referring to mesh methodologies we then add that for triangular and tetrahedral grids, since the flow is never aligned with the grid, you will generally obtain more accurate results by using the second-order discretization. For quad/hex grids, you will also obtain better results using the second-order discretization, especially for complex flows.
      For a simple flow that is aligned with the flow gradients (e.g., laminar flow in a rectangular duct modeled with a quadrilateral or hexahedral grid), the numerical diffusion will be naturally low, so you can generally use the first-order scheme instead of the second-order scheme without any significant loss of accuracy.
      For most cases, you will be able to use the second-order scheme from the start of the calculation. In some cases, however, you may need to start with the first-order scheme and then switch to the second-order scheme after a few iterations. For example, if you are running a high-Mach-number flow calculation that has an initial solution much different than the expected final solution, you will usually need to perform a few iterations with the first-order scheme and then turn on the second-order scheme and continue the calculation to convergence.
      In summary, while the first-order discretization generally yields better convergence than the second-order scheme, it generally will yield less accurate results, especially on tri/tet grids, concluding that strong emphasis should be taken to control convergence. STAR-CCM+ allows for various of direct ways to do just that by Modifying Algebraic Multigrid Parameters, Setting FAS Multigrid Parameters and Modifying Multi-Stage Time-Stepping Parameters for example.

  • QUICK:

    The quadratic upstream interpolation for convective kinematics scheme of B. P. Leonard (“A stable and accurate convection modelling procedure based on quadratic upstream interpolation” – 1979) is one of the most popular techniques for remedying the cell-Péclet/Reynolds problem.
    The form of the linearized steady-state Burger’s equation after applying second-order one-sided differences for the QUICK approximation is as follows:
    difference equation_2
    For the one-dimensional domain shown in the figure the Φ value at a control volume face is approximated using three-point quadratic function passing through the surrounding nodes and one other node on upstream side:
    We observe and conclude the following important notions:


    • Although the original intention by B. P. Leonard to produce a method with third order accuracy this was actually with respect to a location in the grid stencil that did not correspond to any grid point. It is most important to note that in contrast to what is frequently thought QUICK is second-order accurate.
    •  We note the switching applied in order to keep the difference into the flow direction only in QUICK it is referred to as upwind biased as the stencil extends to both downstream and upstream of a specific grid point.
    • As for the second-order upwind approximation neither non-negativity of the difference equation coefficients nor diagonal dominance hold for this scheme, implying that it does not satisfy a maximum principle in the sense described earlier but in contrast to the former it has a cell-Péclet/Reynolds restriction:
      difference equation_2
    • Although the above restriction there is a general favorable experience regarding the cell-Péclet/Reynolds problem.  Nevertheless, although being  a very accurate scheme, in regions with strong gradients, overshoots and undershoots can occur  lead to stability problems in the

Practical guidelines

  • When the flow is aligned with the flow gradients (for example, laminar flow in a rectangular duct modeled with a quadrilateral or hexahedral mesh) the first-order upwind discretization may be acceptable, otherwise, (as is the case when the flow crosses the mesh lines obliquely),  first-order discretization increases the numerical discretization error through numerical diffusion. For triangular
    and tetrahedral meshes, since the flow is never aligned with the flow gradients, more accurate results shall be obtained by using the second-order discretization. Keep in mind (especially for complex flows) that even for hex meshes, you will also obtain better results using the second-order discretization.
  • For most cases, it is possible to use the second-order scheme from the start of the calculation. There  are some cases, however, where it is preferable to perform an initial few iterations with the first-order scheme and then switch to the second order (high-Mach-number flow for example).
  • The QUICK scheme may provide better accuracy than the secondorder scheme for rotating or swirling flows. The QUICK scheme is applicable to hexahedral
    meshes. An alternative scheme, the 3rd order MUSCL scheme, may be used on all types of meshes.
    Keep in mind howeever that the second-order upwind scheme is sufficient in most cases and the QUICK scheme will not provide significant improvements in accuracy.

We shall close the discussion about the cell- Reynolds problem with a summary table regarding the approximations presented and reference to their cell-Reynolds restriction:

difference equation_2

Back to “All About CFD…” Index

MSI Engineering Software – Mechanical Engineering Software Distributor (STAR-CCM+, FloEFD, FloTherm, MSC Nastran, etc…)
TAZ-Engineering: Thermal Management and CFD Consultancy

Leave a Reply