. .ri.‘
9:694.
V4411

134'. a

w! .i
.l 311
1 zllhl‘

 

 

 

p"
‘-r7_:'I
'C {{2/07

This is to certify that the
dissertation entitled

DISCRETE-ERROR-TRANSPORT EQUATION
FOR ERROR ESTIMATION IN CFD

presented by

YUEHUI QIN

has been accepted towards fuifillment
of the requirements for the

Ph.D degree in Mechanical Enmeering

 

 

W

Major Profefiéor’s Signature

W K/ 200%
V

Date

MSU is an Affirmative Action/Equal Opportunity Institution

 

LIBRARY '
Michigan State
University

 

 

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

6/01 cJCIRC/DateDue.pss-p.15

DISCRETE-ERROR-TRANSPORT EQUATION
FOR ERROR ESTIMATION IN CFD

By

Yuehui Qin

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY
Department of Mechanical Engineering

2004

ABSTRACT

DISCRETE-ERROR-TRANSPORT EQUATION
FOR ERROR ESTIMATION IN CFD

By

Yuehui Qin

With computational fluid dynamics (CFD) becoming more accepted and more widely
used in industry for design and analysis, there is increasing demand for not just more
accurate solutions, but also error bounds on the solutions. One major source of error is
from the grid or mesh. A number of methods have been developed to quantify errors in
solutions of partial differential equations (PDEs) that arise from poor-quality or
insufficiently fine grids/meshes. For PDEs of interest to CF D, it has been shown that the
error at one location could be generated elsewhere and then transported there, and thus is
not a function of the local mesh quality and the local solution. So, a transport equation
for error is needed to understand the generation and evolution of errors. Error transport
equations have been developed for finite-element methods but not for finite-difference
(FD) and finite-volume (FV) methods.

In this study, a method is developed for deriving error-transport equations for
estimating grid-induced errors in solutions obtained by using FD and FV methods. The
error-transport equations derived are discrete in that they depend only on the FD or FV
equations and are independent of the PDEs that the FD or FV equations are intened to
represent. The usefiilness of the DETEs developed was evaluated through test problems
based on four one-dimensional (l-D) and two two-dimensional (2-D) PDES. The four 1-

D PDEs are the advection-diffiision equation, the wave equation, the inviscid Burgers

equation, and the steady Burgers equation. The two 2-D PDEs are the 2-D advection-
diffusion equation and the system of Euler equations. For PDEs that are not linear,
linearization procedures were proposed and examined. For all test problems based on 1-
D PDEs, the residual is modeled by the leading term of the remainder in the modified
equation for the FD or FV equation. The residual was also modeled by using functional
relationship suggested by data mining, where actual residuals generated by the numerical
experiments were fitted by using least-square minimization. For all test problems, grid-
independent solutions were generated to assess how well the residuals are modeled and
how well grid-induced errors are predicted by the DETEs.

Results obtained show that if the actual residuals are used, then the DETEs can
predict the grid-induced errors perfectly. This is true for all test problems evaluated,
including those based on PDEs that are nonlinear and have time derivatives and for test
problems with weak solutions. Results obtained also show that the leading terms of the
modified equation is useful in modeling the residual if the grid spacing or cell size is
sufficiently small so that the leading terms are bounded, a condition that is often not
satisfied in practice. The usefulness of data mining in constructing residuals show the
power-law to produce better fit than local linear least square of smoothness, resolution,
aspect ratio and solution gradient. However, a more extensive database is needed before
this approach can be expected to yield a more generally applicable models for the
residual. The usefulness of Euler DETE in predicting grid-induced errors in the Navier-
Stokes solutions was also examined. Results obtained show that error predicted by Euler
DETE matches very well with the actual error for the hi gh-Reynolds-number Navier-

Stokes solutions.

ACKNOWLEDGEMENTS

I would like to express my deepest gratitude to my advisor Dr. Tom I-P. Shih for
his support and encouragement throughout my graduate studies, for his help in my

professional and personal growth, and for his insightful advices on the academic research.

I would also like to thank my co-advisor, Dr. Farhad J aberi, my committee
members Dr. Matt Mutka and Dr. Harold Schock, for their many valuable comments and
discussions on this dissertation. My thanks also go to Dr. 2.]. Wang for his help on the
numerical algorithms for the Euler equations on arbitrary-shaped elements. I appreciate

the help and friendship of my fellow graduate students in the CFD laboratory.

I am indebted to my father Xudian Qin, my mother Tairong Liu, my brother
Guodong Qin and my sister Yanhui Qin, for their continuous support during all these
years of studies. I am especially gratefiil to my husband Yongxiang (Daniel) Li for his
understanding and love during the past few years. His support and encouragement was in

the end what made this dissertation possible.

iv

TABLE OF CONTENTS

LIST OF TABLES ........................................................................................................... viii
LIST OF FIGURES ........................................................................................................... ix
NOMENCLATURE ........................................................................................................ xiii
CHAPTER 1 ....................................................................................................................... 1
INTRODUCTION .............................................................................................................. 1
1.1 Background 1
1.2 Literature Survey ................................................................................................ 4
1.2.1 Multiple-Grids Error Estimation ................................................................. 4
1.2.2 Single-Grid Algebraic Error Estimation ..................................................... 7
1.2.3 Single-Grid PDE Error Estimation ........................................................... 14
1.3 Motivations and Objectives .............................................................................. 18
CHAPTER 2 ..................................................................................................................... 20
DEVELOPMENT OF DISCRETE-ERROR-TRANSPORT EQUATION ...................... 20
2.1 Error Transport Concept ................................................................................... 20
2.1.1 Problem Description ................................................................................. 20
2.1.2 Grid-Independent Solution Generation ..................................................... 21
2.1.3 Imperfect Grid and Error Transport .......................................................... 22
2.2 Discrete-Error-Transport Equation (DETE) ..................................................... 24
2.2.1 Description of Model Equation ................................................................. 24
2.2.2 Numerical Method of Solution ................................................................. 24
2.2.3 Discrete-Error-Transport Equation (DETE) ............................................. 26
CHAPTER 3 ..................................................................................................................... 37
EVALUATION OF DISCRETE-ERROR-TRANSPORT EQUATION ON MODEL
EQUATIONS .................................................................................................................... 37
3.1 Derivation of Discrete-Error—Transport Equations ................................................ 38
3.1.1 Advection-Diffusion Equation (Linear and Steady) ................................ 38
3.1.2 Wave Equation (Linear and Unsteady) ..................................................... 40
3.1.3 Inviscid Burgers Equation (Nonlinear and Unsteady) .............................. 42
3.1.4 Steady Burgers Equation (Nonlinear and Steady) .................................... 45
3.2 Modeling the Residual in the Discrete-Error-Transport Equation .................... 49
3.2.1 Advection-Diffusion Equation .................................................................. 50
3.2.2 Wave Equation .......................................................................................... 51
3.2.3 Inviscid Burgers Equation ......................................................................... 51
3.2.4 Steady Burgers Equation ........................................................................... 52
3.3. Results and Discussion ..................................................................................... 53
3.3.1 Advection-Diffusion Equation .................................................................. 54
3.3.2 Wave Equation .......................................................................................... 5 8

3.3.3 Inviscid-Burgers Equation ........................................................................ 59

3.3.4 Steady-Burgers Equation .......................................................................... 62
3.4. Summary of Chapter 3 ...................................................................................... 64
CHAPTER 4 ..................................................................................................................... 83
ANALYSIS AND MODELING OF RESIDUAL ............................................................ 83
IN THE DISCRETE-ERROR-TRAN SPORT EQUATION ............................................ 83
4.1 Test Problems and Derivation of Discrete-Error—Transport Equations ............ 83
4.2. Behavior of the Residual in the Discrete-Error-Transport Equation ................ 88
4.3. Results of the Discrete-Error-Transport Equation in Estimating Errors ........... 93
4.4. Modeling the Residual in the Discrete-Error—Transport Equation .................... 94
CHAPTER 5 ................................................................................................................... 109
ERROR ESTIMATION FOR PROBLEMS GOVERNED BY EULER EQUATIONS 109
5.1 Description of Governing Equations .............................................................. 109
5.2 Finite-Volume Method of Solution ................................................................. 111
5.3 Discrete-Error—Transport Equations (DETEs) ................................................ 115
5.4 Test Problems .................................................................................................. 117
5.4.1 Test Problem 1: Oblique Shock Problem ................................................ 117
5.4.2 Test Problem 2: Inviscid Flow around Airfoil ........................................ 118
5.5 Treatment of Boundary Conditions and Convergence Criteria ...................... 119
5.5.1 Zero flux/flow boundary conditions ....................................................... 119
5.5.2 Free stream boundary conditions ............................................................ 119
5.5.3 Convergence Criteria .............................................................................. 121
5.6 Results ............................................................................................................. 123
5.6.1 Oblique Shock Problem. ......................................................................... 123
5.6.2 Airfoil Problem. ...................................................................................... 124
5.7 Summary of Chapter 5 .................................................................................... 127
CHAPTER 6 ................................................................................................................... 138
ERROR ESTIMATION FOR PROBLEMS GOVERNED BY NAVIER-STOKES
EQUATIONS .................................................................................................................. 138
6.1 Description of 2-D Navier-Stokes Equations ................................................. 138
6.2 Euler DETE for Error Estimation of Navier-Stokes Solutions ....................... 140
6.3.1 Clean Airfoil Problem. ............................................................................ 142
6.3.2 Iced Airfoil Problem. .............................................................................. 143
6.4 Results ............................................................................................................. 144
6.4.1 Clean Airfoil Problem. ............................................................................ 144
6.4.2 Iced Airfoil Problem. .............................................................................. 145
6.5 Summary of Chapter 6 .................................................................................... 147
CHAPTER 7 ................................................................................................................... 165
CONCLUSIONS ............................................................................................................. 165
APPENDIX A ................................................................................................................. 170
Publications Associated with this Work ......................................................................... 170

vi

BIBILOGRAPHY ........................................................................................................... 1 72

vii

LIST OF TABLES

Table 2. 1 Uniform grids used to get grid-independent solution for the 2-D chamber

problem ..................................................................................................................... 31
Table 3.1 Summary of Cases Studied for 1-D Advection—Diffusion Equation ............... 65
Table 4.1 Summary of Parameters Investigated for the 1-D Test Problem ..................... 98

Table 4.2 Summary of Parameters Investigated in Sequence 1 for the 2-D Test Problem
................................................................................................................................... 99

Table 4.3 Summary of Parameters Investigated in Sequence 2 for the 2-D Test Problem
................................................................................................................................. 100

Table 4.4 Summary of Parameters Investigated in Sequence 3 for the 2-D Test Problem
................................................................................................................................. 101

Table 4.5 Summary of Parameters Investigated in Sequence 4 for the 2-D Test Problem
................................................................................................................................. 102

Table 4.6 Summary of Parameters Investigated in Sequence 5 for the 2-D Test Problem
................................................................................................................................. 102

Table 4.7 Summary of Parameters Investigated in Sequence 6 for the 2-D Test Problem
................................................................................................................................. 104

viii

LIST OF FIGURES

Figure 2.1 Geometry and flow conditions of the 2-D chamber problem ......................... 31
Figure 2.2 Grid-independent solution for the 2-D chamber problem .............................. 32
Figure 2.3 Perturbed grids at location (h,h) for the 2-D chamber problem ..................... 33
Figure 2.4 Error distribution in the flow domain for the 2-D chamber problem with
perturbed grid ............................................................................................................ 34
Figure 2.5 Interpolating from coarse to fine mesh ........................................................... 35
Figure 2.6 Residual (R8) computed for Eq. (2.14). .......................................................... 35
Figure 2.7 Residual (RC) computed for Eq. (2.17) ........................................................... 36
Figure 3.1 Actual and modeled residual for advection-ditfusion equation: .................... 67
Figure 3.2 Actual and modeled residual for advection-diffusion equation: .................... 68

Figure 3.3 Actual and predicted solution error (e,) for advection-diffusion equation: 70

Figure 3.4 Actual versus modeled residual for advection-diffusion equation: S = 1. ..... 71

Figure 3.5 Actual versus predicted error for advection-diffusion—equation: S = 1. .......... 71
Figure 3.6 Error = F(gn'd spacing) for advection-diffusion equation: ............................. 72
Figure 3.7 Error = F(smoothness) for advection-diffusion equation: .............................. 74
Figure 3.8 Actual versus modeled residual for wave equation. ....................................... 75
Figure 3.9 Actual versus predicted error for wave equation ............................................ 75
Figure 3.10 Results for wave equation: I=201, Co=0.8 at n = 100 and n = 200 .............. 76
Figure 3.11 Actual versus modeled residual for inviscid Burgers equation: ................... 77
Figure 3.12 Actual versus predicted error for inviscid Burgers equation: ........................ 78
Figure 3.13 Actual versus modeled residual for inviscid Burgers equation: .................... 79
Figure 3.14 Actual versus predicted error for inviscid Burgers equation: ........................ 80

ix

Figure 3.15 Actual versus modeled residual for steady Burgers equation: ...................... 81

Figure 3.16 Actual versus predicted error for steady Burgers equation: .......................... 82
Figure 4.1 Schematic of removing two grid points after x = xr for the 1-D test problem.
................................................................................................................................... 98
Figure 4.2 Residual computed by using Eq. (4.13) for the 1-D test problem. ................. 98
Figure 4.3 Residual computed by using Eq. (4.14) for the 2-D test problem with
Sequence 1 (Table 4.2). ............................................................................................ 99
Figure 4.4 Schematic showing where grid lines are removed to study the effects of S > 1
and S < l for the 2-D test problem .......................................................................... 100
Figure 4.5 Residual computed by using Eq. (4.14) for the 2-D test problem with
Sequence 2 (Table 4.3): .......................................................................................... 100
Figure 4.6 Residual computed by using Eq. (4.14) for the 2—D test problem with
Sequence 3 (Table 4.4): .......................................................................................... 101
Figure 4.7 Residual computed by using Eq. (4.14) for the 2-D test problem with
Sequence 4 (Table 4.5): .......................................................................................... 102
Figure 4.8 Residual computed by using Eq. (4. 14) for the 2-D test problem with
Sequence 5 (Table 4.6): .......................................................................................... 103
Figure 4.9 Residual computed by using Eq. (4.14) for the 2-D test problem with
Sequence 6 (Table 4.7). .......................................................................................... 104
Figure 4.10 Actual absolute error calculated using Eq. (4.12) for the 2-D test problem:
................................................................................................................................. 105
Figure 4.11 Predicted absolute error using Eq. (4.11), with actual residual calculated
using Eq. (4.10) for the 2-D test problem: .............................................................. 105
Figure 4.12 Actual residual (RA) and residual modeled (RM) by Eq. (4.18) for the 1-D
test problem ............................................................................................................. 106
Figure 4.13 Actual residual (RA) and residual modeled (RM) by Eq. (4.19) for the 1-D
test problem ............................................................................................................. 106
Figure 4.14 Actual residual (RA) and residual modeled (RM) by Eq. (4.20) for the l-D
test problem ............................................................................................................. 107
Figure 4.15 Definition of 0t and B for Eq. (4.21) for the 2-D test problem. ................ 107

Figure 4.16 Actual residual (RA) and residual modeled (RM) by Eq. (4.25) for the 2-D

test problem ............................................................................................................. 108
Figure 4.17 Error in predicted residual: (RA-RM)/RA for the 2-D test problem. ......... 108
Figure 5.1 Schematic of Oblique Shock Problem. ......................................................... 128
Figure 5.2 Computational domain and local mesh for the airfoil problem. ................... 129
Figure 5.3 Schematic of fi’eestream boundary conditions ............................................. 129
Figure 5.4 Grid-independent solution in density for the oblique shock problem: ......... 130
Figure 5.5 Exact Residual Rs in the continuity equation for the oblique shock problem:

321 x161case .......................................................................................................... 130
Figure 5.6 Exact Residual R1; in the continuity equation for the oblique shock problem:

161 x 81 case ........................................................................................................... 131
Figure 5.7 Exact Residual R8 in the continuity equation for the oblique shock problem:

81 x 41case .............................................................................................................. 131
Figure 5.8 Exact Residual R8 in the continuity equation for the oblique shock problem:

41 x 21 case ............................................................................................................. 132
Figure 5.9 Actual absolute error in solution of density for the oblique shock problem:

321 x 161 case ......................................................................................................... 132
Figure 5.10 Actual absolute error in solution of density for the oblique shock problem:

161 x 81 case ........................................................................................................... 133
Figure 5.11 Actual absolute error in solution of density for the oblique shock problem:

81 x 41 case ............................................................................................................. 133
Figure 5.12 Actual absolute error in solution of density for the oblique shock problem:

41 x 21 case ............................................................................................................. 134
Figure 5.13 Grid-independent solution in density for the inviscid airfoil problem: ...... 134

Figure 5.14 Exact Residual Rs in the continuity equation for the inviscid airfoil problem:
385 x 129 case ......................................................................................................... 135

Figure 5.15 Exact Residual Rg in the continuity equation for the inviscid airfoil problem:
193 x 65 case ........................................................................................................... 135

xi

Figure 5.16 Exact Residual R3 in the continuity equation for the inviscid airfoil problem:
97 x 33 case ............................................................................................................. 136

Figure 5.17 Actual absolute error in solution of density for the inviscid airfoil problem:
385 x 129 case ......................................................................................................... 136

Figure 5.18 Actual absolute error in solution of density for the inviscid airfoil problem:
193 x 65 case ........................................................................................................... 137

Figure 5.19 Actual absolute error in solution of density for the inviscid airfoil problem:
97 x 33 case ............................................................................................................. 137

IMAGES IN THIS DISSERTATION ARE PRESENTED IN COLOR.

xii

O

DX

e,E

LI)

LDN

<<1<C'-1

:3‘

NOMENCLATURE

aspect ratio

advection speed
diameter of chamber inlet or exit duct
grid spacing

solution error

total energy
differential operator
difference operator
difference operator
pressure

residual

resolution factor
resolution factor
smoothness factor
Reynolds number
temperature
x-component velocity

y-component velocity

velocity vector
velocity magnitude

grid spacing in x-direction

xiii

Ay

At

Greek

grid spacing in y-direction
grid spacing in x-direction
grid spacing in y-direction

time-step size

dynamic viscosity
kinematic viscosity
density

ratio of specific heats

Subscripts and Superscripts

ab

GI, gi

1°C

absolute error
coarse-grid
grid-independent solution
cell index

modeled

time step

relative error

xiv

CHAPTER 1
INTRODUCTION

1.1 Background

Computational fluid dynamics, commonly referred to as CFD, is concerned with
obtaining numerical solution to partial differential equations that govern fluid flow
problems by using the computer. The equations governing fluid flow problems are the
continuity, momentum (Navier-Stokes), and energy equations as well as additional

equations such as species concentrations if there is combustion.

CFD began to develop in the mid 1960's and its first successes came to prominence in
the 1970's. The CFD-service industry started in the 1980's and expanded significantly in
the 1990's, with the development of high-speed low-cost computers and the
improvements in numerical algorithms. There has been considerable growth in the
development and application of CFD in all aspects. In design and development, CFD
programs are now considered to be standard numerical tools, widely utilized within all
kinds of industries from aerospace, automotive, materials processing, electronics

appliances to bioengineering and biochemistry.

With CFD becoming more accepted and more widely used in industry for design and
analysis, there is increasing demand for not just more accurate solutions, but also error

bounds on the solutions. The accuracy of a CFD solution is determined by the difference

between the flow problem being studied and the numerical solutions. The error in each
CFD solution can be attributed to three sources (see, e.g., Roache (1998), Oberkarnpf et
al.(1995), and Cosner et al.(2004)). The first is inadequate modeling of physics that are
not resolved by first principles. For any thermal-fluids phenomenon, the physical model
should account for all the main physics features relevant to the problem of interest.
Roache (1998) summarized the possible error sources in the model as the physical
phenomenon not completely understood and the uncertainties of parameters used in the
model. In order for CFD to simulate a problem and minimize the model errors, the

physics and the use of physical models need to be better understood.

The second is non-physical effects such as numerical diffirsion and dispersion that
result when the governing partial differential equations (PDEs) are discretized into
algebraic equations on a discrete domain by finite—difference, finite-volume, or finite-
element methods. Because the governing equations are usually complex and nonlinear,
errors are easily created when using the numerical methods to approximate these

equations.

The third source of error is from a poor quality or an insufficiently fine grid or mesh.
The grid or mesh is the discrete spatial domain used to generate the approximate
numerical solution. It must represent the geometry and enable the physical model and the

numerical algorithm to resolve the relevant flow physics to the required accuracy.

For most users of CFD, especially those who do not have access to the source code,

the mesh and the time-step size are the only part of the solution procedure that the user
has full control. The importance of the mesh cannot be over emphasized. The mesh must
represent the geometry with sufficient detail and enable the algebraic analog of the
governing PDEs to resolve the relevant flow physics. For complicated steady and
unsteady, three-dimensional problems, the number of grid points or cells that can be used
in a mesh is restricted by either the available computer resource or a need to have a
practical tum-around time in computing a solution. With a constraint on the number of
grid points or cells, accuracy demands grid points to be placed in regions where they are
most needed to resolve the geometry and flow physics (e.g., by r- or h-refinement).
Unfortunately, this non-uniform distribution can create what are referred to as poor-
quality cells, which can induce considerable errors in the computed solutions (see Carey
(1997)). Also, even when using the solution to adapt the grid as optimally as possible,
the solution generated is often far from being grid independent because of resource

constraints.

Thus, the questions are (1) what can one believe in the computed solutions? (2) What
is the error bound on the solutions from a poor-quality or an insufficiently fine grid? This

dissertation addresses these issues.

1.2 Literature Survey

A number of investigators have developed and evaluated ways to quantify errors in
solutions of PDEs that arise from poor-quality or insufficiently fine grids/meshes.
Roache (1998) reviewed and classified all of these methods into two categories: methods
based on multiple grids and methods based on a single grid. A brief review of the
methods based on multiple grids is given is section 1.2.1. Single-grid error-estimation
methods can be classified as algebraic or PDEs, which are reviewed in sections 1.2.2 and

1.2.3, respectively.

1.2.1 Multiple-Grids Error Estimation

Richardson Extrapolation is by far the most popular multiple-grid error estimation
method for finite difference and finite volume methods (Richardson, 1910, 1927). It is a
technique for obtaining a higher-order estimation from a series of lower-order discrete
values. Richardson Extrapolation involves the solution of the governing equations on
three grids of different size. The development is based on a Taylor series expansion of
the solution as a function of the mesh size. The first step in estimating the error using
Richardson Extrapolation is to determine the order of convergence of the solver. Then, a
higher-order estimate can be computed from the lower-order discrete values by
extrapolation method. Finally, the error level in the coarse mesh can be estimated by the
difference of the lower-order solution on the coarse mesh and the extrapolated higher-

order solution. This process can continue with more levels of refinement. If the grids are

fine enough, Richardson Extrapolation can be quite reliable, and it has been successfully
used in a wide range of applications, (Roache (1972), Ferziger (1988, 1989, and 1993),

and Demuren and Wilson (1994)).

The intricacies and pitfalls of Richardson Extrapolation have also been extensively
investigated (Celik & Zhang (1995), Celik & Karatekin (1997), Coleman et a1. (2001),
Eca & Hoekstra (2002)). Some of the issues include specifying a characteristic mesh size
parameter, identifying regions of smooth monotonic convergence, and achieving the

