Download Numerical Methods - Computational Fluid Dynamics - Lecture Notes and more Study notes Fluid Dynamics in PDF only on Docsity!
Numerical Methods for Computational Fluid Dynamics
Introduction
Before considering the numerical analysis of the differential equations of fluid mechanics it is
necessary to examine in a general way the methods used for the numerical analysis of differential
equations. Any problem posed in terms of a differential equation must include a set of boundary
conditions.
The boundaries may be closed or open in a mathematical sense. If there is an open boundary, it
is not necessary to specify boundary conditions on that open boundary. In physical problems an
open boundary occurs most often when time is the variable and conditions need not be specified
at a future time.
In the general classification of second-order ordinary differential equations, the standard form of
the equations is given as follows.
2 2
2
2
Fu G x
u E t
u D x
u C x t
u B t
u A
In this equation, the coefficients A to G may be zero, constant, functions of x and y only, or
functions of x, y, and u. In addition, the coefficients A to E may be functions of ∂u/∂t or ∂u/∂x as
well. Here we have used x and t as independent variables, implying space and time. However,
any two physical or mathematical variables can be used. General second-order partial differential
equations are classified as parabolic, elliptic or hyperbolic depending on the following relationship
among the coefficients:
If B^2 – 4 A C > 0, the equation is hyperbolic.
If B^2 – 4 A C = 0, the equation is parabolic.
If B^2 – 4 A C < 0, the equation is elliptic.
These partial differential equations (PDE) have nothing to do with conic sections. It is just that
the conditions for the different behavior of the PDEs are the same as the equations that
determine the kind of curve that results from a plot of the equation At 2 + Bxt + Cx^2 + Dt + Ex + G
= 0.
You should be able to verify the classifications of common partial differential equations that are
given below.
The wave equation, (^2)
2 2 2
2
x
u c t
u
, is hyperbolic.
The conduction equation, (^2)
2
x
u
t
u
, is parabolic.
Laplace’s equation, 2 0
2
2
2
y
u
x
u , is elliptic.
Poission’s equation, 2 ( , , )
2
2
2 Gx y u y
u
x
u
, is elliptic.
The classification of these partial differential equations has implications for the kinds of boundary
conditions needed. Elliptic problems require a closed boundary with boundary conditions
specified at all points of the boundary. Hyperbolic and parabolic equations have an open
boundary. The coordinate for this boundary is time and the nature of the open boundary matches
our common sense that we cannot have future results affecting earlier ones, except in science
fiction. In addition, hyperbolic equations provide a set of space curves known as characteristics.
The solutions of hyperbolic equations propagate along these characteristic curves and this can
provide an important constraint on numerical solutions of these equations.
Problems with more than two independent variables have a behavior that is a mixture of the
classifications above. For example the conduction equation in two space dimensions,
2
2
2
2
y
u
x
u
t
u
has elliptic behavior in the x-y plane, but parabolic behavior in the x-t and
y-t planes. Fluid dynamics equations can have a variety of behaviors. The steady Navier-Stokes
equations are elliptic equations, requiring boundary conditions at all parts of the region defining
the flow. The transient Navier-Stokes equations behave in mixed ways. Subsonic flows have
elliptic behavior in the space dimensions, but parabolic behavior in the time direction. Supersonic
flows have a hyperbolic nature in the time domain, but elliptic behavior in the space domains only.
The differences in the behavior of the differential equations require differences in the numerical
approaches to solving different flow environments.
Of Various kinds of boundary conditions can arise in fluid dynamics and related problems. If the
dependent variable is specified at the boundary, the boundary condition is called a Dirichlet
boundary condition or a boundary condition of the first kind. If the gradient of the dependent
variable is specified, the boundary condition is called a Neumann boundary condition or a
boundary condition of the second kind. When the boundary conditions are specified as a
relationship between the dependent variable and its gradient, the boundary conditions are
classified as being of the third kind or mixed.
The three different types of boundary conditions can be illustrated in a heat transfer problem
where we are solving a differential equation for temperature. A specified boundary temperature
would be a Dirichlet boundary condition. A specified boundary heat flux, which is proportional to
the temperature gradient by Fourier’s Law, would be a Neumann or second-kind boundary
condition. Convection boundary conditions are mixed – or third-kind – boundary conditions. In
such a boundary condition, the external convective heat flux is set equal to the internal
conductive heat flux. This gives a mixed boundary condition such as h (^) conv(T∞ - T) = -k∂T/∂x at x =
All numerical methods convert the differential equation and boundary conditions into a set of
simultaneous linear algebraic equations. The differential equation provides an accurate
description of the dependent variable at any points in the region where the equation applies. The
numerical approach provides approximate numerical values of the dependent variable at a set of
discrete points in the region. The values of the dependent variable at these points are found by
solving the simultaneous algebraic equations.
Two fundamentally different approaches are used to derive the algebraic equations from the
differential equations: finite differences and finite elements. In the finite-difference approach,
numerical approximations to the derivatives occurring in the differential equations are used to
replace the derivatives at a set of grid nodes. In the finite-element approach, the region is divided
into elements and an approximate behavior for the dependent variable over the small element is
For steady-state problems, the n superscript is omitted. For problems with only one or two space
dimensions two or one of the directional subscripts may be omitted. The general use of the
notation remains. The subscripts (and superscript) on the dependent variable represent a
particular point in the region, (x i , yj , z k, t n ) where the variable is defined.
Finite-difference Expressions Derived from Taylor Series
The Taylor series provides a simple tool for deriving finite-difference approximations. It also gives
an indication of the error caused by the finite difference expression. Recall that the Taylor series
for a function of one variable, f(x), expanded about some point x = a, is given by the infinite
series,
3 3
3 2 2
2
x a dx
d f x a dx
d f x a dx
df f x f a x a x a x a
[2-5]
The “x = a” subscript on the derivatives reinforces the fact that these derivatives are evaluated at
the expansion point, x = a. We can write the infinite series using a summation notation as
follows:
(^)
0
n
n
x a
n
n x a dx
d f
n
f x [2-6]
In the equation above, we use the definitions of 0! = 1! = 1 and the definition of the zeroth
derivative as the function itself. I.e., d 0 f/dx^0 | (^) x=a = f(a).
If the series is truncated after some finite number of terms, say m terms, the omitted terms are
called the truncation error. These omitted terms are also an infinite series. This is illustrated
below.
Termsused Truncation error
x a dx
d f
n
x a dx
d f
n
f x nm
n
xa
n
m n
n
n
xa
n
n
(^)
(^0) (^1)
[2-7]
In this equation the second sum represents the truncation error, εm , from truncating the series
after m terms.
1
n m
n
xa
n
n
m x a dx
d f
n
[2-8]
The theorem of the mean can be used to show that the infinite-series truncation error can be
expressed in terms of the first term in the truncation error, that is
1 1
1 ( - ) ( 1 )!
m
x
m
m
m x a dx
d f
m
[2-9]
Here the subscript, “x = ξ”, on the derivative indicates that this derivative is no longer evaluated at
the known point x = a, but is to be evaluated at x = ξ, an unknown point between x and a. Thus,
the price we pay for reducing the infinite series for the truncation error to a single term is that we
lose the certainty about the point where the derivative is evaluated. In principle, this would allow
us to compute a bound on the error by finding the value of ξ, between x and a, that made the
error computed by equation [2-9] a maximum. In practice, we do not usually know the exact
functional form, f(x), let alone its (m+1) th^ derivative.
In using Taylor series to derive the basic finite-difference expressions, we start with uniform one-
dimensional grid spacing. The difference, ∆xi , between any two grid points is the same and is
given the symbol, h. This uniform grid can be expressed as follows.
∆xi = xi – xi-1 = h or x (^) i = x 0 + ih for all i = 0,…,N [2-10]
Various increments in x at any point along the grid can be written as follows:
xi+1 – xi-1 = xi+2 – xi = 2h xi-1 – xi = xi – xi+1 = –h x (^) i-1 – xi+1 = xi – xi+2 = –2h [2-11]
Using the increments in x defined above and the notation fi = f(xi ) the following Taylor series can
be written using expansion about the point x = x (^) i to express the values of f at some specific grid
points, xi+1 , x i-1 , xi+2 and xi-2. The conventional Taylor series expression for f(x) in equation [2-5]
can be adapted for use in finite differences by writing an expansion equation about a particular
grid point, x = xi , to determine the value of f(x) at another grid point, x i+k. From equation [2-10],
we see that xi+k = xi + kh so that fxi+k ) = f(xi + kh). The difference, x, in the independent variable,
x, between the evaluation point, x (^) i + kh, and the expansion point, x (^) i , is equal to kh. Using x (^) i = a
as the expansion point and kh as x allows us to rewrite equation [2-5] as shown below.
3 3
3 2 2
2
kh dx
d f kh dx
d f kh dx
df f x kh f x xxi xxi x xi
i i [^2 -^12 ]
The next step is to use the notation that f(x (^) i + kh) = fi+k , and the following notation for the nth
derivative, evaluated at x = x i.
i (^) i xx i
n
n n i x x
i xx
i dx
d f f dx
d f f dx
df f
2 ' '' [2- 13 ]
With these notational changes, the Taylor series in equation [2-12] can be written as follows.
3 '''
2 ' ''
kh f
kh f (^) i k fi fikh fi i [2- 14 ]
Finite-difference expressions for various derivatives can be obtained by writing the Taylor series
shown above for different values of k, combining the results, and solving for the derivative. The
simplest example of this is to use only the series for k = 1.
3 '''
2 ' '' 1
h f
h f (^) i fi fih fi i [ 2 - 15 ]
We can rearrange this equation to solve for the first derivative, f’ (^) i ; recall that this is the first
derivative at the point x = xi.
and points forward (in the +x direc-tion) of that point. The expression in equation [2-19] is called
a backwards difference for similar reasons.
An expression for the first derivative that has a second-order error can be found by subtracting
equation [2-18] from equation [2-15]. When this is done, terms with even powers of h cancel
giving the following result.
5 '''''
3 ' ''' 1 ^ 1
h f
h f (^) i fi fih fi i [2- 20 ]
Solving this equation for the first derivative gives the following result.
1 1 2
4 '''''
2 ' 1 1 ''' Oh h
h f f f
h f h
f f f (^) i i i i i i i
^ [2- 21 ]
The finite-difference expression for the first derivative in equation [2-21] is called a central
difference. The point at which the derivative is evaluated, x i , is central to the two points (xi+1 and
xi-1) at which the function is evaluated. The central difference expression provides a higher order
(more accurate) expression for the first derivative as compared to the forward or backward
derivatives. There is only a small amount of extra work (a division by 2) in getting this more
accurate result. Because of their higher accuracy, central differences are usually preferred in
finite difference expressions.
Central difference expressions are not possible at the start of end of a boundary. It is possible to
get higher order finite difference expressions for such points by using more complex expressions.
For example, at the start of a region, x = x 0 , we can write the Taylor series in equation [2-14] for
the first two points in from the boundary, x 1 and x 2 , expanding around the boundary point, x 0.
3 ''' 0
2 '' 0
' 1 ^0 0
h f
h f f f h f [2- 22 ]
3 ''' 0
2 '' 0
' 2 ^0 0
h f
h f f f h f [2- 23 ]
These equations can be combined to eliminate the h 2 terms. To start, we multiply equation [2-22]
by 4 and subtract it from equation [2-23].
3 ''' 0
2 '' 0
' 0 0
3 ''' 0
2 '' 0
' 2 1 0 0
h f
h f f h f
h f
h f f f f h f
This equation can be simplified as follows
3 ''' 0
' 2 ^1 0 0
h f f f f h f [2-24]
When this equation is solved for the first derivative at the start of the region a second order
accurate expression is obtained.
2 ''' 0
' 2 1 0 0 Oh h
h f f f f h
f f f f
[2-25]
A similar equation can be found at the end of the region, x = xN , by obtaining the Taylor series
expansions about the point x = xN , for the values of f(x) at x = x N-1 and x = xN-2. This derivation
parallels the derivation used to obtain equation [2-25]. The result is shown below.
2 ' 2 1 ''' Oh h
h f f f f h
f f f f (^) N N N N N N N N
^ [2-26]
Equations [2-25] and [2-26] give second-order accurate expressions for the first derivative. The
expression in equation [2-25] is a forward difference; the one in equation [2-26] is a backwards
difference.
The evaluation of three expressions for the first derivative is shown in Table 2-1. These are (1)
the second-order, central-difference expression from equation [2-21], (2) the first-order, forward-
difference from equation [2-16], and (3) the second-order, forward-difference from equation [2-
25]. The first derivative is evaluated for f(x) = e x. For this function, the first derivative, df/dx = e x.
Since we know the exact value of the first derivative, we can calculate the error in the finite
difference results.
In Table 2-1, the results are computed for three different step sizes: h = 0.4, h = 0.2 and h = 0.1.
The table also shows the ratio of the error as the step size is changed. The next-to-last column
shows the ratio of the error for h = 0.4 to the error for h = 0.2. The final column shows the ratio of
the error for h = 0.2 to the error for h = 0.1.
Table 2-
Tests of Finite-Difference Formulae to Compute the First Derivative – f(x) = exp(x)
x f(x)
Exact f'(x)
h = .4 h = .2 h = .1 Error Ratios f’(x) Error f’(x) Error f’(x) Error (h=.4)/ (h=.2)
(h=.2)/ (h=.1)
Results using second-order central differences
0.6 1.8221 1. 0.7 2.0138 2.0138 2.0171 0. 0.8 2.2255 2.2255 2.2404 0.0149 2.2293 0.0037 4.
0.9 2.4596 2.4596 2.4760 0.0164 2.4637 0.0041 4. 1.0 2.7183 2.7183 2.7914 0.0731 2.7364 0.0182 2.7228 0.0045 4.02 4.
1.1 3.0042 3.0042 3.0242 0.0201 3.0092 0.0050 4. 1.2 3.3201 3.3201 3.3423 0.0222 3.3257 0.0055 4. 1.3 3.6693 3.6693 3.6754 0.
1.4 4.0552 4.
Results using first-order forward differences
0.6 1.8221 1.8221 2.2404 0.4183 2.0171 0.1950 1.9163 0.0942 2.15 2. 0.7 2.0138 2.0138 2.4760 0.4623 2.2293 0.2155 2.1179 0.1041 2.15 2.
0.8 2.2255 2.2255 2.7364 0.5109 2.4637 0.2382 2.3406 0.1151 2.15 2. 0.9 2.4596 2.4596 3.0242 0.5646 2.7228 0.2632 2.5868 0.1272 2.15 2. 1.0 2.7183 2.7183 3.3423 0.6240 3.0092 0.2909 2.8588 0.1406 2.15 2.
1.1 3.0042 3.0042 3.3257 0.3215 3.1595 0.1553 2. 1.2 3.3201 3.3201 3.6754 0.3553 3.4918 0.1717 2. 1.3 3.6693 3.6693 3.8590 0.
1.4 4.0552 4.
Results using second-order forward differences
0.6 1.8221 1.8221 1.6895 0.1327 1.7938 0.0283 1.8156 0.0066 4.69 4.
calculations were done with an excel spreadsheet which has about 15 significant figures. Figure
2-1 shows the effect of step size on error for a large range of step sizes.
For the large step sizes to the right of Figure 2-1, the plot of error versus step size appears to be
a straight line on this log-log plot. This is consistent with equation [2-17]. If we take logs of both
sides of that equation and solve for n, we get the following result.
log( ) log( )
log( ) log( )
log
log
2 1
2 1
1
2
1
2
h h
h
h
n
[2- 27 ]
Equation [2-27] shows that the order of the error is just the slope of a log(error) versus log(h) plot.
If we take the slope of the straight-line region on the right of Figure 2-1, we get a value of
approximately two for the slope, confirming the second order error for the central difference
expression that we are using here. However, we also see that as the step size reaches about
10 -5, the error starts to level off and then increase. At very small step sizes the numerator of the
finite-difference expression becomes zero on a computer and the error is just the exact value of
the derivative.
Figure 2-1. Effect of Step Size on Error
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E-
1.E+
1.E+
1.E-17 1.E-15 1.E-13 1.E-11 1.E-09 1.E-07 1.E-05 1.E-03 1.E-
Step Size
Error
Final Observations on Finite-Difference Expressions from Taylor Series
The notes above have focused on the general approach to the derivation of finite-difference
expressions using Taylor series. Such derivations lead to an expression for the truncation error.
That error is due to omitting the higher order terms in the Taylor series. We have characterized
that truncation error by the power or order of the step size in the first term that is truncated. The
truncation error is an important factor in the accuracy of the results. However, we also saw that
very small step sizes lead to roundoff errors that can be even larger than truncation errors.
The use of Taylor series to derive finite difference expressions can be extended to higher order
derivatives and expressions that are more complex, but have a higher order truncation error.
One expression that will be important for subsequent course work is the central-difference
expression for the second derivative. This can be found by adding equations [2-15] and [2-18].
4 ''''
2 '' 1 ^ 1
h f
h f (^) i fi fi fi i [2- 28 ]
We can solve this equation to obtain a finite-difference expression for the second derivative.
2
1 1
2 '''' 2
'' (^11) Oh
h
h f f f f h
f f f f (^) i i i i i i i i
^ [2- 29 ]
Although we have been deriving expressions here for ordinary derivatives, we will apply the same
expressions to partial derivatives. For example, the expression in equation [2-29] for the second
derivative could represent d 2 f/dx^2 or ∂^2 f/∂x^2.
The Taylor series we have been using here have considered x as the independent variable.
However, these expressions can be applied to any coordinate direction or time.
Although we have used Taylor series to derive the finite-difference expressions, they could also
be derived from interpolating polynomials. In this approach, one uses numerical methods for
developing polynomial approximations to functions, then takes the derivatives of the
approximating polynomials to approximate the derivatives of the functions. A finite-difference
expression with an n th^ order error that gives the value of any quantity should be able to represent
the given quantity exactly for an n th^ order polynomial.*
The expressions that we have considered are for constant step size. It is also possible to write
the Taylor series for variable step size and derive finite difference expressions with variable step
sizes. Such expressions have lower-order truncation error terms for the same amount of work in
computing the finite difference expression.
Although accuracy tells us that we should normally prefer central-difference expressions for
derivatives, we will see that for some convection terms special one-sided differences, known as
upwind differences will be used.
In solving differential equations by finite-difference methods, the differential equation is replaced
by its finite difference equivalent at each node. This gives a set of simultaneous algebraic
equations that are solved for the values of the dependent variable at each grid point.
- (^) If a second order polynomial is written as y = a + bx + cx (^2) ; its first derivative at a point x = x 0 is given by the
following equation: [dy/dx]x=x0 = b + 2cx 0. If we use the second-order central-difference expression in equation [2-21]to evaluate the first derivative, we get the same result as shown below:
0
0
2 0
2 0
2 0
2 0
2 0 0
2 0 0 0 0
( ) ( ) [ ( ) ( )]
0
b cx h
bh cxh
h
bh c x xh h cx xh h
h
a bx h c x h a bx h c x h
h
y x h y x h
dx
dy
x x
In the finite-element solution process, the values of φ 1 , φ 2 , 3 , and φ 4 in the element shown will
also be present in the equations for other elements. Consider, for example, the element to the
right of the element shown. This element will share a common boundary with the element shown.
That common boundary will include the unknowns that we have called φ 2 , and φ 3 for the element
shown. When we consider all the elements in the region, we have to assemble the equations for
each element in such a way that we identify the nodes that are common to more than one
element. For example, the value of φ 2 in the element shown here will be at the lower-left corner
of the element to the right and the upper-right corner of the element below the one shown. It will
also occur in the upper-left corner of the element that shares the corner, but no other boundary,
with the point for φ 2. The process of converting the element-by-element analysis to an analysis
of all elements (and all nodes) for the solution domain is called the assembly process for the
calculation.
There are a variety of ways to obtain the final finite-element equation. The original approach to
finite elements used a variational principle. In this approach some physical function, such as
energy, was set to a minimum. Such an approach is limited to linear differential equations. An
alternative approach, known as the method of weighted residuals, is applied to a simple problem
starting on page 17, operates as follows. The differential equation is written so that the right-hand
side is zero. The polynomial approximations, such as equation [2-32], are substituted into this
differential equation. The result is then integrated over the element, usually with some kind of
weighting function, and set equal to zero. Thus, instead of having the result be zero at every
point in the element, as it would be for the differential equation, the average result, integrated
over the element is set to zero. The result is a set of algebraic equations that involve the
unknown values such as φ 1 , φ 2 , φ 3 , and φ 4. The result for each element is assembled with the
results of neighboring elements to provide a set of simultaneous algebraic equations that can be
solved for all unknown values of φ.
Application of Finite-Differences to a Simple Differential Equation
To illustrate the use of finite differences and finite elements for solving differential equations we
will consider a simple one-dimensional problem shown in equation [2-33]. This equation
describes one-dimensional conduction heat transfer with a heat source proportional to the
temperature. If the equation for the heat source is bT, with b>0, the a 2 term in the equation below
equals b/k, where k is the thermal conductivity.
aT withT T atx andT T atx L dx
d T 0 A 0 B
2 2
2 [2-33]
To solve this problem using a finite-difference approach, we replace the differential equation by a
finite-difference analog. We will use a uniform, one-dimensional grid in the x direction. If we
have N+1 nodes on the grid, numbered from zero to N, the step size, h, will equal (L - 0)/N. At
any grid node, i, the finite-difference representation of the differential equation can be written as
follows, using equation [2-29] to convert the second derivative into a finite difference form.
2
1 ^ 1
i
i i i aT
h
T T T
[2-34]
The boundary conditions for the temperature provide the values of temperature at the two ends of
the grid: T 0 = TA and T N = TB. If we had a boundary condition expressed as a gradient or mixed
boundary condition, the boundary temperature would not be known, but we would have an
additional finite-difference equation to solve for the unknown boundary temperature. The system
of equations given by [2-34], with the known boundary conditions, can be written in detail as
follows. (The numbers in parentheses in the left margin represent an equation count.)
(1) (a 2 h 2 -2)T 1 + T 2 = -TA
(2) T 1 + (a 2 h 2 -2)T 2 + T 3 = 0
(3) T 2 + (a 2 h 2 -2)T 3 + T 4 = 0
[ (N – 8) additional similar equations with T 4 through T N-4 multiplied by (a 2 h 2 -2) go here. ]
(N-4) TN-4 + (a 2 h 2 -2)T (^) N-3 + TN-2 = 0
(N-3) TN-3 + (a 2 h 2 -2)T N-2 + TN-1 = 0
(N-2) TN-2 + (a 2 h 2 -2)T (^) N-1 = -TB
We see that we have to solve N-2 simultaneous linear algebraic equations for the unknown
temperatures T 1 to T (^) N-1. The set of equations that we have to solve is a simple one since each
equation is a relationship among only three temperatures. *^ None of the other temperatures on
the grid are unknowns here. We cannot solve the equation until we specify numerical values for
the boundary temperatures, T A and T B, and the length, L. We also need to specify a value for N.
This will give the value of h = (L – 0)/N.
Before specifying a, T A, TB, L, and N, we want to examine the exact solution of equation [2-33],
which is given by equation [2-35]. (You can substitute this equation into equation [2-33] to show
that the differential equation is satisfied. You can also set x = 0 and x = L in equation [2-35] to
show that this solution satisfies and the boundary conditions.)
sin( ) cos( ) sin( )
cos( ) ax T ax aL
T T aL T B^ A A
[2-35]
We can obtain the exact solution of equation [2-35] without specifying values for a, T A, TB, and L.
However, we have to specify these values (as well as a value of x) to compute a numerical value
for T. In the numerical solution process we cannot obtain an analytical form like equation [2-35];
instead, we must specify numerical values for the various parameters before we start to solve the
set of algebraic equations shown above. For this example, we choose a = 2, T A = 0, TB = 1, and
L = 1; we can then determine the solution for different values of N. Table 2-2 compares the exact
and numerical solution for N = 10.
The error in Table 2-2 is defined as the absolute value of the difference between the exact and
the numerical solution. This error is seen to vary over the solution domain. At the boundaries,
where the temperatures are known, the error is zero. As we get further from the boundaries, the
error increases, becoming a maximum at the midpoint of the region.
Table 2- Comparison of Exact and Numerical Solution for Equation [2-33] i x Exact Solution
Numerical Solution
Error
0 0 0 0 0
1 0.1 0.21849 0.21918 0. 2 0.2 0.42826 0.42960 0.
3 0.3 0.62097 0.62284 0. 4 0.4 0.78891 0.79115 0.
5 0.5 0.92541 0.92783 0. 6 0.6 1.02501 1.02739 0.
- (^) These equations are said to form a tridiagonal matrix. In a matrix format only the main diagonal and the
ones above it and below it have nonzero coefficients. A simple algorithm for solving such a set of equations is given in the appendix to these notes.
sin( )
cos( ) cos( ) cos( ) sin ( ) sin( )
2 2
aL
kaT T aL T aL T aL T aL aL
ka q
A B xL A B A
^ [2-40]
The heat flow at x = 0 is entering the region; the heat flow at x = L is leaving the region. The total
heat leaving the region, q x=L – q x=0 must be equal to the total heat generated. Since the heat
source term for the differential equation was equal to bT, the total heat generation for the region,
Q gen tot , must be equal to the integral of bT over the region. Using the exact solution in equation
[2-35], we find the total heat generated as follows.
ax T ax dx aL
T T aL Q bTdx b
L
A
B A
L
gen tot
0 0
sin( ) cos( ) sin( )
cos( ) [2-41]
Performing the indicated integration and evaluating the result at the specified upper and lower
limits gives the following result.
1 cos( ) sin( ) sin( )
cos( ) aL T aL aL
T T aL
a
b Q (^) A
B A gen tot [2-42]
We can simplify this slightly by using a common denominator of sin(aL) and using the same trig
identity, sin 2 x + cos 2 x = 1, used previously.
1 cos( )
sin( )
cos( ) cos( ) cos( ) sin ( ) sin( )
2 2
aL a aL
bT T
T T aL T aL T aL T aL a aL
b Q
B A
gentot B A B A A
[2-43]
The result for the total heat generated can be compared to the total heat flux for the region, found
by subtracting equation [2-39] from equation [2-40].
sin( )
cos( )
sin( )
cos( ) 0 aL
T T aL ka aL
kaT T aL q q
A B B A xL x
^ [2-44]
To make the comparison between the heat flux and the heat generated we need to use the
definition of a^2 =b/k in the paragraph before equation [2-35]. If we make this substitution in
equation [2-44] for the net heat flux and do some rearrangement, we obtain equation [2-45],
which confirms that the net heat outflow equals the heat generated, as shown in equation [2-43].
1 cos( )
sin( )
cos( ) cos( ) sin( )
0 aL a aL
bT T T T aL T T aL aL
ka q q
B A x L x A B B A
^ [2-45]
In the finite difference result we need to address the computation of the gradients at x = 0 and x =
L. To do this we will use the second-order expressions for the first derivative given in equations
[2-25] and [2-26].
0
h
T T T
dx
dT
x
[2-46]
h
T T T
dx
dT (^) N N N
xL
[2-47]
The exact and numerical values for the boundary temperature gradients are compared for both h
= 0.1 and h = 0.01 in Table 2-3. For both gradients, we observe the expected relationship
between error and step size for a second-order method: cutting the step size by a factor of ten
reduces the error by a factor of 100 (approximately).
Table 2- Comparison of Exact and Numerical Boundary Gradients in Solution to Equation [2-33] Location of Gradient
Exact Solution
Step Size
Numerical Solution
Error
x = 0 2.
h = 0.1 2.2357 0.
h = 0.01 2.1999 0.
x = L = 1 -0.
h = 0.1 -0.9332 0.
h = 0.01 -0.9155 0.
Application of Finite Elements to a Simple Differential Equation
Here we apply a finite-element method, known as the Galerkin method, *^ to the solution of
equation [2-33]. In applying finite elements to a one-dimensional problem, we divide the region
between x = 0 and x = L into N elements that are straight lines of length h. (For this one-
dimensional problem, the finite element points are the same as the finite difference points.) For
this example, we will use a simple linear polynomial approximation for the temperature. We give
this approximation the symbol, T ˆ^ , to distinguish it from the true temperature, T. If we substitute
our polynomial approximation into the differential equation, we can write a differential equation
based on the approximation functions. The Galerkin method that we will be using is one method
in a class of methods knows as the method of weighted residuals. In these methods one seeks
to set the integral of the approximate differential equation, times some weighting function, w (^) i (x), to
zero over the entire region. In our case, we would try to satisfy the following equation for some
number of weighting functions, equal to the number of elements.
aT dx i N dx
d T w x
L i^01 ,
0
2 2
2
^
[2-48]
In this approach, we are using an integral approximation. If we used the exact differential
equation, we would automatically satisfy equation [2-48] since the exact differential equation is
identically zero over the region.
In this one-dimensional case, the geometry is simple enough so that we do not have to use the
dimensionless ξ-η coordinate system. We have the same one-dimensional grid in this case that
we had in the finite difference example: a set of grid nodes starting with x = x 0 = 0 on the left side
and ending at x = xN = L on the right. (We will later consider the region between node i and node
i+1 as a linear element. We can write our approximate solution over the entire region in terms of
- (^) As noted above, the method of weighted residuals is just one approach to finite elements. This approach is
used here, rather than a variational approach, because variational approaches cannot be applied to the nonlinear equations that occur in fluid mechanics problems.
tridiagonal matrix equations to be solved for the unknown values of T (^) i. The desired set of
algebraic equations is obtained by substituting equation [2-49] into [2-48] and carrying out the
indicated differentiation and integration. The details of this work are presented below.
In finite elements it is necessary to distinguish between a numbering system for one element and
a global numbering system for all nodes in the region. In our case, the global system is a line
with nodes numbered from 0 to N. Our elements have only two nodes; we will use Roman
numerals, I and II, for the two ends of our element. We will continue to use the Arabic numbers 0
through N for our entire grid. Over a single element, then, we write the approximate temperature
equation [2-49] as follows.
I II II I
I II II I
II I I II II I x x x x x
x x T x x
x x T x T x T x T
ˆ( ) ( ) ( ) [2-51]
The second equation in [2-51] recognizes that there are four parts to the definition of the basis
functions in equation [2-50], only two of which are nonzero. Here we consider only the second
nonzero part of φI and the first nonzero part of φI+1. We will substitute the element equation in [2-
51] – instead of the general equation in [2-49] – into equation [2-48] to derive our desired result.
Equation [2-48] is a general equation for the method of weighted residuals. The Galerkin method
is one of the methods in the class of weighted residual methods. In this method, the basis
functions are used as the weighting functions. Since the individual elemental weighting function
is defined to be zero outside the element, we need only apply equation [2-48] to one element
when we are considering the weighting function for the element. For the Galerkin method then,
where the weighting function is the basis function, we can write equation [2-48] for each element
as follows.
aT dx i I II dx
xII dT
x i^
1
2 2
2
^
^ [2-52]
At this point, we apply integration by parts to the second derivative term in equation [2-52]. This
does two things for us. It gets rid of the second derivative, which allows us to use linear
polynomials. If we did not do so, taking the second derivative of our linear polynomial would give
us zero. It also identifies the boundary gradient as a separate term in the analysis. (In a
multidimensional problem, we would use Green’s Theorem; this has a similar effect, in multiple
dimensions, to integration by parts in one dimension.) The general formula for integration by
parts is shown below.
b
a
b
a
b
a
udv uv vdu [2-53]
We have to manipulate the second derivative term in equation [2-50] to get it into the form shown
in [2-53]. We do this by noting that the second derivative is just the derivative of the first
derivative and we can algebraically cancel dx terms as shown below.
2
1 1
2
(^2) x
x i
x
x i
x
x i^ dx
dT dx d dx
dT
dx
d dx dx
d T II
I
II
[2-54]
The final integral in equation [2-54] has the form of the integration by parts formula in equation [2-
53]. We can identify u = φi , and v =^ dT^ ˆ^ dx .Using these definitions of u and v in equation [2-53]
gives the following result.
dx dx
dT
dx
d
dx
dT dT dx
d
dx
dT dx dx
d T II
I
II
I
II
I
II
I
II
I
x
x
i
x
x
i
x
x
i
x
x
i
x x (^) i^
2
2
[2-55]
With this result, we can rewrite equation [2-52] as follows.
dx a T dx i I II dx
dT
dx
d
dx
dT II
I
II
I
II
I
x
x i
x
x
i
x
x
i^0 ,
2 (^)
[2-56]
We need to evaluate this equation using the linear basis functions in equation [2-50] that we have
selected for this analysis. The basis functions have the following first derivatives for the element
under consideration. Once the derivatives of the basis functions are known we can compute the
derivatives of the approximate temperature. In taking the derivatives of the approximate
temperature polynomial, we consider the nodal values, T I and T II, to be constants.
II I
II I
II I
II
II I
II I II
I I
II I
II
II I
I
x x
T T
x x
T
x x
T
dx
d T dx
d T dx
dT
dx x x
d and dx x x
d
[2-57]
We can substitute the temperature polynomial and the basis functions, and their derivatives, into
equation [2-56]. We will show the details for φ = φI, on a term-by-term basis, starting with the first
term in equation [2-56]. The basis function, φI, is zero at the upper limit of this evaluation, x = x (^) II)
so we only have the lower limit where φI = 1. Rather than substituting the interpolation
polynomial for the approximate temperature gradient, we replace this term by the actual gradient.
We can then handle boundary conditions that use this gradient. If we do not have a gradient
boundary condition, we can use the resulting equations to compute the gradient.
1 1
II I xx x x
II
II I xx
II
x
II I x
II
x
x
I dx
dT
dx
dT
x x
x x
dx
dT
x x
x x
dx
dT
x x
x x
dx
dT
II
II
I
II
I
[2-58]
The middle integral with two derivative terms becomes the simple integral of a constant
II I
II I
x
II I x
x II I
x II I
II I
II I
x
x
I
x x
T T
x x x
T T
dx x x
T T
x x
dx dx
dT
dx
d
II
I
II
I
II
I (^)
^ 2 ( )
[2-59]
The final term in equation [2-56] requires the most work for integration. The last step in this
integration is left as an exercise for the interested reader.