numerical solutions in the asymptotic range.

Celik & Zhang (1995) applied Richardson Extrapolation to turbulent pipe flow and
turbulent flow past a backward facing step. They found that careful consideration was
needed in the application of this method. Richardson extrapolation gave good results in
calculating the apparent order of the numerical procedure only in the regions where flow
variables converged monotonically. In those regions, the grid independent results can be

obtained and the discretization error bounds can be calculated.

Richardson Extrapolation was extended to include estimation of uncertainty using the
Grid Convergence Index (GCI) with factor of safety approach (Roache, 1998). Wilson &
Stern (2001, 2002) used the generalized Richardson Extrapolation for multiple input
parameters to get detailed verification and validation for RANS simulation of a naval

surface combatant.

 

Ilinca, et al. (2000) and Zhang, ct al. (2001) compared different error estimation
techniques for finite-volume solutions, including solution reconstruction, Richardson
Extrapolation and solutions of error equations. They found that for smooth subsonic or
supersonic flows, Richardson Extrapolation provides the more accurate error estimates.
But for flows with discontinuities, Richardson Extrapolation always overpredicts the

exact error, because it assumes smoothness of the solution with the grid refinement.

Freitas, C.J., et al. (2003) reviewed the work on the specification of standards for the
reporting of numerical uncertainty in archival publications and the methods for
quantifying numerical uncertainty. They also presented a five-step procedure for the
numerical uncertainty estimation based on the Richardson Extrapolation and the Grid
Convergence Index (Roache, 1993). The proposed five—step procedure for uncertainty
estimation is as follows:

Step 1: Define a representative grid size

Step 2: Obtain solutions on three sets of grids with significantly different

resolutions

Step 3: Calculate the apparent order of the method

Step 4: Obtain the higher-order solutions by extrapolation using lower-order

solutions

Step 5: Calculate the estimated extrapolated error and the fine grid convergence

index
They also noted that Richardson Extrapolation is far from being perfect, especially for the

regions, where the predicted variables may not exhibit a smooth, monotonic dependence

 

on grid resolution, and for the time-dependent calculations, the non-smooth response can

also be a fimction of time and space.

Another issue in the application of Richardson Extrapolation is in calculating the
apparent order of the numerical method. For PDE problems with inexact or varying
convergence order, the order of a numerical method often depends on space location and
is not accurately satisfied on different levels of grids used in the extrapolation. Shyy and
Garbey (2003) proposed a more robust method of finding the order of a numerical

method by solving a least square minimization problem on the residual.

In conclusion, methods based on multiple grids require solutions to be generated on a
series of increasingly finer grids (at least three, sometimes as many as four to six different
grids to make sure that the three sets of gird used are in the asymptotic range). These
methods can give definitive statements on errors, but is more expensive in requiring
solutions at finer meshes, which can be prohibitive for realistic engineering problems. If
Richardson’s extrapolation is used, then the meshes must be sufficiently fine for the

truncated Taylor series to be bounded.

1.2.2 Single-Grid Algebraic Error Estimation

The algebraic single-grid error-estimation methods can further be categorized into

error “indicators” and error “estimators”. The error “indicators”, also referred to as grid-

quality measures, are developed to evaluate the quality of the meshes. Many error

 

indicators have been developed to evaluate the quality of the meshes and to guide the

process of solution adaptive grid generation.

Most grid—quality measures developed so far only account for the geometry of the
cells in a mesh. Many measures were developed to evaluate the shapes of two-
dimensional elemcnt. The aspect ratio and angular measure were widely used to evaluate
the triangle qualities (Mitchell et a1. (1971)). The measures of triangles were then applied
to quadrilaterals (Robinson (1987, 1988)). Robinson summarized a set of shape or
distortion parameters for the four-node quadrilateral, which are aspect ratio, warpage,
skew angle and tapers. Among these measures, aspect ratio is described as the ratio of
maximum side to minimum side of a planar quadrilateral mesh. The warpage of a
quadrilateral is measured by its deviation from a flat plane. Skew angle and tapers are

distortion parameters and are measured by their deviation from square.

Liu and Joe (1994)) and Lo (1997) extended the grid-quality measures to the tlrree-
dimensional tetrahedral elements. Robinson (1994) extended his two-dimensional
measures to three-dimensional hexahedral elements, including aspect ratios, skew ratios,

tapers and warpages. Some recent studies worked on the geometric measures for 3D

hexahedral, pyramid, and wedge type elements (Owen (1999), Wa and Chen (2000)).

Field (2000) reviewed the geometric measures and proposed the following criteria for
a good measure:

0 An ability to detect all the de-generate elements

o Non-dimensionality, which can make the triangle/tetrahedral and
quadrilateral/hexahedral have the similar dimensionless value.

0 Boundedness, which can limit the measure with arbitrarily large values

0 Normalization, which will let the value of measure take positive value between

zero and one for convenient and general comparison.

Compared to the research on geometric measures, relative fewer investigators worked
on grid-quality measures that can account for both geometry and solution. The early
work for this study was done by Zlamal (1968) and Babuska and Aziz (1976). They
studied the minimum and maximum angles in different flow conditions and presented the
requirements for the triangular element angle based on the approximate finite-element
solutions. Krizek (1992) extended the angle requirement to 3-D tetrahedral elements.
Though a good triangulation should have as few anisotropic triangles as possible, Rippa
(1992) showed that the highly distorted element aligned with the flow direction might be
very effective for the strongly anisotropic solution situations. He also set up a connection
between the shape of the triangle and the behavior of the approximated solution. It was
extended to include adaptation to align the grid with a vector field by Brackbill (1993)
and Knupp (1996). Sorokin (1994) and Lepape and Jacquotte (1994) developed grid-
quality measures which combined geometry and solution gradient. Berzins (1998, 1999
and 2000) showed that error indicators which can account for both geometry and solution

might perform better than the geometric-only measures in evaluating the mesh qualities.

In conclusion, grid-quality measures can be classified into two broad groups: those

that involve information about the mesh geometry only, and those that involve
information about the mesh and the solution generated on the mesh. Each of these two
groups can further be classified as local (involving only one cell) and macro (involving
the cell and adjacent neighboring cells). Grid-quality measures that involve geometry
only and are local include edge lengths, volume size, aspect ratio and orthogonality or
skewness. Grid-quality measures that involve only geometry but are macro include
smoothness. Grid quality measures that involve geometry and solution and are local
include those that measure resolution such as y+ value of the first few grid points normal
to the wall. Grid quality measures that involve geometry and solution and are macro
include smoothness in conjunction with a solution variable. Comprehensive reviews on
error indicators used to evaluate the mesh can be found in Carey (1997) and Roache

(1998)

The usefulness of grid-quality measures for evaluating 2-D or 3-D structured and
unstructured meshes and its encouraging evidence have been reviewed as in the above.
However, the grid-quality measures do not directly correlate solution errors with grid
qualities, i.e., they do not give an estimation of the numerical error of the solutions; hence
the use of term “Surrogate Estimators” and “single grid error indicator” by Roache

(1998).

Many algebraic error “estimators” have been developed for the quantification of the
solution errors. These algebraic error estimators assume that the error at a grid point or

cell is a function of the grid in question and the solution there. Most of work on single-

10

grid error estimation has been done in relation to the solution-adaptive mesh refinement

(Peraire et al. (1976) , Zienkiewicz & Zhu (1987, 1992), Sonar (1993), Oden (1994),

Mackenzie et al. (1994), Carey (1997), Venditti & Darrnofal (2000)).

Within the fiamework of finite element methods, numerous error estimators that

account for both the geometry and the solution have been developed by a number of

investigators. These single-grid error estimation strategies can be classified into the

following different techniques (Oden et al. (1989), Verfurth (1994)):

Residual based methods

This error estimation method can give a strict upper bound for the solution error
in the energy norm and is very popular in recent error estimation in the finite
element method. (Ainsworth and Oden (1992), Strouboulis and Hague (1992)).
Duality methods

This approach was developed based on the Synge’s ‘hypercircle’ in mathematics
and Ekeland’s study on convex analysis and variational problems (Synge (1957),
Ekeland and Teman (1976)). The duality theory of convex optimization is used to
compute the error bounds.

Subdomain-residual methods

This method was proposed by Babuska and Rheinboldt (1978) for bilinear finite-
element approximations. For this method, the local error in an element can be
obtained from a patch of elements surrounding the given element. This method is
very similar to the well-known residual method and the main difference lies in the

fact that they used the norm of negative Sobolev spaces and the error estimates

ll

 

are given in an asymptotic form where the element size approaches zero.

0 Interpolation methods
For this method, the interpolation theory of finite element is applied to generate
the estimation of local error for each element based on the leading term of the
truncation error.

0 Using solution of local problem
For this approach, the original problem was transformed to a similar and simpler
locally discrete problem. The approximate solution of the original problem is
then computed on the locally simplified problem and used to get the estimated
error (Bank and Weiser (1985), Oden et al. (1989), Bank and Welfert (1990)).

0 Sharp a priori error estimation method
This error estimation method is often called “sharp a priori error estimator”
(Eriksson (1985), Eriksson and Johnson (1988)). In this method, the higher-order
difference schemes are used based on the computed numerical results to improve

the accuracy in predicting the derivatives in the priori error estimators.

The methods mentioned above relate mainly to the linear elliptic problems. For the
area of nonlinear problems or parabolic and hyperbolic problems, the algorithms of error
estimation are less well understood. Some error estimation methods for parabolic and
hyperbolic equation using finite element methods can be found in Davis and Flaherty
(1982), Eriksson and Johnson (1991), Adjerid et al. (1999), Johnson et al. (1995), Suli

(1998) and Barth and Larson (1999).

12

 

The algebraic single-grid posteriori error estimation for finite-difference and finite
volume methods is still in its early development stage compared with the numerous
developed error estimation approaches for finite element method. Haworth et al. (1993)
presented a local error estimation approach based on the cell-to-cell imbalances in kinetic
energy and momentum and used it on a transient flow problem in an axis-symmetric IC
engine flow. Jasak and Gosman (1998, 1999) proposed a residual error estimation
method for the finite volume method. This method borrowed the idea from residual
method for the finite element method, in which the error was estimated from the

inconsistency between control volume integration and the face interpolation value.

The single-grid algebraic error estimators typically only consider gradients of the
scalar fields. The CFD solution has vectors and tensors, and the different components of
the vector and tensor fields should be considered collectively instead of individually. For
example, a highly anisotropic mesh can be optimal if all cells are aligned with the flow

direction, which requires a link between the cell geometry and a vector quantity.

Shih et al. (2000), Gu et a1. (2000) and Gu & Shih (2001) proposed grid-quality
measures that also account for the vector and tensor nature of the flow field and link the
solution to the geometry and size of each cell in a grid. Among these measures, one is
based on a length scale representing the curvature of velocity field in three-dimensional
and an effective cell dimension in the direction of change. The second is based on the
directional change in the velocity vector. The third is based on the spectral radius of the

deformation tensor and a cell dimension along its eigenvector. The fourth measure

13

‘ blnu1—1

 

‘I

gauges the smoothness of cell dimensions along the flow direction. The measures
developed were evaluated by establishing and interrogating a database made up of a few
test problems with complex flow patterns. The evaluation showed the measures
developed to have correlations with the solution error defined with respect to the grid-

independent solution.

However, a serious deficiency of the single-grid algebraic error estimation methods is

that they assume the error at a grid or cell is a function of the local grid and local solution.

Babuska et al. (1994 and 1997) recognized that solution error is not a function of the
local grid and local solution and proposed the need for the error transport equation. The

single-grid PDE error estimation methods are reviewed in the next section.
1.2.3 Single-Grid PDE Error Estimation

Single-grid PDE error estimators, first proposed by Babuska et al. (1994, 1995 and
1997), recognize that errors once generated can be transported by advection and diffusion
to other parts of the flow field. Thus, solution error is not a function of the local grid and
local solution. So, a transport equation is needed to understand the generation and

evolution of errors.
A method for deriving transport equations for error was presented by Babuska et al.

(1994, 1995 and 1997) for finite-element methods. Ferziger (1993), Van Straalen et al.

(1995) and Zhang et al. (1997, 2000 and 2001) applied it to finite-difference and finite-

14

3‘ 1'?- m-I-rfl

 

volume methods. Their method proceeds as follows. Consider a differential operator L
operating on dependent variable U:

L(U) = f (H)
where f contains the non-homogeneous terms. If UI, is an approximate solution, then its
substitution into Eq. (1.1) will produce a residual R since it will not satisfy Eq. (1 .1); i.e.,

L(Ua)—f=R (1.2)
If Eqs. (1.1) and (1.2) are linear or linearized, then subtracting Eq. (1.2) from Eq. (1.1)
yields the following transport equation for error:

L(e) = -R (13)
Where

e = U - U. (1.4)
is the solution error. Since the error in the solution is defined with respect to the exact
solution of the PDE, the error-transport equation given by Eqs. (1.3) and (1.4) can

account for both errors from the grid and errors from the numerical scheme.

The main difficulty in implementing the error-transport equation given by Eq. (1.3) is
how to model the residual, which represents the local source of error generation. Zhang
et al. (1997) proposed a method for approximating the residual based on the flux
difference across the grid cell interfaces for a particular numerical method based on

Roe’s upwind flux-differencing scheme.

The modified equation of Warming & Hyett (1974) has been used to guide the

modeling of the residual. The modified equation is derived by first expanding each term

15

of a difference scheme in a Taylor series and then eliminating time derivatives higher
than first order by the algebraic manipulations. The modified equation is the PDE truly
solved by the FD/FV equation. The difference between the PDE and the modified
equation represents the error in the numerical scheme and the mesh/time-step size in

representing the PDE.

Zhang et al. (2000 and 2001) showed that modeling the residual by using the leading
terms of the truncation error in the modified equation can give excellent results. But they
studied only cases where the grid spacing or cell sizes are small. So far, only meshes
with equally-spaced grid points have been considered (i.e., have only investigated errors
induced by inadequate resolution, but not mesh quality). Also, most studies have been

based on 1-D model equations.

Recently, Celik et al (2003) modeled the residual by expanding terms in the FD/FV
equations about the cell center and keeping the leading terms. This modeling approach is
essentially the same as using the modified-equation except time-derivatives are not

replaced by spatial derivatives.

It should be noted that Eqs. (1.3) and (1.4) are only valid for finite-expansion
methods such as finite element and spectral, but not for collocation-type methods such as
finite-difference (FD) and finite-volume (FV). This is because for FD and FV methods,
the differential operator L in Eq. (1.2) is replaced by a discrete operator so that

subtracting Eq. (1.2) from Eq. (1.1) will not yield Eqs. (1.3) and (1.4). This

16

 

 

 

 

inconsistency was also noted by Roache (1998). To rectify this inconsistency, Roache
(1998) suggested all FD and FV solutions at grid points or cells to be made into
continuous functions. However, Roache (1998) recognizes that the resulting continuous

function is non-unique and will depend on the interpolants used.

Thus, an error transport equation based on the discrete finite-difference (FD) and
finite-volume (F V) framework is needed to understand the generation and transportation
of solution errors in CFD. More effective methods for modeling the residual in the error
transport equation are also necessary for more accurate and effective estimation of the

solution errors.

17

 

 

 

1.3 Motivations and Objectives

Among the methods developed so far to quantify grid-induced errors in solutions of
PDEs, the methods based on multiple grids require solutions to be generated on a series
of increasingly finer grids and can give definitive statements on errors, but are more
expensive in requiring solutions at finer meshes, which can be prohibitive for realistic
engineering problems. The single-grid algebraic error estimation methods have a serious
deficiency, i.e., the assumption that the error at a grid or cell is a function of the local grid
and local solution, which is shown to be wrong. The error transport equations developed
so far only work for finite-expansion methods such as finite element and spectral

methods, but not for finite-difference (FD) and finite-volume (FV) methods.

Therefore, the objectives of this dissertation are as follows:

1. Develop an error-transport equation that can be used to estimate grid-induced
errors in finite-difference (F D)/finite-volume (FV) solutions.

2. Evaluate the usefulness of the discrete-error-transport equation in predicting grid-
induced errors.

3. Use the error-transport equation to study the effects of both mesh
spacing/resolution and mesh smoothness on solution accuracy.

4. Examine the usefulness of the modified equation in modeling the residual.

5. Examine the alternative approach for modeling of the residual in the discrete

error-transport equation

18

 

 

 

The remainder of this dissertation is organized as follows: Chapter 2 presents the
concept and derivation of a discrete-error-transport equation for estimating grid-induced
errors in FD and FV solutions. The approach to evaluate the discrete-error-transport
equation is also presented in this Chapter. Chapter 3 evaluates the usefulness of the
discrete-error—transport equation in predicting grid-induced errors and the modified
equation in modeling the residual by applying them to a few test problems. An
alternative approach for modeling of the residual in the discrete-error-transport equation
is examined in Chapter 4. A system of discrete-error-transport equations (DETEs) is
developed for the set of Euler equations and evaluated on the test problems governed by
the Euler equations in Chapter 5. Chapter 6 evaluates the usefulness of the Euler DETEs
to estimate the grid-induced errors for solutions obtained using Navier-Stokes solvers.

The conclusion and future work of this research are discussed in Chapter 7.

l9

 

 

 

CHAPTER 2

DEVELOPMENT OF DISCRETE-ERROR-TRANSPORT
EQUATION

In this chapter, a discrete-error-transport equation (DETE) that can be used with
finite-difference/finite-volume methods is developed. The concept of error transport is
illustrated by a 2—D laminar flow problem. The concept and derivation of the DETE is
demonstrated by a simple model equation: the one-dimensional advection-diffusion

equation.

2.1 Error Transport Concept

In this section, the description and formulation of a 2-D test problem will be
presented. The details on how to get the grid-independent solution are described next.
This is followed by the illustration of how the error is transported in the flow domain for

the case with an imperfect grid.

2.1.1 Problem Description

The problem chosen is a 2-D square chamber (300 x 300 mm) with two inlet ducts

and one exit duct, a schematic diagram of which is shown in Figure 2.1. The inlets are

two straight ducts (200mm x 50 mm). The exit duct has the same dimension as that of

the inlet ducts. The geometry is symmetric about the x-axis.

20

 

 

The fluid in the chamber is water with density p = 998.2 kg/m3 and dynamic viscosity

p: 1.003 x 10'3 kg/(m.s). The inlet velocity is uniform with u = 0.0048 m/s so that the
Reynolds number Re = 222 = 240 and the flow remains laminar. The pressure outlet is
u

fixed with the static pressure equal to the standard atmosphere pressure P0. The no-slip

boundary conditions are employed at the walls.

Because of the symmetry in this problem, only half of the domain is simulated. The
steady-state solution of this incompressible laminar flow is obtained using FLUENT.
The governing equations for this flow problem are the continuity and momentum

equations (full Navier-Stokes equations). Because u is a constant and does not depend

on the temperature, the energy equation is decoupled from the continuity and momentum
equations and is not used in this study. The details of the governing equations, boundary
conditions, numerical schemes and control parameters can be found in the FLUENT

manual.
2.1.2 Grid-Independent Solution Generation

The detailed process to obtain the grid-independent solution for this test problem is as
follows (see Gu (2002) for more details):
1. Generate an initial uniform grid with unit aspect ratio and obtain its solution
(referred to as level k grid/solution).
2. Generate a finer uniform grid with unit aspect ratio and obtain its solution

(referred to as level (k+1) grid/solution)

21

3. Interpolate the level (k+1) solution on the level (k+1) grid onto the level k grid
using bilinear interpolation.

4. Compute the average relative error between level k and level (k+1) solutions by

i Vim —Vk+l,i (2'1)
—k is] Vin-rm

Ere:
N

 

 

 

 

 

where V“ is the velocity vector at the ith cell on level k mesh, Vim is the level (k+1)
solution projected on the level k mesh at the ith cell, and N is the total number of cells in

the level k mesh.

5. If E; is less than the lower bound 8 , then the level (k+1) solution is taken to be

the grid-independent solution. Otherwise, repeat steps 2 to 4.

The series of uniform meshes used to generate the grid-independent solution of this
test problem is summarized in Table 2. 1. The lower bound 8 is taken to be 0.5% for this
test problem and the solution on level 8 grid is taken to be the grid-independent solution.
Figure 2.2 shows the grid-independent solution represented by streamlines and contour
plots of x- and y- velocities. Because of the symmetric geometry and boundary

conditions, the flow pattern is symmetric about the x-axis.

2.1.3 Imperfect Grid and Error Transport

With the grid-independent solution generated, we now perturb the grid at one

22

location, making it of poor quality. Except for this location, the grid elsewhere is perfect.
The perturbed grid is shown in Figure 2.3. The solution is also obtained for this

imperfect grid using FLUENT.

The absolute error and the relative error of the ith cell on the imperfect grid are

 

 

defined as follows:
Eabj = Von 'Vi‘ (2'2)
Von “Vi, (2'3)
Ere,i = ——-—
VGIJI

 

where Von is the velocity vector at the ith cell on the grid-independent mesh, and Vi is

the solution on the imperfect mesh at the im cell.

The absolute error and relative error between the grid-independent solution and the
imperfect-grid solution are shown in Figure 2.4. It can be seen that the error is not
confined to the region where the cells are of poor quality, but rather it gets transported
throughout the flow domain. Since the grid is perfect elsewhere, it is very clear that error
is not a function of local grid and local solution. If error is not a function of local grid
and local solution, an error transport equation is needed to estimate the error distribution
in the flow domain. The concept and derivation of the discrete-error-transport equation
(DETE) for the finite-difference and finite-volume methods are presented in the next

section using a simple model equation.

23

 

2.2 Discrete-Error-Transport Equation (DETE)

2.2.1 Description of Model Equation

The model equation used to present and evaluate the error-transport equation is the

following one-dimensional (l-D) advection-diffusion equation:

 

dU d’U (2-4)
C—=v
dx dx2
Or
(1 d (2.5)
L =0, L=— C—v—
(U) dx( dx)

In the above equation, C, the advection speed, and v , the kinematic diffusivity,
are constants. If the boundary conditions at x = x0 = O and x = t’ = 1 are U(xo) = 1 and

U( l) = 2, respectively, then the exact solution of Eq. (2.4) is

x—x0)_l (2.6)

 

exp(Re
U(X) - U(Xo) _
U(t’) — U(x o) - exp(Re) —l

 

 

where Re is the Reynolds number, defined by

Re=C€/v (2.7)

2.2.2 Numerical Method of Solution

To present and evaluate the error-transport equation, we consider only second-order

central differencing. When the mesh is uniform (i.e., Ax = 8/0 - 1) = constant, where I =

24

number of grid points), Eq. (2.4) becomes:

Um - Ur r Um - 2U; + Ui-l (2°83)

C———'—=v 2
2Ax Ax

 

i=2, 3, I, Ax =e/(1-1)
U1=1 (2.8b)

UI = 2 (2.8c)

Equation (2.8) forms a linear system of equations with a tridiagonal coefficient
matrix, which can be solved very efficiently by the Thomas algorithm. Since central

differencing is used, the coefficient matrix loses diagonal dominance when the cell

 

Reynolds number, Re A = C Ax / v exceeds 2.

Though the exact solution exists in Eq. (2.6), a series of solutions were generated by
using uniform grids with increasingly finer mesh to obtain the grid-independent solution.

A grid-independent solution was obtained for the following parameters:

C= 10 m/s (2-9)
v=1 (m/s)2
Z =1 m

The total number of grid points used was 1,001, so that Ax = 0.001 and Re A = 0.01.
This solution coincides with the exact solution given by Eq. (2.6), (the maximum percent
error is less than 0.005%). This grid-independent solution was used to evaluate the
methods presented here. The agreement between the grid-independent solution and the

exact solution is a check on the numerical method of solution.

25

2.2.3 Discrete-Error-Transport Equation (DETE)

In finite-element methods, the solution U is replaced by a finite expansion in terms of
known functions (e.g., U = ZciBi ). Thus, the same differential operator is used in

generating the exact solution and the approximate solution. For such methods, Eqs. (1 .1)

to (1.4) is the correct way to define the residual and the solution error.

In finite-difference and finite-volume methods, the differential operators are replaced
by difference operators as shown by Eq. (2.8). For these methods, Eqs. (1.1) to (1.4) are
not meaningful because the differential operator and the difference operator are not the
same. Thus, a discrete version of the error-transport equation is needed. Its derivation

proceeds as follows.

First, we re—write the difference operator given by Eq. (2.8a) as follows:

LDUi = o (2.10a)

C v 2v C v
L = ————— E"+—E°+ ——-— E" 2.1%
D [ 2h hi] h2 [2h hi] ( )

where h is the grid spacing (i.e., h = Ax ), and E is the shift operator (i.e., EMUi = Um ).

For a fixed C and v , Eq. (2.10) shows that the difference operator is a function of h

(i.e., LD = LD (h) ). Thus, for the grid-independent solution with hg = 0.001m, we have

26

 

 

 

 

L,,(h,)Ug,i =0 (2.11)
For a solution on a coarser mesh h (hc > hg), we have

LD(hc) U,i =0 (2.12)

From Eqs. (2.11) and (2.12), it can be seen that two different discrete error-transport
equations can be derived, depending upon which h is used, hc or hg. If hg is used, then we
have

LD(h,)U,,. =0 (2.13a)
LD(hg)Uc,i =Rg (2.13b)
Subtracting Eq. (2.13b) from Eq. (2.13a) gives

Ln(hg) e.- = R8 (2.14)
61 = U33 — Ucfi' (2.15)

The parameter c; is the error in the solution ch at grid point i, defined with respect to
the grid-independent solution Ugj. If hc is used, then we have
LD(hc) U“,i = Rc (2.16a)
LD(hc)Uc'i =0 (2.16b)
Subtracting Eq. (2.16b) from Eq. (2.16a) gives

Ln(hc) 3i = Re (2.17)

Note that in Eqs. (2.13) and (2.16), LD(ha,i) U13; = 0 if and only if or = B (e.g., Eqs.

(2.13a) and (2.16b).

27

 

 

 

Equations (2.14) and (2. 17) are the two discrete error-transport equations. In practice,
the residual in these two equations must be modeled, a subject to be further discussed in
Chapters 3 and 4. In the remainder of this chapter, we want to show how to compute the
residual if the grid-independent solution is known and to examine, which of the two

error-transport equations is the correct one to use.

In order to calculate the residual for the discrete error-transport equation given by Eq.
(2.14), one must interpolate the coarse-grid solution Uc to the grid-independent grid with
spacing hg. This can be problematic since there are many grid points between two
successive coarse grid points. If linear interpolation is used (see Figure 2.5), then we will
have a residual that is step firnction. If higher-order global or spline polynomials are
used, then the residual computed is a function of the interpolant used, and hence is not
meaningfirl. Numerical experiments performed indicate that this is precisely the case
(see Figure 2.6). Thus, though the discrete error-transport equation given by Eq. (2.14)
is essentially the same as the one defined by Eqs. (1.3) and (1.4), it is inappropriate for

finite-difference and finite-volume methods.

In order to compute the residual for the discrete error-transport equation given by Eq.
(2.17), one must interpolate the fine-grid solution Us to the coarse grid with spacing he.
Since the number of grid points in the grid-independent grid far exceeds the number of
grid points in the coarse grid, interpolation in this case is trivial. In fact, linear

interpolation is adequate, and is the interpolant used here (see Figure 2.7).

28

 

 

 

Thus, the error-transport equation given by Eq. (2.17) is the more appropriate one. It
also turns out that this is the more natural one to use since it is the solution on the coarse
mesh that we seek for an error estimate. Since Ugj is not unique in the strict sense and
may differ from the exact solution because of spurious modes permitted by FD/FV
methods used, Eqs. (2.15) and (2.17) can only account for the grid-induced errors, but

not for errors fi'om the numerical scheme.

The derivation of the DETE also illustrates the approach in implementing the DETE:
a grid-independent solution is computed for the problem. This solution is taken to be the
accurate baseline or reference solution. Next, a series of approximate solutions are
generated based on the imperfect meshes by using the same flow model and the same
numerical scheme. Then the difference between each approximate solution and the grid-
independent solution can be defined as the exact error, which is induced by the imperfect
grid. Then, a linear or linearized discrete-error—transport equation is derived and the
predicted error is obtained by using Eq. (2.17). The predicted error is then compared
with the actual error to evaluate the usefulness of the DETE in estimating grid-induced

errors.

The main difficulty in implementing the DETE is how to obtain the residual in the
error equation. If the grid-independent solution is known, then the residual at every grid
point or cell can be computed by using Eq. (2.16a) and the error can be predicted using

Eq. (2.17). But, since the grid-independent solution U“ is generally unavailable, a

method is needed to model the residual. Zhang et a! (2000, 2001) suggested

29

 

 

 

 

 

modeling the residual by using the leading term of the modified equation. The modified
equation is the PDE truly solved by the FD/FV equation. The difference between the
PDE and the modified equation represents the error in the numerical scheme and the
mesh/time-step size in representing the PDE. Since the discrete-error-transport equation
addresses only grid-induced errors, the modified equation contains too much information.
Nevertheless, this modeling approach is used and evaluated in Chapter 3. An alternative

approach for residual modeling based on data mining is evaluated in Chapter 4.

30

 

 

 

No-Slip Wall

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

D Symmetry Plane
V l'/ h Al)’ x Z/ L
eocrty Pressure
............... —’—.—o—.—.—.—.—.—.—.—. .—.—.—.—.—
Inlet \I 1: Outlet
1~———200 300 200—————~

 

 

 

 

 

 

Figure 2.1 Geometry and flow conditions of the 2-D chamber problem

Table 2. 1 Uniform grids used to get grid-independent solution for the 2-D chamber
problem

 

 

 

 

 

 

 

 

 

Level Grid Spacing (m) Total Cell Number

1 8.333 x 10'3 (1/120) 864

2 6.250 x 10‘3 (1/160) 1,536
3 4.167 x 10‘3 (1/240) 3,456
4 3.125 x 10'3 (1/320) 6,144
5 2.083 x 10‘3 (1/480) 13,824
6 1.563 x 103(1/640) 24,576
7 1.250 x 10'3 (1/800) 38,400
8 1.042 x 10'3 (1/960) 55,296

 

 

 

 

 

31

 

 

V (ms)

0.001
0.000
0.001
0.002
0.003
0.004
0.005

  
   
  
 
   

 

 

(b) flow pattern in the domain

Figure 2.2 Grid-independent solution for the 2-D chamber problem

 

 

‘3
“ll “———1

T W

 

(a) location of perturbed grids

 

(b) local view of the poor-quality cells

Figure 2.3 Perturbed grids at location (h,h) for the 2-D chamber problem

 

E .1 (In-s)

 

 

>.
I

O 3
' 4.151502
>' 3.69E-02
3.23502
0.2 2.77 15-02
1 2.31 E-02

 

 

(b) relative error in velocity as defined in Eq. (2.3)

Figure 2.4 Error distribution in the flow domain for the 2-D chamber problem with
perturbed grid

34

a
m 400

-800

 

A L
interpolant

    

 

 

 

 

Figure 2.5 Interpolating from coarse to fine mesh.

 

 

 

 

 

 

Figure 2.6 Residual (Rg) computed for Eq. (2.14).
Lefi: linear interpolant; Right: cubic interpolant.

35

-. 3 ,,,,,,,,, :1 WI 11“ 1'2 - 3
I ll ‘3 4‘ ‘ ‘3‘
‘ 0'8 _ .L ‘
“.5
L “a 0.4 ~ "l 3
l3 o_o ___.... ll ‘ l '1‘ l“ “‘ 3‘3:
_ 1 ‘ . l - .
0 0.5 1 0 0.5
X X

 

 

 

 

 

0.06 -
00.04-
0:
0.02 -
L
0.00 .
0 0.5 1
x

Figure 2.7 Residual (RC) computed for Eq. (2.17)

36

CHAPTER 3

EVALUATION OF DISCRETE-ERROR-TRANSPORT
EQUATION ON MODEL EQUATIONS

In this chapter, the usefulness of the discrete-error-transport equation in predicting
grid-induced errors is evaluated by applying it to a number of test problems by using four
one-dimensional PDEs: the advection-diffirsion equation (linear and steady), the wave
equation (linear and unsteady), the inviscid Burgers equation (nonlinear and unsteady),
and the steady Burgers equation (nonlinear and steady). For the two PDEs that are not
linear, a non-conservative and a conservative linearization procedure is proposed and
examined. For all test problems, the residual is modeled by the leading term of the
remainder in the modified equation. Also, for all test problems, grid-independent
solutions were generated to assess how well the residuals are modeled and how well grid-
induced errors are predicted. For the ID Advection-Diffusion problem, the effects of
both mesh spacing and mesh smoothness on the residual and the solution error were also

examined

This chapter is organized as follows. First, the discrete-error-transport equations for
the four test problems are derived. Then, modeling of the residuals by the modified
equation for the test problems is presented. This is then followed by an evaluation of the
usefulness of the discrete-error-transport equation and the ability of the modified

equation in modeling the residual for the four test problems.

37

3.1 Derivation of Discrete-Error-Transport Equations

3.1.1 Advection-Diffusion Equation (Linear and Steady)

The model equation of 1-D advection-diffusion equation and its DETE derivation
have already been discussed in Chapter 2. To illustrate the effects of resolution and
smoothness on residual and solution error, the DETE on non-uniform meshes is derived

as follows.

Finite-Difference/Finite-Volume Eguations on Non-Uniform Meshes

To evaluate mesh resolution and mesh smoothness for a 1-D problem, we define the

following smoothness parameter:

x. —x. Ax.
‘=-.:”- <30
i i-l '

 

If S is a constant (and recall that Z =1), then

Ax2 +Ax3+---+Ax, =1
Ax2(1+S+S2+---+SH)=1 (3.2)
Ax2 (1 -s"')/(1 —S) =1

Thus, once the smoothness parameter S and the number of grid points I are specified,
then Ax2 can be computed by using Eq. (3.2) and the location of all grid points are given
by

x, = Ax, (1+S+S2---+S“2) =(1-si-‘)/(1-s‘*'),

(3.3)
i= 1, 2, 3, ...,I

38

To facilitate analysis, the non-uniformly distributed grid points in the physical

domain (x) are mapped to a transformed domain (2‘, ), where all grid points are equally
distributed; i.e.,

éi =(i—1)A§, A§=1/(I—1) (3.4)
Combining Eqs. (3.3) and (3.4) gives the following functional relationship between x and
E:

xi = (1 — S““"‘ )/(1 - s“) (3.5)

The model equation given by Eq. (2.4) is mapped to the transformed domain by

noting

_d__§_d_ 36
d5"- xdé ()

fi'lfi

1
dx

where the metric coefficient fix is obtained by using Eq. (3.5). The model equation in

the transformed domain is
d dU
Cd§”§vc1‘§‘§ at“é (3'7)

By using second-order central differencing and rearranging, the above equation becomes

 

 

 

1 v v v C
__.___ U._ + + U- + U 0
Axi+l+AXi [(C 2 Axi) H (Ax, A in) ' (2- vi“) M]: (3.83)
2

,i=2, 3, ..., L1,
or

LDNU1=O

39

 

 

1 C v v v C v
L =——— ———-——-—-— E"+ -—-—+ E°+ —— E+1
D. _ ‘< 2 Ax) <Ax A ) <2 Ax333) ‘

r i i+l

i= 2, 3, 1-1, (3.8b)

Note that if S = 1, then the above equation reduces to Eq. (2.10a).

Discrete Error-Transport Equation on Non-uniform Meshes

The discrete error-transport equation is obtained by using the procedure outlined by
Eqs. (2.16a) and (2.17), which for Eq. (3.8) becomes

LDNGI) 61 = R1 (3.9)

where h is the coarse mesh, and R, in the above equation still needs to be specified. For

this model equation, its actual value can be computed by R; = LDN(h) Us; with linear

interpolation of the grid-independent solution Us) to the coarse mesh. This computed

value will be used to assess the values from a model of the residual.
3.1.2 Wave Equation (Linear and Unsteady)

The wave equation that describes a wave propagating in the x direction with constant
speed C is given by

6U+C6U_

E 332—0 (3.10)

In this study, C was set equal to 2, and the following sinusoidal waveform is being
advected:

U=1.5+0.5sin(21tx—1t/2),for 03x31 (3.1121)

40

 

 

 

U = O , for x > 1
The exact solution of Eq. (3.10) with the above wave form at time t = 0 is

U = 1.5 + 0.5 sin[21t(x — Ct) —- 1t / 2]

By using the Lax-Wendroff scheme, Eq. (3.10) becomes

U?“ ‘U? +C U11 " Uli-r = Czk 3 U?” “2U? + U11

k 2h 2 h2

 

i=2, 3, ..., I-l; n=0, 1, 2,....
which in operator form is

LD(h,k)U§' = 0, i=2, 3, ...,I-1;n=0, 1,2,...

Czk C re2
+———-—— 13“
112) (Zh 2112)

A c kc2
L =—+ -—- E"‘+
D k Zh 2112) (

 

 

In the above equations, h is the grid spacing; k is the time-step size; and

Uf=U(x,,t"), xi=(i—1)h, t"=nk,AU"=U“”—U"

(3.110)

(3.12)

(3.13)

(3.14a)

(3. 14b)

(3. 14c)

Though the PDE given by Eq. (3.10) only permits boundary conditions (BCs) at the

inflow boundary at x = 0, the following weakly redundant BC was imposed at the

outflow boundary x = t’ = 10: 620 / ax2 = 0. This BC did not create a problem because

computations were always terminated before the waveform given by Eq. (3.12) reached

the outflow boundary.

The discrete-error-transport equation for the FD/FV equation given by Eq. (3.14a) is

given by

41

LD(h,k)e§‘ = R? (3.15)
where L1) is given by Eq. (3.14b). Equation (3.15) is derived in a manner similar to the

derivation of Eqs (2.17) and (3.9).
3.1.3 Inviscid Burgers Equation (Nonlinear and Unsteady)

The inviscid Burgers equation is given by

6U+6(U2 /2) _
at 6x

 

0 (3.16)

Aside fiom not being linear, Eq. (3.16) differs from Eq. (3.10) in two ways. First, the
wave speed is not a constant and varies with the solution U. Second, Eq. (3.16) admits

weak solutions. Even initially smooth solutions can become discontinuous in time.

The following sinusoidal wave form is assumed to have entered the flow domain at

timet=0:
U=O.15+0.0581n(21tX-1t/2),fOI'OSXSI (3.173)
U=0,forx>l (3.17b)

With the above waveform, a weak solution will form. The exact solution to Eq. (3.16)
before the weak solution forms is

U = 0.15+0.05 sin[21t(x—Ut)—1t/2] (3.18)

By using the Lax-Wendroff scheme, Eq. (3.16) becomes

| w

UM—UI' F“ —F“ l n ,, 3.1%
r k 1+ 1+12h1-1=§ [Ai+l/2(Fi+l ( )

“Fin)—A?—r/2(Fin _Fi':1)]

2

5"

42

where

A=U, F=U2/2 (3.1%)
tr U? +Ui1+ n U? + U.“—
Am =——2—‘-,.A,_,,, =——§—‘ (3.190)

Because Eq. (3.19) is not linear in U, it needs to be linearized before the discrete-
error-transport equation can be derived. Two different ways to linearize the FD/FV

equations are described below.

Linearization Method 1
Since Eq. (3.16) can be written as

6U 6U 6 6
——+U—=0 LU=0,L=—+U— ,
at 6x °r at ax (320)

through chain rule, the above equation can be used as a guide to linearization.
Comparing Eq. (3.20) with Eq. (3.10) suggests treating the U in the differential operator
L as a variable coefficient instead of the dependent variable. Applying the Lax-Wendroff

method to Eq. (3.20), treating U in the differential operator L as a variable coefficient,

gives
n+1_ n n _ 3""
Ui Ur +(U?)c UI+1 U14
1 11: 2h (3.21)
=2h—2[((U?)C)Z(U§1t —2U§' +U§'_, )]
01'

43

LD(h,k)U§‘ = 0 (3.22a)

= +(L’Ii ___)__hc +1 -I ___1_L n H _ " 3.22b
Lo k+ —(E "E ) 2h2(Ui)°(E 2+E ) ( )

where terms such as (U,u )c with subscript c are treated as variable coefficients.
By using Eq. (3.22), the discrete-error—transport equation is given by

n“- 9 3' _ 3'
ci 61 +(Un) e1+1 e1--l _1_ k

k , c 2h 2112[((U ) )2 (ell—2e? +e?_ l)]= R“ (3.23a)

 

01'
LD(h,k)e{‘ = R," (3.23b)

where R? is the residual, and L9 is given by Eq. (3.22b).

Linearization Method 2

Though Linearization Method 1 appears reasonable, the F D/FV equation given by Eq.
(3.22) and the discrete-error—transport equation given by Eq. (3.23) are both written in
non-conservative form, so that “shock capturing” is not possible. A linearization that

maintains the conservative form, enabling “shock capturing”, is as follows:

U“ U1,

 

 

 

UP“ _U“+ (— 2i——+')c Ui+1— (— )c U?—
k 2h
1 k +U," U" U!‘ +U;' U" (3'24)
=57“ U” ” 2) U..." -(-'—‘5—- 2 ').U.
U“ 2U:I—l UII Un 2U?| UH
_(__ 1 + Ui _2_c‘U) +( 1 + Ui— 21M “1]

 

01'

44

LD(h,k)U? =0 (3.25a)

A +l_ n -l
LD= I +4lh((U:'.. E (U -..) E )
l k [Jill-I'll:I +1_ 11 0
21.2 [—7—] «Ur... E (U.).E)
1k U?+U?-1 0_ —1
+5h—2[—4——]C((U? ) E (U? _,) E ) (3.25b)

Similar to Eqs. (3.21) and (3.22), terms such as (Uf)c with subscript c are treated as

variable coefficients.

By using Eq. (3.25), the discrete-error-transport equation is given by

 

 

 

e?+'e +(Ui+lc e1+1 -(U?.1)ce?.1
k 4h

1 1‘ U?,,+U" UL. U1+1+Usn U." 11
211—2“ 2 )cC" ..1 (— - _Ie)c e1] (3.26a)
1k U_;'__+U;',U_." ,, U?+U?,U :1, -11"

+_ 2112 —[( 2 2 )cei -( 2 2en)c I 1]

or
LD(h,k)e? = R? (3.26b)

where R? is the residual, and L1) is given by Eq. (3.25b).

3.1.4 Steady Burgers Equation (Nonlinear and Steady)

Consider the steady-state Burgers equation given by

45

a (U2 /2 62 U
“7—267 (327)

with v = 0.1, and BCs: U(xo = 0)=1 and U(€=1)= 2.
By using second-order central differencing, Eq. (3.27) becomes

(U2 /2) - (U2 /2).-. _ V U,,, — 2U3. + UH =

2b h2

1+1

0 (3.28a)

 

i= 2, 3, ..., I-l; U1 =1; U1: 2 (3.28b)
The system of equations given by Eq. (3.28) is tridiagonal and can readily be solved by

the Thomas algorithm.

Since Eq. (3.28) is not linear, linearization is needed before a discrete-error—transport
equation can be derived. Following the approaches presented in section 3.1.3, two

linearization methods are used.

Linearization Method 1

Since Eq. (3.27) can be written as

2
Uggzva U
6x 6x2

 

or LU = 0,

(3.29)

2
L=Ui-v 6
6x 6x2

 

By treating the U in the differential operator L as a variable coefficient and using central

differencing yields

(U1). 24% = v 112 , (3.3021)

 

46

i=2, 3,...,I-1; U1=1; U1=2
or

L,,(11)Ui =0

= (Ui)c

L
D 211

+1 -1 V +I —1
(E —E )-h—2(E -2+E )

The discrete-error-transport equation corresponding to Eq. (3.31) is

 

e. —e. e. —2e.+e.
U- £'—"_'__V 1+1 1 1—1 =R-
( ')° 2h h2 '
01'
LD(h)ei :Ri

where L is given by Eq. (3.31b).

Linearization Method 2
A linearization that maintains the conservative form is

(U1+1/2)cU - (Ui-I /2)c Ui-I _ v Ui+l - 2Ui + Ui—l =

2b h2

HI

0,

 

i=2, 3, ..., I-l; U1 =1; U1=2
or

L,,(11)Ui =0

L. = $101me -(U.-. arm-$03“ -2+ E")

The discrete-error-transport equation corresponding to Eq. (3.34) is

47

(3.30b)

(3.31a)

(3.31b)

(3.32a)

(3.32b)

(3.33a)

(3.33b)

(3.34a)

(3.346)

01'

(Um /2)cei+1 -(Ui-l l2)cei-l _
2h

 

LD(h)ei = R1

where L is given by Eq. (3.34b).

i+l

48

(3.35a)

(3.35b)

3.2 Modeling the Residual in the Discrete-Error-Transport Equation

With the discrete-error-transport equations derived, the next step is to model the
residuals in those equations. The residual represents the local generation and destruction
of errors by the mesh and the time-step size if problem is unsteady. Once the residual is
known, then the grid-induced error at every grid point or cell can be computed by the

discrete-error-transport equation derived in section 3.1.

As noted in Chapter 2, if the grid-independent solution is known, then the residual at
every grid point or cell can be computed by using Eq. (2. 16a), which is generalized below
for steady and unsteady problems based on the derivations in section 3.1:

R? =L,,(h,k)U;;‘i (3.36)
In the above equation, LD(h, k) is the difference operator based on h (h > hg) and k (k >
k8) used to obtain the numerical solution whose grid-induced errors are sought, and LD(h,

k) operates on the grid-independent solution UL . But, since the grid-independent

solution U2, is generally unavailable, a method is needed to model the residual.

The modified equation is used as the residual modeling approach in this chapter. The
modified equation is the PDE truly solved by the FD/FV equation. The difference
between the PDE and the modified equation represents the error in the numerical scheme
and the mesh/time-step size in representing the PDE. Since the discrete-error-transport

equation addresses only grid-induced errors, the modified equation contains too much

49

information. Nevertheless, this modeling approach is used here.

3.2.1 Advection-Diffusion Equation

The modified equation for the finite-d1fference/finite—volume equation given by Eq.

(2.4) on non-uniform meshes is

2
ng_vd U

dx dx 2 = RHS (3.37a)

 

d’U A_l‘,2_vd2U diet. A? -11 d’U d6. 4:2 _v d’é. d_U 4:2

 

 

 

 

 

RHS: .C
§‘[ 62:3 6 d§2 d§2 8 dg" d: 6 d§3 d5, 24
A? d‘U .
—v 12 fix dE,‘ ]+h1gher—order terms (3.37b)

In the above equation, the RHS terms represent the spurious physics introduced by the
numerical method of solution. Following Zhang, et al. (2000, 2001), we use the leading

terms on the RHS as a model for R1 in the error-transport equation given by Eq. (3.9); i.e.,

 

 

 

d’UA_§:_ dzU 62:, A62 _vd’U 6;, Ag2

 

 

 

 

R: .c
' 5”? dg’ 6 var} d? 8 dg’ d: 6
3 (3.38)
2 2 4

_v_d £3. EALWIK i. d U]

d5, d5, 24 12 dg“

On uniform meshes, the above equation becomes:
R — 9 d3U -1 d4U 112 339
i61111312611“i (')

For our model equation, R1 is numerically differenced, and the derivatives are all

replaced by second-order formulas (center when possible, otherwise biased upwind). At

50

 

 

 

grid points next to boundaries, one-sided differencing is used.

3.2.2 Wave Equation

The modified equation for the FD/FV equation given by Eq. (3.14) for the wave

equation is
6 U 6 U
6 t +C-- 6 x =RHS (3.40a)

63U +C2

RHS —(-112 422-6710 —k(h2— JZS4TC11) U+HOT (3.406)

The residual in the discrete-error-transport equation for Eq. (3.14) given by Eq. (3.15) is

modeled by the leading terms of RHS; i.e.,

6U +C2

4 n
R? = 9-(112- c k )— ———1<(112 c21<2)a—U (3.41)
6 6x‘ 3

3.2.3 Inviscid Burgers Equation

The modified equation for the inviscid Burgers equation is more complicated, partly
because of the nonlinearity and partly because different linearization methods have
different modified equations. Because of this complexity and because the modified
equation is an imperfect model of grid-induced errors, the modified equation for the
FD/FV equations for the inviscid Burgers is taken to be the modified equation for the
wave equation (Eq. (3.10)), except replace the constant wave speed C by the variable U.

Thus, the residuals in the discrete-error-transport equations for both linearization methods

51

 

— Eqs. (3.22b) and (3.23b) for Linearization Method 1 and Eqs. (3.25b) and (3.26b) for

Linearization Method 2 — are modeled by

USU U2
6x3

 

 

U
Rr=——w—Umz
' (6( ) 6x4

3.2.4 Steady Burgers Equation

4 n
+7iuw—Umbau‘

BAD

The steady Burgers equation is also not linear. For the reasons mentioned in section

3.2.3, the modified equation for the FD/FV equations for the steady Burgers equation is

taken to be the modified equation for the advection-diffirsion equation (Eq. (3.39)),

except replace the constant advection speed C by the variable U. The residual in the

discrete-error—transport equation for both linearization methods — Eqs. (3.31b) and

(3.32b) for Linearization Method 1 and Eqs. (3.34b) and (3.35b) for Linearization

Method 2 — is

 

52

04»

3.3. Results and Discussion

In this section, we evaluate the discrete-error-transport equation derived in section 3.1
and the modeling of the residual by the leading term of the modified equations derived in
section 3.2 for four l-D test problems: the advection-diffusion equation (linear and
steady), wave equation (linear and unsteady), the inviscid Burgers equation (nonlinear

and unsteady), and the steady Burgers equation (nonlinear and steady).

For each test problem, a grid-independent solution is generated along with several
other solutions that are not grid independent. For solutions that are not grid independent,
the following are examined:

1. How well does the leading term of the modified equation model the residual?
This is assessed by comparing the modeled residual with the exact residual
computed by using Eq. (3.36).

2. How well does the discrete-error-transport equation predict grid-induced errors?
This is assessed by comparing the predicted error by using the exact and the
modeled residual (i.e., solving e? by LD(h,k)e? = R?) with the actual error
computed by comparing the actual solution against the grid-independent solution

(i.e., computing e? by e? = U23, —U? ).

For thel-D Advection-Diffusion problem, the effects of both mesh spacing and mesh

smoothness on the residual and the solution error were also examined.

53

3.3.1 Advection-Diffirsion Equation

This sub-section presents results obtained for the advection-diffirsion equation with
the FD/F V equation given by Eq. (3.8), the discrete-error-transport equation given by Eq.
(3.9), and the modeled residual given by Eq. (3.38). Table 3.1 summarizes the S and I
simulated to evaluate the usefulness of the error-transport equation in predicting solution

CI'I'OI'.

Error-Transport Eguation with Actual and Modeled 3,

Figure 3.1 shows the actual residual (denoted as RA) and the modeled residual

(denoted as RM) as a function of x and ‘c', (the transformed coordinate) for two different

A number of grid points (I = 11, 41) and two different smoothness parameters (S = 0.85,
0.95). By comparing the modeled residual (RM) to the actual residual (RA), several

observations can be made. The first is that RA oscillates as 5 = 1 is approached when

the number of grid points was increased from 11 to 41., but RM does not, The oscillation
is generated in the process of linear interpolation of the fine-grid solution to the coarse
grid, when the grid-spacing of the coarse-grid is actually smaller than or near to that of
the grid-independent grid. Figure 3.2 shows that the oscillation of the actual residual for
the I = 41 and S = 0.95 case has been greatly reduced when a finer- mesh- solution (I =
2001) is used as the grid-independent solution. This confirms the idea put forward earlier
in Chapter 2 that the error transport equation should be defined on the coarse grid and the

interpolation should be from fine to coarse grid. For this test problem, the gradient in the

54

solution U increases as E is increased. When S is near 1 (S = 0.95) and the grid is coarse

(I = 11), the shape of RA and RM is similar to the shape of the solution U. But, as S
departs further from unity and as the grid is refined, they no longer have that
resemblance. When the number of grid points is increased from 11 to 41, RM does
approach RA as expected. These comparisons show that the modeled residual (RM) is

qualitatively similar to the actual residual in behavior and magnitude.

Figure 3.3 shows the actual error (e.; ei =Ug3i —U,) in the solution and the error

predicted by using the error-transport equation (Eq. (3.9)) in which two different R1 were
used: RA and RM. From this figure, it can be seen that if the actual R1 is used in the
error-transport equation, then the error is predicted perfectly. This is interesting
especially in view of the large oscillations observed. The results shown here indicate the

error-transport—equation concept to be correct and useful.

Since the actual residual (RA) is typically unavailable, it must be modeled. Figure
3.3 shows the predicted error from the modeled residual to follow the same trend as the
actual error. Thus, qualitatively, the prediction is quite good. Quantitatively, however,
the prediction is less satisfactory. Since the trend is captured, this approach should be

useful for mesh refinement though not necessarily for error estimation.

Results on Effects of Spatial Resolution (S=1)

Figure 3.4 shows the actual residual (RA) and the residual modeled by using the

55

modified equation (RM) as a function of x for two different grids (I = 11 and I = 101).
From this figure, it can be seen that the modeled residual (RM) approaches the actual
residual (RA) as the number of grid points increases. With I = 101, RM and RA are

almost identical. With I = 11, RM differs from RA though the functional form is similar.

Figure 3.5 shows the grid-induced error predicted by using the discrete-error-
transport equation with the actual and the modeled residual shown in Figure 3.4. Also
shown in the figure is the actual error. From this figure, it can be seen that if the modeled
residual (RM) is used, then the error predicted (denoted as E-RM) improves as the grid is
refined. As shown in Figure 3.5, the agreement is excellent for the I = 101 case. For the

= 11 case, there is some discrepancy between the predicted and the actual error.

Figure 3.6 shows the actual and predicted error for a series of solutions with uniform
grid spacing (S = 1) but different number of grid points (I = 6, 11, 21, 101). From this
figure, it can be seen that the predicted error has the same trend as the actual error when
the grid spacing is sufficiently small. When the mesh is extremely coarse, the predicted
error shows an incorrect oscillation. As expected, the error decreases as grid spacing
decreases. Less obvious is the location of maximum error. As the mesh is refined, the
location of the peak error shifts right. When the grid spacing exceeds 0.2 m, the cell
Reynolds ntunber exceeds 2, and this can cause the system of equations given by Eq.
(3.8) to lose diagonal dominance. Losing diagonal dominance could produce some
spurious oscillations in the solutions if the number of equations (i.e., I) is large. Since I is

small, these oscillations did not lead to errors of the opposite sign.

56

Note that only the absolute error is given. Since the solution is between 1 and 2
(imposed by the boundary conditions to the model equation), the absolute error shown

here is effectively the same as the more meaningful relative error.

Results on Effects of Mesh Spacing and Smoothness

Figure 3.7 shows the actual and predicted error for a series of solutions with the same
number of grid points (I = 21) but different smoothness parameters (S = 0.6, 0.7, 0.8,
0.85, 0.9, 0.95). Note that the results are shown with respect to the physical coordinate x

and the transformed coordinate é. We examined only S less than unity because the

solution to our model equation has the highest curvature as x increases so that we need

finer grid spacing there.

From Figure 3.7, it can be seen that the predicted error has the same trend as the
actual error, except for the incorrect oscillation in the predicted error. When S < 0.9, the
predicted error can be negative when it should be positive. This could be because many
grid spacing were larger than 0.2 m so that the cell Reynolds number were greater than 2.
Though the trend is generally captured by the error-transport equation with the residual
modeled by the modified equation, the predicted error magnitude can be over or under

predicted.

If we examine only the actual error results, we observe the following about the effect

57

of S: Since the total number of grid points is fixed, the maximum grid spacing increases
and the minimum grid spacing decreases as S decreases. Thus, when S is small (e.g.,
0.6), grid spacing is very small near x = 1, where the gradient is the steepest and very
large near x = 0 where the gradient is the least steep. Accordingly, when S is small, error
is low near x = l, but high further away from x = 1. As S is increased from 0.6 to 0.85,
the error decreases with the location of the error peak shifted towards x = 1. As S is
further increased from S = 0.85 to 0.95, the error starts to increase in the region near x =
1 because the grid spacing is no longer fine enough there. Clearly, for a given number of

points, there is an optimum S.

3.3.2 Wave Equation

This sub-section presents results obtained for the wave equation with the FD/F V
equation given by Eq. (3.14), the discrete-error-transport equation given by Eq. (3.15),

and the modeled residual given by Eq. (3.41).

Figure 3.8 shows the actual residual (denoted as RA) and the residual modeled by
using the modified equation (denoted as RM) as a function of x after 11 = 100 time steps
for the I = 101 grid and after n = 1000 time steps for the I = 1001 grid. Since the Courant
number Co = Ck/h was set at 0.8 for both grids, the results shown for the case with n =
100 time steps and I = 101 and the results shown for the case with n = 1000 time steps
and I = 1001 are for the same instant of time. Figure 3.8 shows the modeled residual

(RM) to match the actual residual (RA) nearly perfectly for most of the waveform except

58

at its two ends for the case of I = 1001 and n =1000. The match is worse with the coarser

mesh (1 = 101) though the number of time steps is five times less.

Figure 3.9 shows the grid-induced error predicted by using the discrete-error-
transport equation with the actual and the modeled residual. Also shown in the figure is
the actual error. From this figure, it can be seen that if the actual residual (RA) is used,
then the error predicted by the discrete-error-transport equation (denoted as E-RA) is
identical to the actual error (denoted as E-A). In fact, the agreement is to the last
significant digit. If the modeled residual (RM) is used, then the actual error (EA) and
the predicted error (E-RM) differ in magnitude. However, the location of the error is
tracked perfectly. This is because the waveform and the grid-induced error are both
advected by the same constant speed C. Similar to the results obtained for the advection-
diffusion equation, the error predicted by using the modeled residual improves as the

mesh is refined.

Figure 3.10 illustrates the solution and error evolution over time. Figure 3.10(a)
shows the numerical solution to deteriorate over time, and Figure 3.10(b) shows the error
to accumulate over time. Error predicted by using the actual residual matches perfectly
with the actual error over time. But, the error predicted by using the modeled residual

deteriorates over time.

3.3.3 Inviscid-Burgers Equation

59

This sub-section presents results obtained for the inviscid Burgers equation with the
FD/FV equation given by Eq. (3.19), the discrete-error-transport equation given by Eqs.
(3.22b) and (3.23b) for Linearization Method 1 and Eqs. (3.25b) and (3.26b) for

Linearization Method 2, and the modeled residual given by Eq. (3.42).

When the FD/FV equations are not linear, the difference operator must be linearized
to construct the discrete-error-transport equation. Since the linearized operator may
differ from the original FD/F V operator used to generate solutions, the following
relationship may no longer hold:

LD(h, k) U? = 0 (3.44)
That is, the solution U? computed by the nonlinear FD/FV equations may not satisfy Eq.
(3.44). Thus, Eq. (3.44) is re-written as

LD(h,k) U? = Rc (3.45)
where Rc is equal to zero if the FD/F V equation is linear and can be nonzero if the F D/F V

equation is nonlinear. This implies the actual residual is not equal to

R? = R1 = LD(h, k) U23, (3.46)
Instead, it is equal to
R? = LD(h, k) (U23, — U?) = R8 — Rc (3,47)

For Linearization Method 1, Rc does not equal to zero so that the actual residual is
given by Eq. (3.47). For Linearization Method 2, Rc does equal to zero so that the actual

residual is given by Eq. (3.46).

60

 

 

 

Figure 3.11 shows the actual and the modeled residuals, and Figure 3.12 shows the
actual error and the error predicted by using the actual and the modeled residuals for the
two discrete-error-transport equations obtained by the two linearization methods for the
case where the solution is still smooth (i.e., the weak discontinuous solution has not yet
developed). These figures show the discrete-error-transport equations produced by both
linearization methods to predict grid-induced errors perfectly (except for the last digit) if

the actual residual is used.

Results were also obtained to evaluate the two linearization methods when there is a
weak discontinuous solution in the flow domain. Figure 3.13 shows the values of the
modeled and the actual residuals, and Figure 3.14 shows the predicted and the actual
error. Interestingly, both the conservative and the non-conservative linearization methods
give perfect results if the actual residual is used. In hind sight, this result is expected
because of two reasons. The first is that the error propagation speed, Uc, is not computed
by the discrete-error-transport equation. Rather, it is taken from the solution of the
nonlinear FD/FV equations. Second, the difference between the nonlinear F D/F V

operator and the linearized operator are accounted for by Rc and R3.

Figure 3.11 and Figure 3.13 show the modeled residual to be quite different fi'om the
actual residual. This is in sharp contrast to the results obtained for the wave-equation test
problem with half the number of grid points and the same initial waveform. This
difference can be explained as follows. The first is that a simplified version of the

modified equation was used. The second is that for the inviscid Burger equation, the

61

initial wave form get compressed on fore part and elongated on the back part as it
propagates. Thus, though the grid spacing may initially be adequate, the effective grid
spacing becomes coarser and coarser as the wave propagates into the flow domain. With
the residual being modeled inaccurately, Figure 3.12 and Figure 3.14 show the error

predicted by the discrete-error—transport equation to be quite poor.

3.3.4 Steady-Burgers Equation

This sub-section presents results obtained for the steady Burgers equation with the
FD/FV equation given by Eq. (3.28), the discrete-error-transport equation given by Eqs.
(3.31b) and (3.32b) for Linearization Method 1 and Eqs. (3.34b) and (3.35b) for

Linearization Method 2, and the modeled residual given by Eq. (3.43).

Results were obtained for two different grids, I = 21 and I =51. Figure 3.15 shows
the actual and the modeled residuals, and Figure 3.16 shows the actual error and the
errors predicted by using the actual and the modeled residuals for the two discrete-error-
transport equations obtained by the two linearization methods. Similar to the inviscid-
Burgers-equation test problem, these figures show the discrete-error-transport equations
produced by both linearization methods to predict grid-induced errors perfectly (except
for the last digit) if the actual residual is used. Also similar to the inviscid-Burger-
equation test problem, the modeled residual is quite different from the actual residual.
This is in sharp contrast to the results obtained for the advection-diffusion—equation test

problem with the same Dirichlet boundary conditions. This difference can be explained

62

as follows. The first is that a simplified version of the modified equation was used. The
second is that the solution for the steady Burgers equation is quite different from the
advection-diffusion equation with much sharper gradients next to the boundary at x = I.
With the residual being modeled inaccurately, Figure 3.16 shows the error predicted by

the discrete-error-transport equation to be quite poor.

63

3.4. Summary of Chapter 3

Four l-D test problems — that account for linear, nonlinear, steady, and unsteady
flows - were used to evaluate the usefulness of the discrete-error-transport equations in
estimating grid-induced errors in finite-difference and finite-volume methods. The
residual in the discrete-error-transport equation was modeled by the leading term of the
remainder in the modified equation. For each test problem, solutions were generated for
a number of meshes, including one that is grid-independent. The grid-independent
solution for each test problem is used to assess the accuracy with which the residuals are
modeled by the modified equation and how well the discrete-error—transport equation

predicts grid-induced errors.

Results obtained show that the modified equation is useful in modeling the residual if
the mesh is sufficiently fine. Results obtained also show the discrete-error-transport
equation to be able to predict the error perfectly (up to the last significant digit) if the
actual residual is used. This is true even for the unsteady, nonlinear test problem with a
moving weak discontinuity in the flow domain. The linearization procedure used to
derive discrete-error-transport equations for nonlinear FD/F V equations — whether
conservative or non-conservative — was found to be unimportant because the error-

propagation speed is computed by the FD/FV equations, and not by the discrete-error

transport equation.

Table 3.1 Summary of Cases Studied for 1-D Advection—Diffusion Equation

 

 

S I Ax min (m) Ax max (m)
1.0 6 0.2 0.2
1.0 11 0.1 0.1
1.0 21 0.05 0.05
1.0 1013 0.01 0.01
1.0 10013 0.001 0.001
1.0 2001 0.0005 0.0005

0.95 21 0.0778 0.0294
0.95 41 7.76 x 10'3 0.057
0.9 21 0.0154 0.114
0.9 41 1.67 x 10'3 0.102
0.85 21 7.12 x 10'3 0.156
0.85 41 2.66 x 10“ 0.150
0.8 21 2.92 x 10'3 0.202
0.8 41 3.32 x 10'5 0.2

0.7 21 3.42 x 10“ 0.3

0.6 21 2.44 x 10'5 0.4

 

’Grid-independent solution.

65

0.4

0.3

0.2

0.1

0.0

0.4

0.3

0.2

0.1

0.0

 

0.5

 

 

 

0.0

RA

 

 

(a)1=11

66

 

 

I=11
S=0.95

 

I=11
S=0.95

 

100 —
1 |=41
50 _ S = 0.85

 

0.2 -

0.1 -

 

 

 

 

 

-0.1 —
-100 -
' -0.2 ‘ ‘
0 0 5 1 0 0 5 1
X X
RA
- ..... o ...... RM
100 —
|= 41
50 _ S=0.85

 

 

 

 

(b)I=41

Figure 3.1 Actual and modeled residual for advection-diffusion equation:

for several S and I

RA = actual residual, RM = modeled residual

x: physical domain, é: transformed domain.

67

I= 41
0-04 r S = 0.95

0.02 -

 

 

-0.02 —

 

 

0.02 -

0.00

 

 

-0.02

 

 

Figure 3.2 Actual and modeled residual for advection—diffirsion equation:
S = 0.95 and I =41, using I = 2001 as grid-independent solution,

RA = actual residual, RM = modeled residual

x: physical domain, .5: transformed domain.

68

 

 

 

1
1
.—
IS .v
. _ _ _ o
0 5. 0 6 0 M
o o
ANnc—‘Vw EE
m
A
.
c
.
m
.

 

 

E-RA

 

I=11

S

0.85

.\.
.\
.\
.\

 

0.5

 

(a)1=11

69

e (10“)
e (10“)

 

 

 

e (10“)

 

 

 

 

(b)I=41

Figure 3.3 Actual and predicted solution error (e) for advection-diffusion equation:
for several S and I.
E-A = actual error, E-RM = predicted error with RM, E—RA = predicted error with RA,

x: physical domain, 25,: transformed domain.

70

0.08 -

 

 

 

0 0.5 1
X X

Figure 3.4 Actual versus modeled residual for advection-diffusion equation: S = 1.

RA = actual residual, RM = modeled residual

0.04 =

O 0.02

 

 

 

 

 

0.00

 

Figure 3.5 Actual versus predicted error for advection-diffusion-equation: S = 1.
E-A = actual error, E-RM = predicted error with RM, E-RA = predicted error with RA,

71

 

 

 

 

(a)
0.10 — +6
0.05
0.00
-0.05 1 1
0.0 0.5 1.0
(b)

Figure 3.6 Error = F(grid spacing) for advection-diffusion equation:
S=1, I = 6,11, 21, and 101.
(a) Actual. (b) Predicted (Eqs. (2.15) and (3.9)).

72

0.015 3

0.010

0.005 .3

 

 

0.015 —

0.010 ~

0.005 _

 

0.000 "8% 3
0.0

 

(b)

73

0.015 1?

 

 

 

 

0.005 -
I M
-0I005 Tm" ”‘ ’ 1 ’ a
0.0 0.5 1.0
X
+ 0.6
(C)
—I— 0 7
+ 0.8 A
—X- 0.85
—u— 0.9
e + 0.95
1.0

 

((1)

Figure 3.7 Error = F (smoothness) for advection-diffusion equation:
I = 21 and S = 0.6, 0.7, 0.8, 0.85, 09, 0.95.

(a) Actual e = f(x). (b) Actual e = f(g). (c) Predicted e = f(x). ((1) Predicted e = f(é).

74

0.01

 

 

 

 

 

 

 

n: 0.00 ‘
l=101
Co=0£
'°-4 ‘ n=100 -0.01
6 8 r 6
X x

Figure 3.8 Actual versus modeled residual for wave equation.
RA = actual residual, RM = modeled residual

 

 

 

 

 

Figure 3.9 Actual versus predicted error for wave equation.
E-A = actual error, E-RM = predicted error with RM, E-RA = predicted error with RA

75

2.0 - n: 3°???

1.5 -

 

 

 

 

 

 

..4_l .

 

 

 

 

(b)

Figure 3.10 Results for wave equation: I=201, Co=0.8 at n = 100 and n = 200.

(a) Grid-independent and coarse-grid solution. Us) = grid-independent solution, U; =
solution from coarse mesh. (b) Solution error.

76

 

 

 

 

 

 

 

 

(b) Linearization Method 2

Figure 3.11 Actual versus modeled residual for inviscid Burgers equation:

I=201, 5 = 2.5, Um —k— :05, 11:40,
h h

RA = actual residual, RM = modeled residual, R8 = LD (h,k)U;i , Rc = LD (h,k) U?

77

-0.05

 

 

 

 

-0.05

 

 

 

(b) Linearization Method 2

Figure 3.12 Actual versus predicted error for inviscid Burgers equation:

I=201, 5 = 2.5, Um 5 =05, n=40
h h

E-A = actual error, E-RM = predicted error with RM, E-RA = predicted error with RA 3

78

0.0 MT.“

 

 

-0.1 ~

 

.—
lu-

X
(a) Linearization Method 1

 

'?

 

 

 

_
h

(b) Linearization Method 2

Figure 3.13 Actual versus modeled residual for inviscid Burgers equation:

with weak solution, I=201, E = 2.5, Umax E=0.5, n=400

RA = actual residual, RM = modeled residual, Rg = LD (h,k)U;i , Rc = LD (h,k) U?

79

 

 

 

 

 

 

 

a) -0.05 - .'
05:?
—— E-A
-0.10 - - ----- o ------ E-RM ii 3!
------- E-RA 65.;
2 4 6 8
x
(a) Linearization Method 1
0.01 -
0.00
d)
1
- L BA
0'01 - ..... o ...... E-RM
------- E-RA
2 4 6 8
X

(b) Linearization Method 2

Figure 3.14 Actual versus predicted error for inviscid Burgers equation:
with weak solution, I=201, % = 2.5, um % =05, n=400

E-A = actual error, E-RM = predicted error with RM, E-RA = predicted error with RA

80

 

 

 

 

 

 

(a) Linearization Method 1

 

 

 

 

 

 

(b) Linearization Method 2

Figure 3.15 Actual versus modeled residual for steady Burgers equation:
1= 21 and 51,

RA = actual residual, RM = modeled residual, R8 = LD (h,k) U23, , Rc = LD (h,k) U?

81

e (103)

 

 

 

 

(a). Linearization Method 1

 

 

 

0.03
4 -
A 0.02
1:3
° 05
1', 2 '
o 0.01
0 o.00<
0

 

 

(b) Linearization Method 2

Figure 3.16 Actual versus predicted error for steady Burgers equation:
I: 21 and 51,

E-A = actual error, E-RM = predicted error with RM, E-RA = predicted error with RA

82

CHAPTER 4
ANALYSIS AND MODELING OF RESIDUAL
IN THE DISCRETE-ERROR-TRANSPORT EQUATION

Since the modified equation and Taylor-series expansions are not usefirl in modeling
the residual when the grid is coarse (the typical situation), this chapter proposes an
alternative approach for residual modeling based on data mining. The approach involves
two steps. The first is to study the behavior of the residual by evaluating the “actual”
residual created by a variety of poor quality meshes in a systematic way. The second step
is to model the residual based on the understanding gained on the behavior of the

residual.

The remainder of this chapter is organized as follows. First, the test problems are
summarized and the discrete-error—transport equations derived. Afterwards, results of the
study on residual behavior and its modeling are described. The usefirlness of the DETE

in predicting the grid-induced errors is also demonstrated for the 2-D test problem.

4.1 Test Problems and Derivation of Discrete-Error-Transport Equations

The two objectives of this study are accomplished by using two test problems - the

one-dimensional (l-D) advection-diffusion equation and the two-dimensional (2-D)

advection-diffusion equation.

83

 

1-D Test Problem
The 1-D advection-diffusion equation and its discrete-error-transport equation are
discussed in Chapter 2 and Chapter 3.1.1, which are repeated below to make this chapter

more nearly self-contained.

The model equation is:

d d
LU=0, L=—— c- — 4.1
6‘ v6] ()

In the above equation, C, the advection speed, and v , the kinematic diffusivity, are
constants. In this study, C, v , and I! were set equal to 10, l, and 1, respectively, so that
Re =C [/v =10. The boundary conditions are U(xo)=1 and U(Z) = 2, at x = x0 = 0

and x = f = 1 respective.

By using central-differencing for the advection and diffusion terms, Eq. (4.1)

 

 

becomes
Ui+l —Ui _ U1 _Ui—|
CM=V x1+1 _xi xi “X14
x1+1 _xi-l xi+l/2 —xi—|/2 (4.2)

i= 2, 3,... 1, U1 =U(xo) = 1, U1=U(€)= 2
where U; = U(xi);.x1 is the location of the i’th grid point; X1+1/2 = (x1 + x111)/2; and I is the
total number of grid points. Equation (4.2) allows for variable grid spacing and can be
written in the following operator form:

LD (h) U1 = 0 (4.3a)

where

84

 

 

_c(E*'-E-')_ v 1 +_ __1_ --
LD— h. +113 hm‘h (131 13°) h(13° E')‘, (4.36)

1+1

h. =x. —xi_,, h.

1 1 1+1 =X —xi’ hi+ll2 =xi+ll2 -xi-l/2

i+1

In Eq. (4.3), E is the shift operator (i.e., Ein Ui = Um ).

Ifthe grid spacing h in Eq. (4.3) is refined to hg, where it is small enough to yield the
grid-independent solution U31 ; i.e., LD(hg) Ug’i = 0, then substitution of Us; into Eq. (4.3)
with a coarse mesh h (i.e., h > ha) will yield a residual; i.e.,

LD(h) Us, = R1 (4.4)
Subtracting Eq. (4.3) from Eq. (4.4) yields the discrete-error-transport equation for the
FD/FV equations given by Eq. (4.3), which is

' L,,(11)ei =R1 (4.5a)
where e, is given by

e; = U81 — U1 (4.5b)
e,- represents the grid-induced error since it is defined with respect to the grid-independent

solution and not the exact solution.

2-D Test Problem

The 2-D advection-diffusion equation studied is given by

LU=0, L =513_(C3 =1;-£9—‘4—‘2-‘Cy -V-—a-] (4.6)
X

In the above equation, Cx and Cy, the advection speed in x- and y-directions, and v , the

kinematic diffusivity, are constants. The boundary conditions at x = y = 0 and x = y = L

85

 

.f
1

= 1 are such that the solution is l-D along the vector direction '6, = i +j (i.e., pointing

45 degrees between the x and y coordinates) so that the exact solution is given by

exp(ReiE—J —1
WE) - U(éo) _ f

U(0—U(§.) ’ exp(Re)—1

(4.7a)

 

where § is along 6;; E, o and Z correspond to x = y= 0 and x = y = L = 1, respectively;

U(§6)=1;U(€)=2;and

Re=,/Ci +c‘j 8/1» (4.76)

The square domain of this 2-D problem is replaced by grid lines that are parallel to
the x- and y-coordinate lines. By using central-differencing for the advection and

diffusion terms on such a grid, Eq. (4.6) becomes

 

 

 

 

 

 

Cx Ui+l.j —Ui—l,j “I'Cy Ui,j+l “Um—1 =
xi+l —xi-l 3’34: "' Yj-r
4.8
Ui+l,j “Um __ Ui.j"Ui—1,j U614! _Ui.i _ Ui.j —Ui,j-l ( )
v x141 —xi xi —xi-l +v yj+l 'Yj Y) "Yj-r
xi+ll2 _xi-l/2 Yj+1/2 ”yjuz

where i = 2, 3,... I andj = 2, 3, ..., J; Um = U(xi, yj);.Xi and yj is the location ofthe (i, j)th
grid point; I is the total number of grid lines along x; and J is the total number of grid
lines along y. Equation (4.8) can be written in the following operator form:

LD (h, k) U1.) = 0 (4.93)

where

 

_C,(E;‘-E;')_ v 1 ,_ __1_ _ -
LD’ h. +h1 hi,,,,‘h (E: Ba) 11.032 En]

1+1 i+1 r

86

 

 

Cy(E:;l _E;') V 1 +1 0 1 o _| 4 9
+ km+kj _k,-+,,2 [k (By -Ey)-k—(E’_EY )‘ ( . b)

i+l i

h. =x. —xi_,, h

:x.

1+1

" xi ’ hi+ll2 = x1+1/2 —xi-l/2

1+1

1‘1 =y1"y1-w ki+l =Y1+1‘Y1-a k3,“, =Y1+112 'y1-112

In Eq. (4.9), E is the shift operator in the x and y directions (i.e., E?“ U11 = U and

iin.j

Ej“Ui3j = UM,“ ).

If the grid spacing h and k in Eq. (4.9a) is refined to hg and kg, where they are small
enough to yield the grid-independent solution U33) ; i.e., LD(hg, kg) U8),- = 0, then
substitution of U8),- into Eq. (4.9a) with a coarse mesh h and k will yield a residual; i.e.,

L001, k) U311 = Rid (4-10)
Subtracting Eq. (4.9a) from Eq. (4.10) yields the discrete-error-transport equation for the
FD/FV equations given by Eq.(4.9a), which is

LD(h, k) e?)- =Ri3j (4.11)
where e13,- is given by

61.1 = U311 - U1.1' (4-12)
c; J represents the grid-induced error since it is defined with respect to the grid-

independent solution and not the exact solution.

Though this test problem is very similar to the 1-D test problem, it enables grid aspect

ratio to be investigated in addition to lack of resolution and smoothness.

87

 

4.2. Behavior of the Residual in the Discrete-Error-Transport Equation

For each of the two test problems, a grid-independent solution is generated. For the
1-D test problem, the grid-independent solution was generated by using 1001 grid points.
When compared to the solution obtained by using I = 2001 grid points, the maximum
percent difference in the solutions is less than 0.005%. For the 2-D test problem, the
grid-independent solution was generated by using 401 x 401 grid points. When
compared to the solution obtained by using 801 x 801 grid points, the maximum percent

difference in the solutions is less than 0.0005%.

With the grid-independent solutions established, the actual residual for any poor

quality or coarse grid for the two test problems can be computed by using Eqs. (4.4) and

(4.10); i.e.,
R1=L6(h)U.,1 (4.13)
8.. =60. k) U...- (4.14)

Results from l-D Test Problem

For the 1-D test problem, this is what was done. Start with the grid that yielded the
grid—independent solution, choose a grid point at x,, remove n grid points following x =
x,,(see Figure 4.1 for schematic of removing two grid points) and then compute the
residual by using Eq.(4.13). The xr locations investigated are 0.2, 0.6, and 0.8. Each xr
location represents a different gradient of solution (U) with the gradient at xr = 0.2 being

the lowest and the gradient at xr = 0.8 being the highest. The number of grid points

88

removed is as follows: n = 1, 3, 6, and 9. By removing 11 grid points, both the resolution
and the grid smoothness are changed. A summary of the parameters studied is given in

Table 4.1.

By modifying a “perfect” grid at only one location in a controlled way, it is hoped to
understand the behavior of the residual as a function of lack of resolution and lack of
smoothness. The ideal resolution is the grid spacing corresponding to the grid-
independent solution (A xgj), and the ideal smoothness is unity. Thus, a resolution
parameter can be defined by

RDX1=DX1/Ax8,1, DXi =0.5(xi+l —xi_,) (4.15)
where RDX is equal to one if no grid points are removed, 1.5 when one grid point is
removed, 2 when two grid points are removed, and so on.

The smoothness parameter is defined by

Si = (X1+3/2 - Xi+1/2) / (Xi+l/2 - X140) (4-16)
S > 1 if the grid spacing gets larger along x, and S < 1 if the grid spacing gets smaller

along x.

Results of the residual calculated by using Eq. (4.13) for the imperfect grids
summarized in Table 4.1 are as follows:
1. The residual is essentially zero everywhere except at the grid point just before and
just after where grid points are removed. This is true whether the number of lines
removed is n = 1, 3, 6, or 9 and whether the location is xr = 0.2, 0.6, or 0.8. This

local nature of the residual can be seen in Figure 4.2.

89

2. When the grid spacing goes from small to large (S > 1), the residual is positive,
and when the grid spacing goes from large to small (S < 1), the residual is
negative (see Figure 4.2).

3. The absolute value of the residual increases as the number of grid points removed
is increased (see Figure 4.2). This is because the resolution decreases and lack of
smoothness increases.

4. The absolute value of the residual increases as xr increases since the gradient of U

increases (not shown).

Results from 2-D Test Problem

For the 2-D test problem, a series of different sequences of imperfect grids were
investigated. For each sequence, this is what was done. Start with the grid that yielded
the grid-independent solution, choose grid location(s) xr and/or yr, remove n grid lines

just after x = xr and/or y = y,, and then compute the residual by using Eq. (4.14).

For the first sequence of numerical experiments, only grid lines with constant x are
removed. The xr locations, after which they started to get removed are 0.25, 0.50, and
0.75. Similar to the 1-D test problem, each xr location represents a different gradient of
U with the gradient at xr = 0.25 being the lowest and the gradient at x, = 0.75 being the
highest. The number of grid lines removed is as follows: 11 = 1, 5, and 10. A summary of

the parameters studied for this sequence is given in Table 4.2.

Results of the residual calculated by using Eq. (4.14) for the imperfect grids

90

 

summarized in Table 4.2 are very similar to those found for the 1-D test problem. As
shown in Figure 4.3, the residual is essentially zero everywhere except along the grid line
just before and just after where grid lines are removed. Also, the residual is positive
when the grid spacing goes from small to large (S > 1), and negative when the grid
spacing goes from large to small (S < 1). Finally, the absolute value of the residual
increases as the number of grid lines removed is increased and as the gradient of U

increases.

For the second sequence of imperfect grids, the goal is to study the effects of
increasing'versus decreasing grid spacing (i.e., S > 1 versus S < 1) when the gradient is
the same (see Figure 4.4). The xr locations, where S is greater than or less than unity, are
x = 0.2775, 0.5050, and 0.7650. The number of grid lines removed at these locations is
10, 1, and 5, respectively. A summary of the parameters studied for this sequence is

given in Table 4.3.

Typical results of the residual calculated by using Eq. (4.14) for the imperfect grids
summarized in Table 4.3 are given in Figure 4.5. Basically, the absolute value of the
residual, |R|, is higher when S > 1 than when S < 1 for the same gradient of U. This is
expected since aliasing error occurs only when solution propagates from a fine mesh to a

coarse mesh, but not vice versa.

For the third sequence of imperfect grids, the goal is to study the effects of removing

grid lines in the x and in the y directions. A summary of the parameters studied for this

91

sequence is given in Table 4.4. Typical results obtained are shown in Figure 4.6. From
Figure 4.6, it can be seen that along x or y = constant lines where S > 1 or < 1, absolute
value of residual (|R|) increases with y or x because the gradient of U increases. It can
also be seen that |R| increases when grid lines removed in the x and y directions intersect,

but this effect is local.

For the fourth sequence of imperfect grids, the goal is to study the effects of
insufficient resolution. A summary of the parameters studied for this sequence is given
in Table 4.5. Note that all cells are square, and the smoothness parameter S = 1. Typical
results obtained are shown in Figure 4.7. From Figure 4.7, it can be seen that R increases

with grid spacing and with the magnitude of the gradient of solution (U).

For the fifth sequence of imperfect grids, the goal is to study the effects of aspect
ratio. A summary of the parameters studied for this sequence is given in Table 4.6.
Typical results obtained are shown in Figure 4.8. From Figure 4.8, it can be seen that

residual (R) increases with aspect ratio.

For the sixth sequence of imperfect grids, the goal is to study the effects of two
regions with poor smoothness in close proximity to each other. A summary of the
parameters studied for this sequence is given in Table 4.7. Typical results obtained are
shown in Figure 4.9. From Figure 4.9, it can be seen that the residual generated at each
non-smooth grid location is relatively unaffected by the residuals generated at upstream

or downstream locations.

92

4.3. Results of the Discrete-Error-Transport Equation in Estimating Errors

Since the usefulness of the DETE in estimating the error for the 1-D Advection-
Diffusion equation has already been demonstrated in Chapter 3, only the results for the 2~

D Advection-Diffusion equation will be presented in this section.

With the grid-independent solutions established, the actual error is calculated using
Eq. (4.12) and the error estimation is obtained using the DETE in Eq. (4.11) for the
imperfect grids with skipped grid lines. The actual residual calculated using Eq. (4.10) is

used as the error source for the DETE.

Figure 4.10 and Figure 4.11 show the absolute error for the case with 10 skipped grid
lines at xr = 0.75. Figure 4.10 shows the actual absolute error and Figure 4.11 shows'the
absolute error predicted using DETE with the actual residual used as the error source. It
can be seen that the predicted error matches perfectly with the actual error. This result is

consistent with the conclusions reached for the four l-D test problems in Chapter 3.

93

 

4.4. Modeling the Residual in the Discrete-Error-Transport Equation

In Section 4.2, the behavior of the residual was explored by performing numerical
experiments. Though the database generated is clearly incomplete, this section seeks to

model the residual by using this database.

l-D Test Problem

The smoothness factor used in residual modeling is defined as:

 

 

 

 

 

— DX. DX-
S=max( ‘ , ' ) (4.17a)
Dxi-r Dxi+1
DX. ‘
35(— ‘DX. , if DX,» I)x,_l
where, i =4 DXH ‘ (4.17b)
‘" ——“‘-, if DXH> DXi
‘DXi J

 

The database generated for the 1-D test problem was fitted by using least-square of
the following parameters: RDX for resolution given by Eq. (4.15), S for smoothness
given by Eq. (4.17a), and solution gradient dU/dx. Three firnctional relationships were

tried: linear polynomial, quadratic polynomial, and power law.

Assuming a linear polynomial in RDX, S, and dUi/dx = (U141 — U1-1)/(x.-+1 — x14) gives
the following modeled residual:

Rm = a(sign(S).(S—1))+ b(RDX — 1) + c €12 (4.18)

94

where a = 0.0179, b = 5.3413 x 105, = —3.3366 x 104. Also, sign(S)=1 if S>1, and

sign(S)= -1 if S<=1. Figure 4.12 shows that there is considerable spread with the fit

given by Eq. (4.18).

Assuming a quadratic polynomial gives

R... = a(sign(S)(§-1»+ b(RDx — 1) + c [Signfséfs - ”12

 

(4.19)

(RDX—1)2

+d +e(sign(S)(S—l)(RDX —1))

where a = 0.0162, b = 6.2235 x 10'6, c = —4.3552 x 104, d = 4.3552 x 104, and e =

3.7484 x 10". Figure 4.13 shows that Eq. (4.19) also has considerable spread.

Assuming a power law gives
Ran = const * sign(S) - (31):)“ - (S — 1)“ -(RDX)’ (4.20)

where const = 0.0334, 01 =1.0, B = 1.0012, 7 = —0.0027. Figure 4.14 shows that Eq.

(4.20) fits the data quite well. In fact, the fit is within 0.4%. Recall that if the residual
can be modeled perfectly, then the discrete error-transport equation will be able to predict

the grid-induced error perfectly.

2-D Test Problem
Similar to the 1-D test problem, the database generated for the 2-D test problem was
fitted by using least square minimization of the following parameters (see Figure 4.15):

RDS for resolution, S for smoothness, aU/as for the solution gradient along the flow

95

 

direction, and AR for the aspect ratio. These parameters are defined below, where the

over bar is defined in the same manner as that in the 1-D test problem.

 

 

 

RDS = (ISM/(18811.19 dsm. = \lef + ijz .cos(or —13) (4.21)

_ dsij dsij dsij dsij

S=max( ', ', ' , ’ ) (4.22)
dSi-l,j (15141.1 dsi.j-l dsi,j+l

av_ava_x 4112:

———— +
65 6x 6s 6y 6s (4'23)

AR = Ax/ Ay or Ay/ Ax , chosen so that AR 21 (4.24)

where or is the flow direction; 0 is the cell diagonal angle;

Since the power law fit reasonably well for the 1-D test problem, it was generalized to

the 2-D test problem. This gives
. all a - P 1
'RM = const “' s1gn(S) * (g) (S — 1) (RDS) (4.25)

But, the coefficients needed to be changed to produce best fit. For Eq. (4.25), the
coefficients are const=0.1186, or =0.9998, B =1.0039, and y = -0.0125. The goodness
of the fit is shown in Figs. 35 and 36. Since the error is nearly a constant for a given

smoothness. Equation (4.25) can readily be refitted by a correction factor; i.e.,

 

_const . ,@0-_ B , l
- f(S) *Slgn(S) (as) (S 1) (RDS) (4.26a)
where
f(S) =1-sign(S)*(2.51*S —1.29)*0.01 (4.26b)

With this correction, the fit given by Eq. (4.26) is accurate within 0.5%.

96

4.5 Summary of Chapter 4

The discrete error-transport equation (DETE) has the potential to predict grid-induced
errors in CFD solutions generated by poor quality or coarse grids. But, its usefulness
hinges on how well the residual in the DETE is modeled. This chapter provides
prelirrrinary results obtained by using data mining. First, the behavior of the residual is
analyzed for a variety of poor quality/coarse grids for two test problems: the 1-D and the
2-D advection-diffusion equations. The actual residuals generated were then curve-fit by
using least-square minimization. The power law was found to provide the best fit. But,
more studies are needed before a general model for the residual can be constructed. The
usefulness of the DETE in predicting the grid-induced is also evaluated for the 2-D
Advection-Diffusion equation. Results obtained show that the discrete-error-transport
equation is able to predict the error perfectly if the actual residual is used. This result is
consistent with the conclusions reached for the four l-D test problems studied in Chapter

3.

97

 

Xr

Figure 4.1 Schematic of removing two grid points after x = xr for the 1—D test problem.

Table 4.1 Summary of Parameters Investigated for the 1-D Test Problem

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Location Number of Grid Points
Removed
x,=0.2 n=1,3,6,9
xr=0.6 n=1,3,6,9
xr=0.8 n=1,3,6,9
l
0.05 — 0.1 '3
m 0_00 V "7T. ' ..TI T. .T. m O_O . .. T :.. . 7T:_;:._ .: .
' l
L
-o.05 ~ -0-1 3-
0.75 ‘ 0.80 0.85 0.75 0.80 0.85
X X
(a) (b)

Figure 4.2 Residual computed by using Eq. (4.13) for the l-D test problem.
(a): xr = 0.8 and n = 3.; (b): Bottom: xr = 0.8 and n = 6.

98

Table 4.2 Summary of Parameters Investigated in Sequence 1 for the 2-D Test Problem

 

 

 

 

 

 

 

 

Location Number of Grid Lines
Removed
x.=0.25 n=l,5,10
x,=0.50 n=1,5,10
x,=0.75 n=1,5,10
1.00

RA (10*)
3.42

  

2.93
2.44
1.95
1.47
0.98
.4
> 0.95 ° 9

0.90

0.75

 

X

Figure 4.3 Residual computed by using Eq. (4.14) for the 2-D test problem with
Sequence 1 (Table 4.2).

Top: xr = 0.75 and n = 1; Bottom: xr = 0.75 and n = 10.

99

Table 4.3 Summary of Parameters Investigated in Sequence 2 for the 2-D Test Problem

 

 

 

 

Location Number of Grid Lines
Removed
xr = 0.2775 n = 10
xr = 0.5050 n = 1
xr = 0.7650 n = 5

 

 

 

 

 

 

 

S>1

 

 

Figure 4.4 Schematic showing where grid lines are removed to study the effects of S > 1
and S < 1 for the 2-D test problem.

 

 

 

3.0 .
‘1?
c 2.0 -
‘-
E 1.0»

0.0 ' 4

o 0.5 1
y

Figure 4.5 Residual computed by using Eq. (4.14) for the 2-D test problem with
Sequence 2 (Table 4.3):

x = 0.2775, n =10.

100

Table 4.4 Summary of Parameters Investigated in Sequence 3 for the 2-D Test Problem

 

 

 

 

 

 

 

 

 

 

 

 

Location Number of Grid Lines
Removed

x,=0.25, y,= 0.25 n= 1, 5,10

xr = 0.75, y,= 0.75 n= 1, 5

x,=0.50,y,=0.25 n=1

x,=0.50,y,=0.50 n=1

x,=0.50,y,=0.75 n=l
:l‘ SLR-J i:tit‘“t‘f:€’::i:fj;;i::i’
‘:'171*"1‘ ’7’”? ...;’:’3,:’f”
iiiéflii” 15431 ’ “1“1’ ZEEE3§%:i?
;:t. T 1M1 7771;":“‘:1::lj¢':— JLlPl-XY
..A; . 1.111' fl: 1 1 3';' 1
7..- t';:1;.:;.;:::; 1 1 ‘ [11
H9? » Hi 17: 7‘, t‘l‘ , 1:;l;;1;;;—JLl-XY
- --—— ; 351%:eeitjiiii
1‘19 """" l ,1; , i l::|i ; . I11ittl ‘

ILl-X ILlPl-X

ILl-XY ILlPl-XY

(a) schematic showing the locations where data are plotted

- - e- - JL1P1-XY
0.06 - lL1-X
----- IL1P1-X

 

.r"

 

 
 
   

K 0.00

 

 

-—d-— |L1 -XY
- - O - - lL1P1-XY
".1 -X

----- IL1 P1 -X

 
 
 

0.06

 

 

 

......

-0.06

 

 

0.7

(b) actual residual at locations with skipped grid lines
Figure 4.6 Residual computed by using Eq. (4.14) for the 2-D test problem with

Sequence 3 (Table 4.4):
x, = 0.75, yr = 0.75 and n = 5.

101

Table 4.5 Summary of Parameters Investigated in Sequence 4 for the 2-D Test Problem

 

Number of Grid Points (1 x J)
401x 401
201 x 201
101 x101
21 x 21
11x11

 

 

 

 

 

 

 

 

4.0 -—o—— 201x201

1- ..... . ------ 101x101
---,--- 21x21

«27‘ 3'°'-—-—--<-——. 11x11 I

 

 

Figure 4.7 Residual computed by using Eq. (4.14) for the 2-D test problem with
Sequence 4 (Table 4.5):

as a function of y along x = 0.5.

Table 4.6 Summary of Parameters Investigated in Sequence 5 for the 2-D Test Problem

 

Number of Grid Points (1 x J)
201x 201 201x101 101x 101
101x101 101x 21 21 x 21
101x101 101x11 11x11

 

 

 

 

 

 

 

 

102

1“; '.' ."

 

6.0

—o-— 201 x201
101x101
201x101

.. ..... . ......

 

——o—— 101 X101
- ..... g ...... 21x21

101x21

p—

0.9

---'---

0.6

0.3

 

0.01

Figure 4.8 Residual computed by using Eq. (4.14) for the 2-D test problem with

 

O
I
o
a
a

1'

Sequence 5 (Table 4.6):

as a function of y along x = 0.9.

103

Table 4.7 Summary of Parameters Investigated in Sequence 6 for the 2-D Test Problem

 

 

 

 

 

 

Location Number of Grid Lines
Removed after Each xi
x1=0.25,xz=0.50 n=1,5
x1=0.25,xz=0.75 n= 1,5
x1 = 0.70, xz = 0.75 n = 5
x1 = 0.7250, xz = 0.75 n = 5

 

 

 

 

Figure 4.9 Residual computed by using Eq. (4.14) for the 2-D test problem with
Sequence 6 (Table 4.7).

(104)
4.71
4.1 5
3.58
3.02
2.48
1.89
1.33
0.77
0.21

-0.36

-0.92

-1.48

-2.05

-2.61

-3.17

 

0 .6 0 .8 1 .0
X
Figure 4.10 Actual absolute error calculated using Eq. (4.12) for the 2-D test problem:
xr = 0.75, n =10

(1 0'51
4.71
4.1 5
3.58
3.02
2.46
1 .89
1 .33
0.77
0.21

-0.36

-0.92

-1 .48

-2.05

-2.61

-3.1 7

 

0.6 0 .8 1 .0
X
Figure 4.11 Predicted absolute error using Eq. (4.11), with actual residual calculated
using Eq. (4.10) for the 2-D test problem:
xr = 0.75, n =10

105

0.4

 

E?!

0.2}r

 

R
9
+++6H+
W
OOWOO

 

 

0 0.5 1
x

Figure 4.12 Actual residual (RA) and residual modeled (RM) by Eq. (4.18) for the 1-D
test problem.

 

R
O
+++6+++
+W++
009143-9300

 

' 0 015 1
x

Figure 4.13 Actual residual (RA) and residual modeled (RM) by Eq. (4.19) for the 1-D
test problem.

106

 

J 12X 1
1 + RM

0.2L ED 1

1 6
1 e ‘
I: 01 c Q g 1
69 1
-0.2i : 1
'0'40 05 I

Figure 4.14 Actual residual (RA) and residual modeled (RM) by Eq. (4.20) for the 1-D
test problem.

 

Figure 4.15 Definition of on and B for Eq. (4.21) for the 2-D test problem.

107

 

 

 

0 RA
1 + RM
0.21
1
"1
m 1
-0.2‘
-0.41 ‘
‘ 1
'0'60 0.5 1

Figure 4.16 Actual residual (RA) and residual modeled (RM) by Eq. (4.25) for the 2-D
test problem.

 

 

w ’TA’ m n
’~ kooooocooooooooooooo o 051
o ' ’
E; + 0255
_ 5000000000000000000001 6 03A0
0 1 . 0501
t gogooofiovooooooooooai fl 0M5
I.” 0‘ 0 0.50.10
a: g > 0.75.1
> 05.9009990999999099‘ 4 and
2; O ' 1' 0.75,10
m dQlQQIGQQQQQQQIQOQQQ
a '5’ 1
M .oooocoooooooooocooo

.m-gfi—— __...‘. 4‘7«—~1

0 05 1
Y

Figure 4.17 Error in predicted residual: (RA-RM)/RA for the 2-D test problem.

108

CHAPTER 5

ERROR ESTIMATION FOR PROBLEMS GOVERNED BY
EULER EQUATIONS

So far, the discrete-error transport equation (DETE) has only been developed and
demonstrated for simple model equations such as the advection diffusion equations and
the inviscid Burger equation. Also, it has only been applied on structured grids, where
finite-difference and finite-volume methods can lead to identical algebraic analogs of
PDES. Thus, the objective of this chapter is to develop a DETE for a more complicated

set of equations in the framework of arbitrary elements.

The remainder of this chapter is organized as follows. First, the set of test equations
selected is summarized. Next, the finite-volume method (FVM) enabling solution on
arbitrary elements is derived. Then, the DETE corresponding to this FVM is derived.
This is followed by the description of the test problems and finally the usefiilness of

DETE for the set of selected equations are examined.

5.1 Description of Governing Equations

The set of governing equations (PDEs) selected is the continuity, x- and y-
momentum, and total-energy equations valid for two-dimensional compressible, inviscid

flow of a colorically and thermally perfect gas. These equations can be written as

109

10.21: 66-

 

 

 

 

 

 

+——- 0 ,
at 6x 6y (51a)
rp \ rpu \ rpv \
2 uv
Q=<pu1, F=<pu ”3 >. G=1p , > (5.1b)
pv puv pv +p
.13 J Lu(E+p)J \v(E+p)J
E =—£—+lp(u2 +v2) (5.1c)
'11—] 2

where p is density; u and v are the x- and y-components of the velocity, respectively; B
is mechanical and thermal energy; P is pressure; T is temperature; and y is the ratio of

specific heats.

110

5.2 Finite-Volume Method of Solution

In this section, a cell-centered finite-volume method is used to generate solutions to
Eq. (5.1) in which the mesh can be structured or unstructured and the cells can have
arbitrary shapes such as triangles, rectangles, or other polygons. The method proceeds as

follows (Wang (2002)).

Integrating Eq. (5.1) over an arbitrary i’th cell with volume V1 in a flow domain

replaced by a structured or unstructured mesh gives

f(ig+_a_g+a_6)dv=0 (5.2a)
V. at 6x By
01'
[[23+VoinV=0 (5.2b)
v. at
i=Fi+Gj (5.2c)

where i and E are unit vectors in the x- and y-directions. By invoking Gauss’ theorem,

[Voidv= F-fidS,fi=n.i+n,i (53)
6V

Vi

Eq. (5.2) becomes (assmning cell does not change in time)

6Q. _.
——'V. + f on d8 = 0
at ‘ WI: (5.4a)

— 1
Qi ..—v— [Q dV (5.4b)

ivi

111

If the boundary 6V, that surrounds cell volume V; involves J1 faces (with the faces

being lines in 2-D and planes in 3-D), then Eq. (5.4) becomes

a“.
—&=‘Pi (5.5a)
at
‘1’ “Li [fonds—ml—i ' 3 (55b)
I Vi j=l Sid \Ii j=| I.) I.) .
f 1 f "
i,j = -S—- Onds (5.50)
I.) Sid

where Sij denotes the length of the jth face.

So far, no approximations have been made. In this study, the mean fluxes on the
faces given by Eq. (5.5c) is approximated by the following first-order Lax-Friedrichs

formula, which is TVD preserving,

f... = 1(6.,6.,H)=1 «6m £16m —( H.016. —6.)1/2 (5.6)

 

 

VILIV
In the above equation, the subscripts L and R refer to the left and right cells about face j;

v”, =(vn’L +vn,R)/2 is the average face normal velocity; and c,v =(cL +cR)/ 2 is the

average speed of sound.

The time derivative in Eq. (5.5a) is approximated by the following two-stage Runge-

Kutta scheme:
61“) = 61“ + At ‘11“??? (5.7a)
6"“ =16" +161" +1211 945‘”) (5.71»)
l 2 I 2 I 2 I I

112

where the solution (6" ) at the 11th time level (t") is assumed to be known, and the solution

(6”) at the (n+1)th time level (t"”) is sought.

Since only steady-state solutions are of interest in this study, local time stepping is

implemented to accelerate convergence rate. The time step size at each cell (Ari) is

limited by CFL condition given by

At, = CFL Vi/i(|vn|+ c)j s”. (5.7c)

i=1

This completes the finite-volume method employed if only the solution is sought.
However, if we also seek to estimate grid-induced errors, then it is important to note that

F and G in Eq. (5.2c) can also be written as

F =AQ , G = BQ (5.8a)
where
_ 6F =
6Q
r 0 1 0 o \
L302 +1—11-v2 (3-y)u (1—y)v y-l
2 2 (5 8b)
—uv v u 0 '
(7’1)(u2+vz)u"YB—E' '12—?(3112-I-v2)+-'-Y-E (l—y)uv yu
(3 p

 

K

”I?
BIS
11

113

 

 

 

 

f 0 0 1 o)
-uv v u 0
1§3v2+7—;—1u2 (1—1)u (3—1)v 7-1 (53°)
(y—l)(u2+v2)v——Y-l£ (l—y)uv l:—Y(3v2+u2)+E yv
\ P 2 P j
2
1‘3: ° +1(u2+v2) (5.8d)
9 7-1 2

This is because F and G are homogeneous functions of Q to degree one. Thus, there
are two ways to evaluate F and G, via Eq. (5.1b) or via Eq. (5.8). Both give identical
results. With Eq. (5.8), f(Q) in Eq. (5.6) becomes

f(Q)=Aan+BQny (59)

Equations (5 .8) and (5.9) are needed because the DETEs require linearization.

Incidentally, these equations are also needed if an implicit method was used to

approximate the time derivatives.

114

5.3 Discrete-Error-Transport Equations (DETEs)

The derivation of the discrete-error-transport equations (DETEs) for the finite-
volume equations (FVEs) given by Eqs. (5.5) to (5.9) proceeds as follows. The first step
is to linearize the FVEs. As noted by Qin & Shih (2002), it is important to linearize by
using values of variables computed on the “coarse” mesh for which error estimation is

sought. Thus, the linearized FVEs have the following form:

in'I' 1 :l[fL+fR-(

dt (v) . 2 +0.1).(6a-6.>1(s.,).=o

 

 

 

VILBV

ft = ((A..5)L.n, +(B,.6)L.ny) (5.10)

r, =((A,.6),..nx +(B,.6)R.ny)
The time derivative is approximated by Eq. (5.7). In Eq. (5.10), all variables with

subscript c are evaluated by using the information on the coarse mesh.

It is important to note that though Eq. (5.10) is said to have been linearized, it is
identical to the non-linearized FVEs because of the identities given by Eq. (5.8). Thus,
when the solution Qc obtained on the coarse mesh is inserted into Eq. (5.10), they will
satisfy Eq. (5.10) perfectly with no residuals. Equation (5.10) is considered linearized
only when the grid-independent solution is inserted into it because in that case the
coefficients A, B, and V are not evaluated by using the information on the mesh used to
generate the grid-independent solution. Thus, when the grid independent solution Qg is

inserted into Eq. (5.10), a residual R8 will be produced.

115

Thus, second step is to subtract the Eq. (5.10) with Qc inserted from the Eq. (5.10)

with the Q1; inserted. The resulting equation constitutes a system of DETEs for the FVEs

given by Eqs. (5.5) to (5.9), which is

de. l J'1 , c
'+ X-z-[fL+fR—(

 

 

 

— V +c e —e S. . = R
dt (V1). ..1 ..).( R .m 1.1% .
° 5.11
fL =((Ac e)L'nx +(Bc e)L.ny) ( a)
f; =((:‘\c (2:)...nx +(Bc e)R.ny)
where e is the grid-induced error and is given by
r.p_8 -i5c ,
— _ (on). «33).
C=Q‘—Q° zi ’ (5.11b)

('p‘v). «5).
(E), -(E).

 

 

116

5.4 Test Problems

5.4.1 Test Problem 1: Oblique Shock Problem

The first test problem selected is uniform supersonic flow with freestream Mach
number M1 and pressure P1 over a sharp edge wedge with half angle 0. A schematic
diagram is shown in Figure 5.1. As can be seen in the figure, only half of the wedge is

simulated because of symmetry. For this problem, there is an exact solution given by

 

1 7+1 M?
———= -l tan0 .
tana K 2 )(Mfsinzej ] (512a)
32— = (A) Mf sin2 e {11) (5.1215)
R 7+1 7+1

pl_tan0= (7+1)Mfsin20
pl tana 2+(7—1)Mlzsin20

|V2| _ sin9 [ 2 {11)} (5.12d)

|Vl| — sin(0-01) (7+1)M,2 sin2 0 7+1

 

(5.120)

 

 

 

where a is the angle of the oblique shock with respect to upstream flow direction; 1

denotes the condition upstream of the oblique shock; 2 denotes the condition down
stream of the oblique shock; V1 = M1 1/7RT, (R is the gas constant); V2 is the magnitude

of the velocity downstream of the oblique shock. This solution is valid as long as the

angle of the wedge 0 is sufficiently small so that a bow shock does not form upstream of

the wedge.

The boundary conditions used are specified freestream conditions at the inflow

117

boundary, extrapolated boundary condition at the outflow boundary, zero flux/flow

boundary condition at the symmetry plane and the wedge surface.

5.4.2 Test Problem 2: Inviscid Flow around Airfoil

The problem of interest is the inviscid flow field around a 944 clean airfoil. The
computational domain and the mesh are shown in Figure 5.2. The boundary conditions
used are specified freestream conditions at the inflow boundary and the outflow

boundary, zero flux/ flow boundary condition at the airfoil surface.

118

5.5 Treatment of Boundary Conditions and Convergence Criteria

The treatment of the freestream and zero flux/flow boundary conditions for the
oblique shock and inviscid airfoil problems and the convergence criteria are presented in

this section.

5.5.1 Zero flux/flow boundary conditions

The zero flux/flow boundary conditions are specified as:

P: = p, (5.13a)
p3 = pi (5.1315)
V, = V. — 261111)}; (5.13c)

where i is the interior cell and g is the corresponding ghost cell.

5.5.2 Free stream boundary conditions

The free stream boundary conditions are specified via the characteristic variables:

_ p
W. -57 (5.14a)
W2 :1,“ +31 (5.14b)
7-1
w, =vn ——2i (5.14c)

119

Figure 5.3 gives a sketch of the characteristic freestream boundary conditions, where

g is the ghost cell, i is the interior cell and 00 is the boundary condition.

 

 

If vn<0,
wm=£%=wm=£% (5mm
9, 9..
2e 2e
“’2: =(—vng)+ :1: W200 =(—vm)+ 7.” (5.146)
2e 2e.
w =—v -—§—=w.=—v. ——'— 5.14
33 ( us) 7"] 31 ( m) 'Y-l ( f)

where, p8, p8 , v"g and c8 are the flow variables of the ghost cell, pa, , pa, , vn and

«3

ca, are the free stream conditions, and pi , pi, v"i and ci are the flow variables of the

interior cell. Then the flow variables of the ghost cell (for v,1 <0), p8, p8 , v"g and c8 can

be solved in terms of pm, pi, P... , pi , v vni , cm and ci , using equations (5.14d)-(5.14

nao’

 

 

f).
If vn >0,
P 1
w... =—§-,—=w. =3; (5.141;)
pg pi
2c .
WZg =vng+—_8_=w2i =vni+£ (514h)
7-1 7-1
20 2c0° .
W38 2v“lg -7—81 = w,“o = vmo — _1 (5.141)

120

Then the flow variables of the ghost cell (for vn <0), pg , pg , vng and c8 can be solved in

terms of P... , pi , pa, , pi , vm , V cm and ci , using equations (5.14g)-(5.l4 i).

ni’

5.5.3 Convergence Criteria

All the cases in this chapter are given a maximum number of steps to run and are
considered to converge when the residual can no longer be lowered and the residual has
already dropped at least 5 to 6 orders, relative to its initial value. In evaluating the

convergence, the L1 , L2 and Lu, norm of the residual have been used and are defined as

 

 

 

 

follows:
11 cells
4 Z (1Resjjl-vi)
L =max i=1 (5.158)
‘ jg) TVOl
11 cells
4 2(Resii -v,)
L = i“ , 5.15b
2 ”13'" TVol ( )
j=1.4 '
L0, = rnaleesj'i , (5.15c)
i=l,n_cells

where, Vi is the cell volume, W0] is the total volume of the computational domain, i
refers to the cell number and j refers to the four unknown variables (p , pu , pv and E ).

As shown in Eq. (5.15a)b) and (5.15a)c), the maximum L2 norm of the four unknown

121

variables is chosen as the L 2 norm, and the maximum absolute value of the residual of

the four unknown variables among all the cells is chosen as the L,o norm.

122

5.6 Results

With the Euler equations summarized, the finite-volume method of solution
described, and the discrete-error-transport equations (DETEs) derived, we now examine
the usefulness of the DETE concept for more complex PDEs, which are coupled and

quasilinear.

5.6.1 Oblique Shock Problem.

For the oblique shock problem, the grid independent solution was generated by using
641 x 321 grid points. Since the oblique shock wave is a weak solution involving a sharp
discontinuity, the solution generated on this or any other grid will clearly not be grid
independent, but simply one selected to be sufficiently good. Figure 5.4 shows the
solution in density for the 641 x 321 grid. The solution obtained on this grid is used to
compute the residuals (R8) in the DETEs (Eq.(5.11)) and to calculate the actual absolute
and relative errors, which can be used to assess how well DETEs predict errors on coarse

meshes.

Figure 5.5 to Figure 5.8 show the residual (R8) for the solutions generated on the
following four coarser meshes: 321 x 161, 161 x 81, 81 x 41, and 41 x 21. Only the
residuals for the continuity equation (density) are shown. These residuals were computed
by substituting the grid-independent solution into Eq. (5.10) with the coefficients A, B,

and V computed by using the information on the coarse meshes. Thus, the residuals

123

shown are exact. As expected, these figures show the residual to be high just upstream
and downstream of the oblique shock. The blobs and jaggedness in the figures are due to

the interpolation of the scheme used to plot contours.

Figure 5.9 to Figure 5.12 show the errors in the solutions generated on the following
four coarser meshes: 321 x 161, 161 x 81, 81 x 41, and 41x 21. Only the absolute errors
in the computed solution for density are shown. These errors were computed by using
Eq. (5.11b); that is by subtracting the coarse mesh solution from the grid-independent

solution at the center of each cell in the coarse mesh (e.g., e=pg —pc ). Bilinear

interpolation is used to transfer the grid-independent solution to the center of each cell in
the coarse mesh. Similar to the behavior of the residual, the error is high just upstream
and downstream of the oblique shock. As expected, the error increases as the mesh gets

COflISCI'.

To demonstrate the usefirlness of the DETEs derived, Eq. (5.11) with the exact R8
(e.g., those shown in Figure 5.5 to Figure 5.8) was used to compute the grid-induced
error. These computations show that the DETEs were able to predict perfectly the errors
if the R8 is exact. Thus, this result is similar to what was found when testing the DETE

concept on simpler model equations.

5.6.2 Airfoil Problem.

The grid independent solution was generated by using 769 x 129 grid points for the

124

inviscid airfoil problem. Figure 5.13 shows the solution in density for the 769 x 129 grid.
The solution obtained on this grid is used to compute the residuals (R3) in the DETEs
(Eq. (5.11)) and to calculate the actual absolute and relative errors, which can be used to

assess how well DETEs predict errors on coarse meshes.

Figure 5.14 to Figure 5.16 show the residual (R8) for the solutions generated on the
following three coarser meshes: 385 x 129, 193 x 65, and 97 x 33. Only the residuals for
the continuity equation (density) are shown. These residuals were computed by
substituting the grid-independent solution into Eq. (5.10) with the coefficients A, B, and
V computed by using the information on the coarse meshes. Thus, the residuals shown
are exact. As expected, these figures show the residual increases as the number of grid

cells decreases.

Figure 5.17 to Figure 5.19 show the errors in the solutions generated on the following
four coarser meshes: 385 x 129, 193 x 65, and 97 x 33. Only the absolute errors in the
computed solution for density are shown. These errors were computed by using Eq.
(5.11b); that is by subtracting the coarse mesh solution from the grid-independent

solution at the center of each cell in the coarse mesh (e.g., e: pg --pc ). Bilinear

interpolation is used to transfer the grid-independent solution to the center of each cell in

the coarse mesh. As expected, the error increases as the mesh gets coarser.

To demonstrate the usefulness of the DETEs derived, Eq. (5.11) with the exact R3

(e.g., those shown in Figure 5.14 to Figure 5.16) was used to compute the grid-induced

125

 

error. These computations show that the DETEs were able to predict perfectly the errors
if the R3 is exact. Thus, this result is similar to what was found when testing the DETE

concept on simpler model equations.

126

5.7 Summary of Chapter 5

In this chapter, the discrete error-transport equation (DETE) is derived for the two-
dimensional, compressible, Euler equations solved by a cell-centered finite-volume
method. The usefulness of the Euler DETE in predicting error is examined on two test
problems: the oblique shock problem and the inviscid flow around an airfoil problem.
The results show that the DETEs are able to predict perfectly the errors if the R8 is exact.
Thus, this result is similar to what was found when testing the DETE concept on simpler
model equations. It shows that DETE concept to be applicable to coupled, quasi-linear

PDEs with weak solutions.

127

Outflow

 

Freestream W311

 

 

 

Symmetry

Figure 5.1 Schematic of Oblique Shock Problem.

   
 

    

 

outflow 10 .
\1 51
>- 0
El -5
-10
_15 c ...5...
inflow -15 -10 -5 0 5 10 15

(a) computational domain and mesh

128

 

(b) mesh around the airfoil

Figure 5.2 Computational domain and local mesh for the airfoil problem.

(-vn)

(-vn)-c {—vn)+c

    

(a) v. <0 (b) v. >0

Figure 5.3 Schematic of freestream boundary conditions

129

2.364
2.230
2.097
1.964
1.830
1.697
~ -_ 1.563
1.430
1.296
1.163

>- 0.5

 

0.0

 

Figure 5.4 Grid-independent solution in density for the oblique shock problem:
641 x 321 case

1 .0
4.47E+03

3.42E+03
2.37E+03
1.32E+03
2.63E+02
-7.89E+02
-1.84E+03
-2.89E+03
-3.95E+03
-5.00E+03

 

0.0 . L .
1.0 1.5 2.0

X

Figure 5.5 Exact Residual R3 in the continuity equation for the oblique shock problem:
321 x 161case

130

8.95E+03
6.84E+03
4.74E+03
2.63E+03
5.26E+02
-1 .58E+03
-3.68E+03
-5.79E+03
-7.89E+03
-1 .00E+O4

 

 

0.0 . l .
1.0 1.5 2.0

X

Figure 5.6 Exact Residual Rs in the continuity equation for the oblique shock problem:
161 x 81 case

1 .0
1.72904
1.33E+04
9.47E+03
5.62E+03
1.76E+03
-2.09E+03
>. 0 _5 -5.9515+03
-9.80E+03
-1.37E+04
-1.75E+04
0 .0 - . . 1
1 .0 1 .5 2 .0

X

Figure 5.7 Exact Residual Rg in the continuity equation for the oblique shock problem:
81 x 41case

131

1.80E+04
1.40E+04
1.00E+04
6.00E+03
2.00E+03
-2.00E+03
-6.00E+03
-1.00E+04
-1.40E+04
-1.80E+04

 

 

Figure 5.8 Exact Residual Rg in the continuity equation for the oblique shock problem:
' 41 x 21 case

1.0

>'0.5

    

0.0

 

1.0 l 1.5 A 2.0
x

Figure 5.9 Actual absolute error in solution of density for the oblique shock problem:
321 x 161 case.

132

1.0

>'0.5

0.0

 

 

Figure 5.10 Actual absolute error in solution of density for the oblique shock problem:
161 x 81 case

>'0.5

    

0.0

 

1.0 ' 1.5 A 2.0
x

Figure 5.11 Actual absolute error in solution of density for the oblique shock problem:
81 x 41 case

133

I I I l I I
00000000000
0
(II
D

 

 

Figure 5.12 Actual absolute error in solution of density for the oblique shock problem:
41 x 21 case

 

 

X

Figure 5.13 Grid-independent solution in density for the inviscid airfoil problem:
769 x 129 case

134

 

—1
-1 0 1

X

Figure 5.14 Exact Residual Rg in the continuity equation for the inviscid airfoil problem:
385 x 129 case

   

-1 0 1
X

Figure 5.15 Exact Residual R1; in the continuity equation for the inviscid airfoil problem:
193 x 65 case

135

s
a
lil
E1
-
1!

   

-1 0 1
X

Figure 5.16 Exact Residual R3 in the continuity equation for the inviscid airfoil problem:
97 x 33 case

2.0E-03
1.4E-03
8.4E-04
2.6E-04
-3.2E-04
-6.SE-04
-1.5E-03
-2.1E-03
-2.6E-03
-3.2E-03
-3.6E-03
-4.4E-03
-4.9E-03
-5.5E-03
-6.1E-03
-6.7E-03
-7.3E—03
-7.8E-03
-8.4E-03
-9.0E-03

 

X

Figure 5.17 Actual absolute error in solution of density for the inviscid airfoil problem:
385 x 129 case

136

5.5E-03

3.6E-03

1 .8E-03
-1.1 E-04
-2.0E-03
-3.8E-03
-5.7E-03
-7.6E-03
-9.4E-03
-1.1E-02
-1 .3E-02
-1 .5E-02
-1 .7E-02
-1 .9E-02
-2.1E-02
-2.3E-02
-2.4E-02
-2.6E-02
-2.8E-02
-3.0E-02

 

-1
-1 O 1
X

Figure 5.18 Actual absolute error in solution of density for the inviscid airfoil problem:
193 x 65 case

uniigiiflmra.

 

X

Figure 5.19 Actual absolute error in solution of density for the inviscid airfoil problem:
97 x 33 case

137

CHAPTER 6

ERROR ESTIMATION FOR PROBLEMS GOVERNED BY
NAVIER-STOKES EQUATIONS

The discrete-error transport equation (DETE) has been developed and demonstrated
for simple model equations and the set of Euler equations. This chapter examines the
usefulness of the Euler DETE developed in Chapter 5 in estimating the grid-induced

errors in the solutions obtained by Navier-Stokes solvers.

6.1 Description of 2-D Navier-Stokes Equations

The set of Navier-Stokes governing equations (PDEs) chosen include the continuity,
x- and y-momentum, and total-energy equations valid for two-dimensional compressible,
viscous flow of a colorically and thermally perfect gas. These equations can be written as

6Q+6E +661 _ 6F, +66,

 

 

 

 

 

 

 

 

 

 

 

at 6x 6y 6x 6y (6°13)
rp \ VP“ 7 rpv 7
2 uv
Q=<pu>,Fi=<pu +p 7,Gi=(p 2 1
pv puv pv +p
.13 . .u(E+p)+q.. (V(E+p)+qy,
r0 1 i0 ‘
In I”
Fv=< >, GV=< 1 (6.1b)
tx)’ tYY
Lu'cxx +vrm (“Txy +vrm

138

 

where,

2 4 at
r — — th-u2——
3u( ) Fax

xx-

au 6v
_+_)

tw=ugw 0x

1")0

=—3u(V-V)+2u-al
3 6y

p 1 2 2
E=—+— u +v
7—1 29( )

q.=-k—-

9 3

qy=_k%§ (ans

 

and p is density; )1 is the dynamic viscosity; u and v are the x- and y-components of the

velocity, respectively; B is mechanical and thermal energy; P is pressure; T is

temperature; and 7 is the ratio of specific heats.

139

6.2 Euler DETE for Error Estimation of Navier-Stokes Solutions

By comparing the 2-D Navier-Stokes equation (Eq. (6.1)) with the 2-D Euler equation
(Eq. (5.1)), it is clear that the only difference between them lies in the viscous diffusion
terms. If error is mainly transported by the convection terms, not by diffusion terms, then
the Euler DETE developed in Chapter 5 may be adequate for the estimation of grid-
induced errors in the solutions obtained by Navier-Stokes solvers. This idea is examined

in this chapter.

By substituting the grid-independent Navier-Stokes solution into the linearized Euler

FV operator (developed in Chapter 5) on the coarse mesh, a residual will be generated:

Li. (U2?) = RSS (6.2)

Substituting the coarse-grid Navier-Stokes solution into the same Euler FV operator

on the coarse mesh, a residual will also be generated:

L'i.(U:‘S)=R:“S (6.3)
Subtracting these two will get:
Li (6”) = (R155 - R2“) (6.4)
and
em = Ug‘f - ufis (6.5)

140

 

 

where,

LE is the linearized Euler FV operator on the coarse mesh
NS . . . , .

Us.i rs the grid-independent Navrer-Stokes solution

U35 is the coarse-mesh Navier-Stokes solution

eNS is the grid-induced error in the Navier—Stokes solutions.

 

141

6.3 Test Problems

With the idea and procedures presented in Section 6.2, this section examines the
usefulness of Euler DETE for estimating the grid-induced errors in the solutions obtained
by Navier-Stokes solvers. Two test problems are selected for this purpose: one is the

clean airfoil case and the other one is the iced airfoil case.

6.3.1 Clean Airfoil Problem.

The first test problem chosen is the turbulent flow around a 944 clean airfoil. The
grid-independent and coarse-grid solutions are obtained using WIND software. The

freestream subsonic conditions are as follows: Mach Number M = 0.12, pressure P =

20.5 PSI, temperature T = 549.67 R, Reynolds number Re=3.5 x 106, and the angle of
attack = 4°. Spalan-Allrnaras was chosen as the turbulence model. (See Chi et al. (2004)

for more details.)

The grid scheme used is the multi-block structured mesh (Chi et al (2004)). The grid-
independent solution for this problem is obtained using 913 x 101 (inner block) + 125 x
21 (outer block) mesh. Because most of the physics happens in the inner block, only the
inner block is changed to get the coarse-grids. Two sets of coarse mesh (457 x 101 and
457 x 61 in the inner block) and the corresponding solutions are obtained for this test
problem. Figure 6.1 shows the mesh used to obtain the grid-independent solution. Since

we are only interested in the results of the inner block, the following figures only show

142

 

the results in the inner block. Figure 6.2 shows the grid-independent solution in x-
momentum. Figure 6.3 compares the local fine and coarse meshes around the leading
edge.

6.3.2 Iced Airfoil Problem.

The second test problem chosen is the turbulent flow around a 944 iced airfoil. The
grid-independent and coarse-grid solutions are obtained using WIND software. The

freestream subsonic conditions are as follows: Mach Number M = 0.12, pressure P =

20.5 PSI, temperature T = 549.67 R, Reynolds number Re=3.5 x 106, and the angle of
attack = 4°. Spalart-Allrnaras was chosen as the turbulence model; see Chi et al (2004)

for more details.

The grid scheme used is the multi-block structured mesh (Chi et al (2004)). The grid-
independent solution for this problem is obtained using 941 x 101 (inner block) + 125 x
21 (outer block) mesh. Because most of the physics happens in the inner block, so only
the inner block was changed to get the coarse-grids. Two sets of coarse mesh (471 x 101
and 471 x 61 in the inner block) and the corresponding solutions are obtained for this test
problem. Figure 6.4 shows the mesh used to obtain the grid-independent solution.
Figure 6.5 shows the grid-independent solution in x-momentum., and Figure 6.6

compares the local fine and coarse meshes around the leading edge.

143

 

6.4 Results

6.4.1 Clean Airfoil Problem.

The grid independent solution for the clean airfoil problem was generated using the

913 x 101 (inner block) + 125 x 21 (outer block) mesh. The solution obtained on this

grid is used to compute the residuals ( RES) in the DETEs (Eq. (6.4) and to calculate the

actual absolute and relative errors, which can be used to assess how well Euler DETEs

predict grid-induced errors on coarse meshes.

Figure 6.7 to Figure 6.8 show the residual (RES) for the solutions generated on the

following two coarser meshes: 457 x 61 and 457 x 101. Only the residuals for the
continuity equation (density) are shown. These residuals were computed by substituting
the grid-independent solution into Eq. (6.2). Thus, the residuals shown are exact. As
expected, these figures show the residual is high near the leading edge and the trailing

edge, and it increases as the number of grid cells decreases.

Figure 6.9 to Figure 6.12 show the absolute errors in the solutions generated on the
following two coarser meshes: 457 x 61, and 457 x 101. Only the absolute errors in the
computed solution for x-momentum and y-momentum are shown. The actual errors were
computed by using Eq. (6.5); that is by subtracting the coarse mesh solution from the
grid-independent solution at the center of each cell in the coarse mesh (e. g.,

e =(pu)8 —(pu)c ). Bilinear interpolation is used to transfer the grid-independent

144

 

solution to the center of each cell in the coarse mesh. The errors predicted by Euler
DETE were obtained by Eq. (6.4) with the exact residual (R:s — RES) used to compute
the grid-induced error. As expected, the error is high near the leading edge and the
trailing edge, and it increases as the number of grid cells decreases. Results obtained

show that the Euler DETEs were able to predict very well the errors in the Navier-Stokes

solutions, if the actual residual is used. .

6.4.2 Iced Airfoil Problem.

The grid independent solution for the iced airfoil problem was generated using the

941 x 101 (inner block) + 125 x 21 (outer block) mesh. The solution obtained on this

grid is used to compute the residuals (R313) in the DETEs (Eq. (6.4)) and to calculate the

actual absolute and relative errors, which can be used to assess how well Euler DETEs

predict grid-induced errors on coarse meshes.

Figure 6.13 to Figure 6.14 show the residual (R :8) for the solutions generated on the

following two coarser meshes: 471 x 61, and 471 x 101. Only the residuals for the
continuity equation (density) are shown. These residuals were computed by substituting
the grid-independent solution into Eq. (6.2). Thus, the residuals shown are exact. As
expected, these figures show the residual is high near the leading edge and the trailing

edge, and it increases as the number of grid cells decreases.

Figure 6.15 to Figure 6.18 show the absolute errors in the solutions generated on the

145

 

following two coarser meshes: 471 x 61, and 471 x 101. Only the absolute errors in the
computed solution for x-momentum and y-momentum are shown. The actual errors were
computed by using Eq. (6.5); that is by subtracting the coarse mesh solution from the

grid-independent solution at the center of each cell in the coarse mesh (e.g.,
e = (pu)g —(pu)c ). Bilinear interpolation is used to transfer the grid-independent
solution to the center of each cell in the coarse mesh. The errors predicted by Euler

DETE were obtained by Eq. (6.4) with the exact residual (RI:S — R :13) used to compute

the grid-induced error. Similar to the clean airfoil problem, the error is high near the
leading edge and the trailing edge, and it increases as the number of grid cells decreases.
Results obtained also show that the Euler DETEs are able to predict very well the errors

in the Navier-Stokes solutions, if the actual residual is used. .

146

 

6.5 Summary of Chapter 6

The usefulness of Euler DETE in predicting grid-induced errors in the solutions
obtained by Navier-Stokes solvers has been examined in this chapter using two test
problems: the clean airfoil and the iced airfoil problems. Results obtained show that error
predicted by Euler DETE matches very well with the actual error for the high-Reynolds-

number Navier-Stokes solutions.

147

 

-4OO

 

-500 0 500

 

(a) multi-block computational mesh

     

 

I 1 . I :1"
10 20 30 40
X

-10 O

(b) computational mesh around the airfoil

Figure 6.1 Grid-independent computational mesh (in Inches) for the clean airfoil

problem

148

 

1 i -11.oo

-0.5 0.0 0.5
X

Figure 6.2 Grid-independent solution in x-momentum for the clean airfoil problem

 

 

 

 

 

(a) grid-independent mesh (913 x 101) (b) coarse mesh (457 x 61)

Figure 6.3 Computational mesh around the leading edge for the clean airfoil problem

149

 

. . . _
\\\\\ \\
\\“:\\§A\\\\\\\\\\\‘\m
“q.\\\\\\\\\:\~:\\u::n ...,. ...,
mum—mm ..
“:5“ wl‘h‘m“?
‘WWV' m

..., ,~ ,
7,4." 3 .13, ...,...
I ‘~I; 'um: um-
94:35:555” _m
11’?" ”as”? ' ""“""'"""""
I,/II,;0 ”Aliyah”... mm
3;,“ IazzzaIlllnnmuu
‘11:; mm

15mm:
:7, 5: ‘.I ”ml!
I

I 4I
r: 4f
%’¢I ;III

III?”

 

 

(b) computational mesh around the airfoil

Figure 6.4 Grid-independent computational mesh (in Inches) for the iced airfoil problem

150

 

 

-015 ' 0 ' 0'5
x

(a) around the airfoil (inner block)

0.08
>- 0.00

-0.08

 

-0.08 0.00 0.08

X
(b) around the leading edge

Figure 6.5 Grid-independent solution in x-momentum for the iced airfoil problem

151

 

 

        

       

    

                                
   

  
   

I
III II I; II I
IIII IIII, IrIIII’lov

      

4' u
s. as

          

             

 

     
   

       

          

     
   

            

     

  

s. a I IIIIIIIIII/ I

I a, a, 42/5,

2. tea/424444343

:5 4,3222
.2........... .. gig. ...,...MHaMMaWMIMMaMaMaWa. . 4

2.. .2... ...,... 2,5225%: .

..‘3ssine“.fi.““..“....m.s............2._ Egg. 3’35, 3 . O
as....finena.“3%,”...yaswa._ g5... gag augegag .. E .

 
      

(a) grid-independent mesh (941 x 101)

1.2+

.Mu

111 #1

M11

.
t1)
1 1‘.

+1?
+1 .
» J1t1
1.11er

1

.r

I . 1.11111?

1.. 1 1 i
)1 1.11

1.111111:

Lr.
2 ..

all. 3.1....

. L .1)...
x .

r .1 11.17.14 .

 

 

 

0.04

0.00

.04

-0

(b) coarse mesh (471 x 61)

Figure 6.6 Computational mesh around the leading edge for the iced

oil problem

152

 

-0.2 0.0 0.2
X

 

(a) leading edge

0.05

  
 

>.

0.00 N

0.90 0.95

X
(b) trailing edge

Figure 6.7 Exact Residual R25 in the continuity equation for the clean airfoil problem:
457 x 61 case

153

>. 0.00

-0.20

 

-0.20 0.00 0.20
X

(a) leading edge

 

0.90 0.95
X

(b) trailing edge

Figure 6.8 Exact Residual RES in the continuity equation for the clean airfoil problem:
457 x 101 case

154

0.02

0.00

 

-0.02
-0.02 0.00 0.02 0.04

   

(a) actual error

0.02

0.00

-0.02
-0.02 0.00 0.02 0.04

   

(b) error predicted using Euler DETE

Figure 6.9 Absolute error in x-momentum for the clean airfoil problem: 457 x 61 case.

155

0.05 '-

0.00 L

 

-0.05 - 1 .23

 

 

 

-0.05 I 0.00 I 0.05
x

(a) actual error

0.05 '-

0.00 '-

-0.05

 

'1 r

-0.05 0.00 I 0.05
x

A

 

(b) error predicted using Euler DETE

Figure 6.10 Absolute error in y-momentum for the clean airfoil problem: 457 x 61 case.

156

0.02

0.00

-0.02

 

-0.02 0.00 0.02 0.04

(a) actual error

0.02

“I'EEQ‘”

>. 0.00

-0.02

   

-0.02 0.00 0.02 0.04

(b) error predicted using Euler DETE

Figure 6.11 Absolute error in x-momentum for the clean airfoil problem: 457 x 101 case.

r11
.1 11.1
an
9'?
U0

0.05 -

>. 0.00-— 1C

 

 

-0.05 i 0.00 ' 0.05
x

(a) actual error

0.05 -

>- 0.00- K

-0.05 -

 

 

-0.05 ' 0.00 i 0.05
x

(b) error predicted using Euler DETE

Figure 6.12 Absolute error in y-momentum for the clean airfoil problem: 457 x 101 case.

158

 

 

0.20

>. 0.00

-0.20

 

-0.20 0.00 0.20
X

(a) leading edge

 

0.10
0.05

0.00 >-

-0.05

 

0.90 1.00
X

(b) trailing edge

Figure 6.13 Exact Residual RES in the continuity equation for the iced airfoil problem:
471 x 61 case

159

 

0.20

>. 0.00

-O.20

 

-O.2O 0.00 0.20

X
(a) leading edge

 

0.10

 

0.05
>.
0.00
-0.05
0.90 1 .00
X
(b) trailing edge

Figure 6.14 Exact Residual R;5 in the continuity equation for the iced airfoil problem:
471 x 101 case

160

0.02

>. 0.00

-0.02

 

-0.04
-0.02 0.00 0.02 0.04

X

(a) actual error

0.02

>. 0.00

-0.02

 

-0.04

-0.02 0.00 0.02 0.04
X

(b) error predicted using Euler DETE

Figure 6.15 Absolute error in x-momentum for the iced airfoil problem: 471 x 61 case.

 

 

 

 

 

 

0.10 -
0.05 " _ ,
> .. "1. \‘7 ,/ r
f L, '
0.00 - \
(T:’<:‘A‘\’ \ \6‘
-0.05 — ,> ‘ 7

 

-0.05' 0.00 ' 0.05 ' 0.10 i
x

(a) actual error

 

x 1

 

ii 20.0
' 16.
0.10 - 1:1 12.:
13 8.9
1 5.3
; 1.6
,1, E 1 '
I — i . I i -
o 05 g! i _.
’~ - ~ 117: .
>- if '52;- ’8' :
0.00 - g.
1" . I
‘14.- x “7‘4 ~- I
‘ 15
-0.05 - I

 

 

-0.05‘ 0.00 ' 0.05 ' 0.10 '
x

(b) error predicted using Euler DETE

Figure 6.16 Absolute error in y-momentum for the iced airfoil problem: 471 x 61 case.

162

0.04

>' 0.00

 

-0.04

 

-0.04 0.00 0.04 0.08
X

(a) actual error

1

1‘!’

0.04 a

I

I
>- 0.00
-0.04

   

-0.04 0.00 0.04 0.08
X

(b) error predicted using Euler DETE

Figure 6.17 Absolute error in x-momentum for the iced airfoil problem: 471 x 101 case.

163

0.04

0.00

 

-0.04

 

-0.04 0.00 0.04 0.08
X

(a) actual error

0.04

0.00

-0.04

   

-0.04 0.00 0.04 0.08
X

(b) error predicted using Euler DETE

Figure 6.18 Absolute error in y-momentum for the iced airfoil problem: 471 x 101 case.

CHAPTER 7
CONCLUSIONS

With computational fluid dynamics (CFD) becoming more accepted and more
widely used in industry for design and analysis, there is increasing demand for not just
more accurate solutions, but also error bounds on the solutions. One major source of
error is from the grid or mesh. A number of methods have been developed and evaluated
to quantify errors in solutions of PDEs that arise fi'om poor-quality or insufficiently fine
grids/meshes. For PDEs of interest to CFD, it has been shown that the error at one
location could be generated elsewhere and then transported there, and thus is not a
function of the local mesh quality and the local solution. So, a transport equation for

error is needed to understand the generation and evolution of errors.

This research developed and evaluated discrete-error-transport equations (DETEs)
for estimating grid-induced errors in CFD. The DETE is derived fi'om the discrete finite-
difference/finite-volume fi'amework instead of the continuous finite element framework.
In Chapter 2, the concept of error transport is illustrated by a 2D laminar flow problem.
The concept and derivation of the DETE is demonstrated by a simple model equation: the

one-dimensional advection-diffusion equation.
The discrete-error-transport equation was first developed for a few 1-D test problems

in Chapter 3: the advection-diffusion equation (linear and steady), the wave equation

(linear and unsteady), the inviscid Burgers equation (nonlinear and unsteady) and the

165

 

 

steady Burgers equation (nonlinear and steady). The usefulness of the discrete-error-
transport equations in estimating grid-induced errors in finite-difference and finite-
volume methods was evaluated on these test problems. The residual in the discrete-error-
transport equation was modeled by the leading term of the remainder in the modified
equation. For each test problem, solutions were generated for a number of meshes,
including one that is grid-independent. The grid-independent solution for each test
problem is used to assess the accuracy with which the residuals are modeled by the
modified equation and how well the discrete-error—transport equation predicts grid-

induced errors.

Results obtained show that the modified equation is useful in modeling the residual if
the mesh is sufficiently fine. Results obtained also show the discrete-error-transport
equation to be able to predict the error perfectly (up to the last significant digit) if the
actual residual is used. This is true even for the unsteady, nonlinear test problem with a
moving weak discontinuity in the flow domain. The linearization procedure used to
derive discrete-error-transport equations for nonlinear FD/FV equations — whether
conservative or non-conservative — was found to be unimportant because the error-
propagation speed is computed by the FD/FV equations, and not by the discrete-error

transport equation.

As noted, the major challenge in implementing both the continuous and the discrete

error-transport equations is in modeling the residual. One focus of this study is on the

modeling of the residual in the discrete error—transport equation (DETE) for FD and FV

166

 

'11

methods. Since the modified equation and Taylor-series expansions are not useful in
modeling the residual when the grid is coarse (the typical situation), an alternative
approach for residual modeling based on data mining is explored in Chapter 4. The
approach involves two steps. The first is to study the behavior of the residual by
evaluating the “actua ” residual created by a variety of poor quality meshes in a
systematic way. The second step is to model the residual based on the understanding
gained on the behavior of the residual. The residuals generated by the numerical
experiments were fitted by using least-square minimization. Preliminary results show the
power-law to produce the best fit. However, a more extensive database is needed before

this approach can be expected to yield a more generally applicable model of the residual.

So far, the discrete-error transport equation (DETE) has only been developed and
demonstrated for simple model equations such as the advection diffusion equations and
the inviscid Burger equation. Also, it has only been applied on structured grids, where
finite-difference and finite-volume methods can lead to identical algebraic analogs of
PDES. Thus, a DETE for a more complicated set of equations in the framework of
arbitrary elements is developed in Chapter 5. The set of DETE for the system of Euler
equations has been developed and evaluated on a few test problems governed by the
Euler equations: the oblique shock problem and the inviscid flow past an airfoil problem.
Results obtained show that DETE concept to be applicable to coupled, quasi-linear PDEs

with weak solutions.

The usefulness of Euler DETE in predicting grid-induced errors in the Navier-Stokes

167

 

 

solutions has been examined in Chapter 6. Results obtained show that error predicted by
Euler DETE matches very well with the actual error for the high-Reynolds-number

Navier-Stokes solutions.

It has been shown that the DETE developed in this research can predict the grid-
induced errors perfectly if the actual residual is used in the DETE. Since the actual
residual is generally unknown, the main challenge in using DETE concept is to develop
models for the residual. This research proposed a residual modeling method based on
data mining. Though the results seem promising for the specific test problems studied,

more studies are needed before a general model for the residual can be constructed.

168

 

 

 

 

 

 

APPENDIX ‘

169

APPENDIX A

Publications Associated with this Work '

Qin, Y. , Keller, P.S., Sun, R.L., Hernandez, E.C., Perng, C.-Y., Trigui, N., Han, Z.,
Shen, F.Z., Shieh, T., Shih, T.I-P., “Estimating Grid-Induced Errors in CF D by
Discrete-Error-Transport Equations”, AIAA Paper 2004 -0656, January 2004

Qin, Y. and Shih, T.I.P., “Estimating Grid-Induced Errors by a Discrete-Error-
Transport Equation: Modeling of the Residual”, AIAA Paper 2003-3850, June 2003

Qin, Y. and Shih, T.I.P., “A Method for Estimating Grid-induced Errors in Finite-
Difference and Finite-Volume Methods”, AIAA Paper 2003-0845, January 2003

Qin, Y. and Shih, T.I.P., “A Discrete Transport Equation for Error Estimation in
CFD”, AIAA Paper 2002-0906, January 2002.

170

 

 

 

BIBILOGRAPHY

 

171

BIBILOGRAPHY

Adjerid, S., F laherty, J .E. and Babuska, I., “A Posteriori Error Estimation for the Finite
Element Method-of-Lines Solution of Parabolic Problems,” Mathematical Models and
Methods in Applied Sciences, Vol. 9, No. 2, 1999, pp. 261-286

Adjerid, S., Belguendouz, B., and Flaherty, J .E., “A Posteriori Finite Element Error
Estimation for Diffusion Problems,” SIAM Journal on Scientific Computing, Vol. 21,
1999, pp. 728-746

Ainsworth, M. and Oden, J .T., “ A Procedure for a-Posteriori Error Estimation for h-p
Finite Element Methods,” Computer Methods in Applied Mechanics and Engineering,
Vol. 101, 1992, pp. 73-96

Babuska, I. and Aziz, A.K., “On the Angle Condition in the Finite Element Method,”
SIAM Journal on Numerical Analysis, Vol. 13, No. 2, April 1976, pp. 215-226

Babuska, I. and Rheinboldt, W., “Error Estimates for Adaptive Finite Element
Computations,” SIAM Journal on Numerical Analysis, Vol. 15, No. 4, 1978, pp. 736-
754.

Babuska, I., Strouboulis, T., and Upadhyay, OS, “A Model Study of the Quality of a
Posteriori Error Estimator for Linear Elliptic Problems: Error Estimation in the Interior
of Patchwise Uniform Grid of Triangles,” Computer Methods in Applied Mechanics and
Engineering, Vol. 114, 1994, pp. 307-378.

Babuska, I., Strouboulis, T., Upadhyay, CS, and Gangaraj, S.K., “A Posterior
Estimation and Adaptive Control of the Pollution Error in the h-Version of the Finite
Element Method,” Intl. Journal for Numerical Methods in Engineering, Vol. 38, 1995,
pp. 4207-4235.

Babuska, 1., Strouboulis, T., Gangaraj, SK, and Upadhyay, C.S., “Pollution Error in the
h-Version of the Finite Element Method and Local Quality of the Recovered
Derivatives,” Computer Methods in Applied Mechanics and Engineering, Vol. 140, 1997,
pp. 1-37.

Babuska, I., Ilrlenburg, F., Strouboulis, T., and Gangaraj, S.K., “A Posterior Error
Estimation for Finite Element Solutions of Helmholtz Equation—Part II: Estimation of
the Pollution Error,” Intl. Journal for Numerical Methods in Engineering, Vol. 40, 1997,
pp. 3883-3900.

Bank, RE. and Weiser, A., “Some a-Posteriori Error Estimators for Elliptic Partial
Differential Equations,” Mathematics of Computation, Vol. 44, 1985, pp. 283-301

Bank, RE. and Welfert, D.B., “Some A Posteriori Error Estimates for the Stokes

172

 

 

.
ii

Equations: A Comparison,” Computer Methods in Applied Mechanics and Engineering,
Vol. 87, 1990, pp. 323-340

Barth, T.J. and Larson, M.G., “A Posteriori Error Estimation for Discontinuous Galerkin
Approximations of Hyperbolic Systems,” 1St International Conference on DG Methods,
May 24-26, 1999, NASA Report NAS-99-010

Berzins, M., “A Solution-Based Triangular and Tetrahedral Mesh Quality Indicator,”
SIAM Journal of Scientific Computing, Vol. 19, No. 6, 1998, pp. 2051-2060

Berzins, M., “Mesh Quality — Geometry, Error Estimates or Both?” Engineering and
Computers, Vol. 15, 1999, pp. 236-247

Berzins, M., “Solution-Based Mesh Quality Indicators for Triangular and Tetrahedral
Meshes,” International Journal of Computational Geometry & Applications, Vol. 10(3),
2000, pp. 333-346

Brackbill, J .U., “An Adaptive Grid with Directional Control,” Journal of Computational
Physics, Vol. 108, 1993, pp. 38

Carey, G.F., Computational Grids: Generation, Adaptation, and Solution Strategies,
Taylor & Francis, Washington, DC, 1997.

Celik, 1., Chen, C.J., Roache, P.J., and Scheuer, G., Editors, Symposium on
Quantification of Uncertainty in Computational Fluid Dynamics, FED-Vol. 158, ASME
Fluids Engineering Division, Summer Meeting, Washington DC, 1993.

Celik, I. and Zhang, W.-M., “Calculation of Numerical Uncertainty Using Richardson
Extrapolation: Application to Some Simple Turbulent Flow Calculations,” ASME
Journal of Fluids Engineering, Vol. 117, September 1995, pp. 439-445.

Celik, I. And Karatekin, 0., “Numerical Experiments on Application of Richardson
Extrapolation with Nonuniforrn Grids,” ASME Journal of Fluids Engineering, Vol. 119,
1997, pp. 584-590

Celik, I., Hu, G., and Badeau, A., “Further Refinement and Bench Marking of a Single-
Grid Error Estimation Technique,” AIAA Paper 2003-0628, Jan. 2003.

Chi, X., Zhu, B., Shih, T.I-P., Addy, HE, and Choo, Y.K., “CFD Analysis of the
Aerodynamics of a Business-Jet Airfoil with Leading-Edge Ice Accretion,” AIAA Paper
2004-0560, Aerospace Sciences Meeting, Reno, January 2004.

Coleman, H.W., Stern, E, Di Mascio, A., and Campana, E., “The Problem with
Oscillatory Behavior in Grid Convergence Studies,” ASME Journal of Fluids
Engineering, Vol. 123, 2001, pp. 438-439.

173

 

 

 

Cosner, R.R., Oberkampf, W.L., Rahaim, CR, and Shih, T.I-P., “AIAA Committee on
Standards for Computational Fluid Dynamics — Status and Plans,” AIAA Paper 2004-
0654, Aerospace Sciences Meeting, Reno, January 2004.

Davis, SF. and Flaherty, J .E., “An Adaptive Finite Element Method for Initial-Boundary
Value Problems for Partial Differential Equations,” SIAM Journal on Scientific and
Statistical Computing, Vol. 3, 1982, pp. 6-27

Demuren, A.O., and Wilson, R.V., “Estimating Uncertainty in Computations of Two-
Dimensional Separated Flows,” ASME Journal of Fluids Engineering, Vol. 116, June,
1994, pp. 216-220

Eca, L. and Hoekstra, M., “An Evaluation of Verification Procedures for CFD
Applications,” 24th Symposium on Naval Hydrodynamics, Fukuoka, Japan, July 8-13,
2002

Ekeland, I. And Teman, R., “Convex Analysis and Variational Problems,” Amsterdam:
North-Holland Pub. Co., 1976

Eriksson, K., “Improved Accuracy by Adapted Mesh-Refinements in the Finite Element
Method,” Mathematics of Computation, Vol. 44, 1985, pp. 321-343

Eriksson, K. and Johnson, C., “An Adaptive Finite Element Method for Linear Elliptic
Problems,” Mathematics of Computation, Vol. 50, 1988, pp. 361-383

Eriksson, K. and Johnson, C., “Adaptive Finite Element Methods for Parabolic Problems
I: A Linear Model Problem,” SIAM Journal on Numerical Analysis, Vol. 28, 1991, pp.
43-47

Ferziger, J .H., “A note on numerical accuracy,” Special Note, International Journal for
Numerical Methods in Fluids, Vol. 8, pp. 995-996, 1988

Ferziger, J .H., “Estimation and Reduction of Numerical Error”, Paper presented at the
forum on “Methods for Estimating Uncertainty Limits in Fluid Flow Calculations,”
ASME Winter Annual Meetings, San Francisco, CA, December 10-15, 1989

Ferziger, J.H., “Estimation and Reduction of Numerical Error,” in Symposium on
Quantification of Uncertainty in Computational Fluid Dynamics, F ED-Vol. 158, ASME
Fluids Engineering Division, Summer Meeting, Washington DC, 1993, pp. 1-8.

Field, D.A., “Qualitative Measures for Initial Meshes,” International Journal for
Numerical Methods in Engineering, Vol. 47, 2000, pp. 887-906

Freitas, C.J., Ghia, U., Celik, 1., Roache, P., and Raad P., “ASME’s Quest to Quantify
Numerical Uncertainty,” AIAA Paper 2003-0627, J an, 2003

174

 

 

fry-"1.";- =
i: .

 

Gu, X., Schock, H.J., Shih, T.I-P., Hernandez, E.C., Chu, D., Keller, PS, and Sun, R.L.,
“Grid-Quality Measures for Structured and Unstructured Meshes,” AIAA Paper 2001-
0652, Jan. 2001.

Gu, X. and Shih, T.I-P., “Differentiating between Error Source and Error Location in
Solution-Adaptive Mesh Refinement,” AIAA Paper 2001-2660, CFD Conference, June
2001.

Gu, X., “Grid-Quality Measures for Error Estimation and Solution-Adaptive Mesh
Refinement in CFD”, Ph.D Dissertation, 2002, Michigan State University

Haworth, D.C., El Tahry, SH. and Huebler, M.S., “A Global Approach to Error
Estimation and Physical Diagnosis in Multidimensional Fluid Dynamics,” International
Journal of Numerical Methods in Fluids, Vol. 17 , No. 1, 1993, pp. 75-97.

Ilinca, C., Zhang, X.D., Trepanier, J.-Y., and Carnarero, R., “A Comparison of Three
Error Estimation Techniques for Finite-Volume Solutions of Compressible Flows,”
Computer Methods in Applied Mechanics and Engineering, Vol. 189, 2000, pp. 1277-
1294

Jasak, H. and Gosman, A.D., “Local Problem Error Estimate in Finite Volume
Discretization,” Submitted to Computer Methods in Applied Mechanics and Engineering,
1998

Jasak, H. and Gosman, A.D., “Residual Error Estimate for the Finite Volume Method,”
Submitted to International Journal for Numerical Methods in Fluids, 1999

Johnson, C., Rannacher, R. and Boman, M., “Numerics and Hydrodynamics Stability
Theory: Towards Error Control in CFD,” SIAM Journal on Numerical Analysis, Vol. 32,
1995, pp. 1058-1079

Johnson, C. and Szepessy, A., “Adaptive Finite Element Methods for Conservation Laws
Based on a Posterior Error Estimates,” Communications on Pure and Applied
Mathematics, Vol. 48, No. 3, 1995, pp. 199-234

Knupp, P., “Mesh Generation using Vector-Fields,” Journal of Computational Physics,
Vol. 119, 1995, pp. 142

Krizek, M., “On the Maximum Angle Condition for Linear Tetrahedral Elements,” SIAM
Journal on Numerical Analysis, Vol. 29, 1992, pp. 513-520

Lepape, M.C. and Jacquotte, O.-P., “Acceleration Techniques for a 3D Structured
Multiblock Mesh Optimization Method,” Numerical Grid Generation in Computational
Field Simulation and Related Fields, N.P. Weatherill, P.R. Eiseman, J. Hauser and J .F.
Thompson (Eds), p. 95, Proceedings of the 4th International Grid Conference, Pineridge

175

 

 

Mm. u
l

‘I

‘4

V

Press Limited: Swansea Wales (UK), 1994

Liu, A. and Joe, B., “Relationship between Tetrahedron Shape Measures,” BIT, Vol. 34,
1994, pp. 268-287

L0, S.H., “Optimization of Tetrahedral Meshes based on Element Shape Measures,”
Computers & Structures, Vol. 63, No. 5, 1997, pp. 951-961

Mackenzie, J ., Sonar, T., and Suli, E., “Adaptive Finite Volume Methods for Hyperbolic
Problems,” Mathematics of Finite Elements and Applications, edited by J .R. Whiteman,
John Wiley and Sons, New York, 1994, Chapter 19.

Mitchell, A.R., Phillips, G. and Wachspress, E., “Forbidden Shapes in the Finite Element
Method,” Journal of the Institute for Mathematics and Applications, Vol. 8, 1971, pp.
260-269

Oberkampf, W.L., Blottner, F.G., and Aeshlirnan, D., “Methodology for Computational
Fluid Dynamics Code Verification and Validation,” AIAA Paper 95-2226, June 1995.

Oden, J .T., Demkowicz, L., Rachowicz, L., and Hardy, O., “Toward a Universal h-p
Adaptive Finite Element Strategy. Part I: Constrained Approximation and Data
Structure,” Computer Methods in Applied Mechanics and Engineering, Vol. 77, 1989,
pp. 79-112

Oden, J .T., Demkowicz, L., Rachowicz, L., and Westennann, I.., “Toward a Universal h-
p Adaptive Finite Element Strategy. Part II: A Posteriori Error Estimation,” Computer
Methods in Applied Mechanics and Engineering, Vol. 77, 1989, pp. 113-180

Oden, J.T., “Error Estimation and Control in Computational Fluid Dynamics,”
Mathematics of Finite Elements and Applications, edited by J .R. Whiteman, John Wiley
and Sons, New York, 1994, Chapter 1.

Owen, S.J., “Constrained Triangulation: Application to Hex-Dominant Mesh
Generation,” Proceedings of 8th International Meshing Roundtable, South Lake Tahoe,
CA, U.S.A., pp. 31-41, October 1999.

Peraire, J., Vahdati, M., Morgan, K., and Zienkiewicz, O.C., “Adaptive Remeshing for
Compressible Flow Calculations,” J. of Computational Physics, Vol. 22, 1976, pp.131-
149.

Richardson, L.F., “The Approximate Arithmetical Solution by Finite Differences of
Physical Problems Involving Differential Equations, with an Application to the Stresses
In a Masonary Dam,” Transactions of the Royal Society of London, Ser.A, Vol. 210,
1910, pp.307-357

Richardson, LR, and Gaunt, J .A., “The Deferred Approach to the Limit,” Philos. Trans.

176

 

 

 

R. Soc. London Ser. A., Vol. 226, 1927, pp. 299-361

Rippa, 8., “Long and Thin Triangles can be Good for Linear Interpolation,” SIAM
Journal of Numerical Analysis, Vol. 29, 1992, pp. 257-270

Roache, P.J., Computational Fluid Dynamics, Hermosa Publishers, PO. Box 8172,
Albuquerque, N.M., 1972

Roache, P.J., “A Method for Uniform Reporting of Grid Refinement Studies,”
Proceedings of Quantification of Uncertainty in Computational Fluid Dynamics, ASME
Fluids Engineering Division Spring Meeting, ASME Pub]. No. FED-Vol. 158,
Washington DC, June, 1993

Roache, P., Verification and Validation in Computational Science and Engineering,
Hermosa Publishers, Albuquerque, New Mexico, 1998.

Robinson, J ., “Some New Distortion Measures for Quadrilaterals,” Finite Elements in
Analysis and Design, Vol. 3, 1987, pp. 183-197

Robinson, J ., “Distortion Measures for Quadrilaterals with Curved Boundaries,” Finite
Elements in Analysis and Design, Vol. 4, 1988, pp. 115-131

Robinson, J ., “Quadrilateral and Hexahedron Shape Parameters,” Finite Elements in
Analysis and Design, Vol. 16, 1994, pp. 43-52

Shih, T.I-P., Gu, X., and Chu, D., “Grid-Quality Measures and Error Estimates,” in
Numerical Grid Generation in Computational Field Simulations, B.K. Soni, J .F.
Thompson, J. Hauser, and RR. Eiseman, editor, Mississippi State University, ISSG,
2000, pp. 799-808.

Shyy, W., and Garbey, M., “A Least Square Extrapolation Method for Improving
Solution Accuracy of PDE Computations,” J. of Computational Physics, Vol. 186, 2003,
pp.l-23.

Sonar, T., “Strong and Weak Norm Refinement Indicators Based on the Finite Element
Residual for Compressible Flow Computations,” Impact of Computing in Science and
Engineering, Vol. 5, 1993, pp. 111-127.

Sorokin, A.M., “Structured and Unstructured Grid Generation: An Evolutionary
Approach,” Numerical Grid Generation in Computational Field Simulation and Related
Fields, N.P. Weatherill, P.R. Eiseman, J. Hauser and J.F. Thompson (Eds), p. 287,
Proceedings of the 4th International Grid Conference, Pineridge Press Limited: Swansea
Wales (UK), 1994

Stern, R, Wilson, R.V., Coleman, H.W., and Paterson, E.G., “Comprehensive Approach

177

 

 

to Verification and Validation of CFD Simulations --- Part 1: Methodology and
Procedures,” Journal of Fluids Engineering, Vol. 123, 2001, pp. 793-802

Strouboulis, T. and Hague, K.A., “Recent Experiences with Error Estimation and
Adaptivity 1: Review of Error Estimators for Scalar Elliptic Problems,” Computer
Methods in Applied Mechanics and Engineering, Vol. 97, 1992, pp. 399-436

Suli, E., “A Posteriori Error Analysis and Adaptivity for Finite Element Approximations
of Hyperbolic Problems,” In An Introduction to Recent Developments in Theory and
Numerics for Conservation Laws, Vol. 5 of LNCSE, pp. 112-194, Springer-Verlag,
Heidelberg, 1998

Synge, J .L., “The Hypercircle in Mathematical physics: a Method for the Approximate
Solution of Boundary Value Problems,” Cambridge [Eng] University Press, 1957

Van Straalen, B.P., Simpson, RB, and Stubley, G.D., “A Posteriori Error Estimation for
Finite-Volume Simulations of Fluid Flow Transport,” Proceedings of the 3rd Annual
Conference of the CFD Society of Canada, Vol. 1, Baniff, Alberta, June 1995.

Venditti, DA. and Darrnofal, D.L., “Adjoint Error Estimation and Grid Adaptation for
Functional Outputs: Application to Quasi-One-Dimensional Flow,” Journal of
Computational Physics, Vol. 164, 2000, pp. 204-227.

Verfurth, R., “A Posteriori Error Estimates for Nonlinear Problems. Finite Element
Discretization of Elliptic Equations,” Mathematics of Computation, Vol. 62, No. 206,
1994, pp. 445-475

Verfurth, R., “A Posteriori Error Estimation and Adaptive Mesh-Refinement
Techniques,” Journal of Computational and Applied Mathematics, Vol. 50, 1994, pp. 67-
83

Wa, K. and Chen, Z., “A Simple and Effective Mesh Quality Metric for Hexahedral and
Wedge Elements,” Proceedings of 9‘h International Meshing Roundtable, Sandia National
Laboratories, October 2000, pp. 325-333

Warming, R. F. and Hyett, B. J ., “The Modified Equation Approach to the Stability and
Accuracy Analysis of Finite-difference Methods,” J. Comp. Physics, Vol. 14, 1974, pp.
159-179.

Wilson, R.V., Stern, R, Coleman, H.W., and Paterson, E.G., “Comprehensive Approach
to Verification and Validation of CFD Simulations --- Part 2: Application for Rans
Simulation of a Cargo/Container Ship,” Journal of Fluids Engineering, Vol. 123, 2001,
pp. 803-810

Wilson, RV. and Stern, R, “Verification and Validation for RANS Simulation of a Naval
Surface Combatant,” AIAA Paper 2002-0904, Jan. 2002.

178

 

 

Zhang, X.D., Trepanier, J .-Y., and Camarero, R., “An a Posterior Error Estimation
Method Based on Error Equations,”, AIAA Paper 1997-1889., 1997, pp. 383-397

Zhang, X.D., Trépanier, J .-Y., and Camarero, R., “A Posteriori Error Estimation for
Finite-Volume Solutions of Hyperbolic Conservation Laws,” Computational Methods in
Applied Mechanics and Engineering, Vol. 185, 2000, pp. 1-19.

Zhang, X.D., Pelletier, D., T répanier, J .-Y., and Camarero, R., “Numerical Assessment of
Error Estimators for Euler Equations,” AIAA Journal, Vol. 39, No. 9, 2001, pp. 1706-
1715.

Zienkiewicz, O.C. and Zhu, J .Z., “A Simple Error Estimator and Adaptive Procedure for
Practical Engineering Analysis,” International J. for Numerical Methods in Engineering,
Vol. 24, 1987, pp. 337-357.

Zlamal, M., “On the Finite Element Method,” Numerical Math., Vol. 12, 1968, pp. 394-
409

179

 

‘

Q‘

.l‘