Advanced Transport Phenomena Module 8 Lecture 36 Outline of Numerical Solution Methods

Dr. R. Nagarajan Professor Dept of Chemical Engineering IIT Madras

1

Outline of Numerical Solution Methods 2

Method of Finite Differences (MFD) for the Numerical Solution of Partial Differential (Field) Equations (PDEs) and Ancillary Boundary Conditions (BCs) 3

MFD  For direct solutions of coupled, nonlinear PDEs governing the fields of momentum density, energy density, and species mass densities of interest  Often for domains (and flow inlet/outlet locations) of complicated geometry 4

MFD  Analysis of transport processes beset by the difficulty in writing down the governing PDEs and BCs  Situation improved markedly owing to:  availability of high-speed/large-memory digital computers, and electronic  development of efficient

numerical

methods based on converting the governing PDEs and BCs to a coupled system of

algebraic

equations amenable to digital computation, often using matrix methods.

 One class of such methods: "method of finite differences" (MFD) 5

MFD  Basic Idea:  Give up the notion of attempting to determine the dependent field quantity

u(x,t ) —

say, at "every" point within the domain

D

of independent variables.

 Accept instead the values of

u

at discrete "node" points within the domain; i.e., seek a finite number of values: u i,j,k,n , where the indices denote discrete spatial locations in

i,j,k (= D,

0, 1, 2, ...) and n (=0,1,2,...) denotes discrete locations in time domain 6

MFD 

Question:

How can the numerical values of u i,j,k,n be found such that the original PDE + BC + IC can be acceptably approximated?



Steps:

1.

Write the dimensionless PDE + BC + IC for the dependent variable(s) u in the most appropriate orthogonal coordinate system. The field equation is generally of the form: 

u t

N

 t 

N

 7

MFD where, in most problems arising in transport theory,

N

(u) will be a nonlinear (or quasi-linear) second-order (in space variables)

differential

operator Boundary conditions (BC) are frequently of the form

a

B  

v

B



b u

B B



c

B



0

8

MFD where and

v

is a length coordinate normal to the boundary B

a, b, c

are coefficients independent of

u

itself.

 BC includes the commonly encountered special cases

a =

0

(u

B specified) and

b =

0 (specified "flux ” ).

 Initial conditions (usually u(

x

, 0)) must also be specified.

9

MFD 2. Lay down a "mesh" in both time and physical (transformed and/or normalized) space variables and accept as a solution the discrete values u i,j,k,n , of the field density at the "node" points (provided the conditions implied by Steps 3, 4, and 6 below are satisfied).

10

MFD 3. Replace the PDE by an appropriate "difference analog" which will provide

algebraic

interrelations between the u i,j,k,n.

.

 Difference analogs (which can be obtained by truncating Taylor series, integration formulae, or variational methods) are not unique; i.e., they exhibit important differences with respect to the local magnitude (truncation) and/or growth (or damping) of inevitable random numerical errors ("stability").

11

MFD  According to the so-called "equivalence theorem" of PDE-numerical analysis, a numerical method which locally approximates the PDE derivatives with mesh refinement ("consistency"), and which is free from the exponential amplification of random numerical errors (i.e., a globally "stable" method) will converge on the solution to the originally posed PDE-BVP.

12

MFD 4. Replace the BC by a difference analog; in all, sufficient algebraic equations must now be available to calculate all of the required u i,j,k,n , frequently by "marching" in t-space.

 The difference analog of the BC may be complicated if the boundary B of the domain

D

is not a simple shape (e.g., not a coordinate surface).

13

MFD 5. Solve the resulting algebraic system of equations for the u i,j,k,n using calculations.

matrix methods/notation to organize the  Solution procedure simplest if the equations are linear and the analog expressions for

N

(u) at most link variables at three adjacent mesh points.

14

MFD  Solution procedure trivial if an "explicit method" is used since each u i,j,k,n+1 , is calculated directly from the "known" (previously determined) values u i.j.k,n ;  however, such methods are usually the mesh is sufficiently fine.

unstable

unless  "Implicit" methods, i.e., FD algorithms for which one must solve implicit algebraic equations for the nodal unknowns, are stable (with respect to the growth of cumulative errors) for any mesh size, even if the source term in the DE is extremely sensitive to the dependent variable (so-called "stiff" equation).

15

MFD 6. Repeat the above process with a finer or coarser mesh; accept the result only if the corresponding u-values (and/or their integrated consequences) are insensitive to mesh size.

 Mesh size cannot be so small that inevitable roundoff errors cause large relative errors in computing first and second differences of values at adjacent mesh points  Ordinarily computation time (and cost) increase rapidly with mesh fineness, and cost limitations set in before roundoff errors do.

16

MFD

Comments

 There are now many excellent books available on numerical methods for solving PDEs (the selection might well be made based on the relevance of the illustrations to the particular interests of the reader). An excellent introductory reference for the MFD is still Von Rosenberg, D. (1969).

17

MFD  For many commonly encountered PDE-BVPs, pre programmed "subroutines" are available on call. An example is Harwell Subroutine DP F LA (Crank-Nicolson method) for quasi-linear "parabolic" problems of the form: u t =a(u, …)u xx + bu x + cu + d

,

with linear two-point boundary data. (Here, the subscripts

t,

xx, and x denote

derivatives

with respect to the indicated variable).

18

MFD  For problems involving irregularly shaped domains and/or steep local gradients, an alternative "discretization" procedure, called the method of finite

elements

(MFE) is also in widespread use; however, the ability to rapidly compute and use "body-fitted" orthogonal coordinates has also permitted the MFD to be effectively applied to irregularly shaped domains (see, e.g., Thompson (1982)).

et al.

19

MFD  More generally, one can use MFD to solve for the fields of interest

and

“simultaneously" an optimal (nearly orthogonal) coordinate grid, "adaptively" distributing the grid points in accordance with a criterion that accounts for the locations of maximum gradients in the dependent variable (see, e.g., Shyy, Tong, and Correa (1985)).

20

MFD  Regions of a flow field in which convection dominates diffusion (momentum, energy, and/or mass) are approximately described by PDEs

without

second order derivative terms (called locally "hyperbolic" PDEs). Within such regions there are local directions (called characteristic directions in the space of independent variables)

along

which the dependent variables can be shown to satisfy ODEs.

21

MFD  This property is numerically exploited in the "method-of characteristics" (MC) (see, e.g., Abbott (1966)). In effect, of all possible local coordinate systems, an optimal set (called "characteristic coordinates") is calculated and then used to advance the solution. When applicable, the MC is usually preferable to the MFD, especially for highly nonlinear problems. However, usually boundary layers (within which diffusion

cannot

be neglected) prevent global applicability of a MC.

22

Method of Finite Elements (MFE) for the Numerical Solution of PDEs on Domains of Complicated Shape

23

MFE 

Steps (Steady-State Illustration)

1.If available, first write a valid

variational principle

for the governing PDE and BCs on the domain of interest.

 For linear problems, the corresponding Lagrange functional (whose extremum is to be found) usually involves an integrand (Lagrange "density") at most quadratic in the dependent variable

u

and its space derivatives.

2. Subdivide the domain into convenient "finite elements" (e.g., tetrahedra for three dimensions) with dependent variables u ijk at the

nodes

as the basic unknowns.

 Use some interpolation

everywhere else.

procedure to define

u

24

MFE 3. Use the variational principle, together with the fact that each finite element contributes additively to the Lagrange integral, to derive a set of the nodal values u ijk

.

algebraic

equations interrelating  This leads to sparse matrix equations for the u ijk 4. The solution to these algebraic equations defines the MFE-discretized solution to the original continuum boundary-value problem 25

MFE

Advantages (Relative to Finite Differences)

 Ease with which shapes and sizes of the finite elements (defined by the positions of the node points) can be chosen.

 This allows elliptic PDEs within complex geometries to be handled, allowing for steep local gradients.

 Ease of handling

boundary conditions via

the basic variational integral, and less truncation error for gradient conditions.

26

MFE

Remarks

 Often a variational principle is

not

available, in which case a "weighted PDE-residual" procedure is used to generate similar algebraic equations (i.e., a weighted local error of the trial solution is "minimized" within each finite element).

 In such cases, the trial function nodal values u ijk and ("local coordinates");

u

is defined

interpolation functions via

the of position  the domain residual within each finite element is made orthogonal to each interpolation formula.

27

MFE  When the mesh, domain, geometry, and trial functions are particularly simple, MFE becomes equivalent to the method of finite differences (MFD)  As in finite-difference methods, nonlinearities are handled by linearization at each stage of the matrix computations, using convergent

iterative

methods.

28

Method of Weighted Residuals (MWR) for the Approximate Solution of Partial Differential (Field) Equations and Ancillary Conditions (BCs, ICs)

29

MWR  Well suited to obtain approximate semi-analytical or efficient numerical solutions to transport problems in which:  Only integrated properties are needed (as opposed to local accuracy);  The solution is not expected to exhibit steep local spatial gradients;  The spatial domain possesses a high degree of symmetry.

30

MWR 

Strategy

:  Consider the problem of finding PDE: u t =

N

(u) (where

u(x,t)

satisfying the

N

(u) is typically a nonlinear, second-order spatial differential operator) on some spatial domain,

V,

with specified initial conditions (ICs) and boundary conditions (BCs).

 Introduce an approximate solution u(x, satisfies the BCs, and contains

N specified t)

which "basis-," trial-, or shape-functions of the spatial coordinates (defined over the entire domain

V)

and

N

x

as yet

undetermined

functions of time.

31

MWR  By satisfying the PDE (expressing local conservation) as accurately as possible (globally, not necessarily locally), it is possible to obtain allowing the

N N

ODEs (and their ICs) previously undetermined functions of time to be found, where

N

may be as small as 1.

32

MWR  Satisfying global conservation conditions is usually discussed in terms of certain (weighted-) integrals of the local “domain-(interior) residual” or local error:

N

  

u t '

33

MWR [ ε v would clearly vanish everywhere if the approximating function u(x,

t)

were indeed exact.]  A subclass of such methods (called the "least-squares" method) simply minimizes the integral.

v   v 2

dV

34

MWR  More generally, one can force any

N

"weighted-domain residuals" to vanish; that is, v 

w

 i v     where the weighting functions w i (

x

) can be smooth functions of position, unrelated to the trial functions, chosen to emphasize locations where accuracy is felt to be especially important.

35

MWR  Alternatively, the weighting functions can be the same as the trial spatial functions themselves, and members of a complete set of functions in the domain

V

(Galerkin method), or even Dirac delta functions centered about preselected points in the domain [equivalent to simply demanding that the domain residual (interior error ε v ) vanish only at these selected ("collocation"-) points x i

(i =

1, 2, ... , avoiding time-consuming quadratures].

N),

thereby 36

MWR While such WR-methods can be computationally efficient (especially when considerable knowledge of the spatial dependence and symmetry of the solutions is "built in" to the trial spatial functions), in general they suffer from a lack of information about their convergence properties (e.g., when

N

increases) and, hence, accuracy. Usually it is wise to demonstrate at least one of the following properties: 37

MWR  Additional terms (incrementing

N)

lead to only a small change in the calculated quantity of interest (QOI).

 Alternative but qualitatively similar trial functions (spatial dependencies), give only a small change in the QOI.

 Your method, when applied to a well-known special case (for which an exact solution may be available) gives acceptable accuracy for a similar QOI.

38

MWR  A pedagogically interesting simple example is provided by the transient linear (heat-) diffusion. In this case, the ultra-simple approximating function:

T

w



T

w



T

0

  

x

h

1 0 x

h

x

 

h     39

MWR  (involving only

one

undetermined function of time; the thermal penetration depth,

viz.,

δ h (t), which vanishes at

t =

0),   combined with a macroscopic energy conservation condition over the "entire" domain 0 ≤x ≤ δ h (t), leads to an ODE for δ h

= 2 (αt) 1/2 .

δ h (t) with the analytic solution:  While local relative errors in the temperature profile are excessive, the corresponding wall

heat fluxes

are found to be (under-)predicted by only 11.4 percent (cf. exact solution)!

40

MWR

Comments

 Theoretical basis for systematically weighting functions w i {x} is as follows: choosing the  If the weighting functions w i (x) are members of a "complete" set of functions on the domain V, then the only function w i

ε v (u) ε v (u)

for which vanishes is, in fact,

ε v every (u) =

0.

weighted integral of  In any practical realization of MWR, only the first few weighted residuals of

ε v

are forced to vanish; hence the associated u(x,

t)

will be an

approximate

solution of the PDE-BVP.

41

MWR  There are deep interrelations between these variants, usually demonstrated for simple linear problems.

 The Galerkin method [1915] has been investigated most extensively, and has the strongest theoretical foundation (since it has the property that, if the exact solution is contained in the trial (basis-)functions, the Galerkin method will find it [Finlayson (1972)]) 42

MWR  Also used "locally" in implementing the method of finite elements (MFE) when:  (a) a variational principle does not exist for the PDE BVP under investigation (the usual situation for nonlinear transport problems), and  (b) MFE is desirable because of the complexity of the overall domain

V

and/or the expected presence of strong local spatial nonuniformities.

43

MWR  In such cases, instead of being defined over the entire domain

V,

the trial (shape-) functions are defined only within each finite

element

(into which the entire domain has been judiciously subdivided) and can be considered to vanish elsewhere.

 Thus, the Galerkin method can be used to approximately satisfy the PDE "over" each finite element, and the "local solution" can be explicitly rewritten in terms of relevant adjacent nodal values

u

i.j,k,n

of the dependent variable.

44

MFD/ MFE/ MWR: Exercises

 For what sort of energy/mass/momentum transport problems would recourse be made to numerical methods (e.g., MFD)?

 What is meant by (a) truncation error? (b) round-off error? (c) convergence? (d) "stability"?

How do these considerations influence the choice of algorithm and mesh spacing? What other factors govern the choice and uniformity of mesh spacing?

45

MFD/ MFE/ MWR: Exercises

 How are difference analog equations obtained? Outline

one

such method. Is there more than one "correct" difference analog of a particular PDE + BC + IC? Write two possible difference analogs for the transient heat conduction equation in one dimension, one being "explicit" and one being "implicit," and discuss their relative merits. How is the Schmidt (graphical) method related to the explicit analog given above?

46

MFD/ MFE/ MWR: Exercises

 Why are matrix methods convenient for organizing and carrying out numerical calculations of the discretized (algebraic) equations? Why is it desirable to linearize

(via

estimation methods) nonlinear terms at each stage of the calculation?

47

MFD/ MFE/ MWR: Exercises  Should the "primitive" differential equations be attacked directly, or is it desirable to first introduce: (a) non dimensionalized, normalized variables, (b) "stretched" (distorted) coordinates to magnify regions in which gradients are expected to be large? If strategy (b) is used, does this influence the nature of the differential equation to be discretized?

48

MFD/ MFE/ MWR: Exercises  In a one-dimensional transient problem (say possible to discretize only one variable (say

y) y, t)

is it and leave the other variable

(t)

continuous? What methods could be used to solve the resulting numerical problem?

 What is the basic difference between problems that are mathematically "parabolic" (containing a time-like or "one-way" variable) and "elliptic" (containing all "two-way" variables)? How does this influence the numerical procedures adopted?

49

MFD/ MFE/ MWR: Exercises  Suppose one obtains a discretized solution (so that all relevant dependent variables are known at the "mesh points"). How elsewhere?

could

(Optional:

these variables be estimated What is a "spline" function, and how could it be applied for interpolation purposes?)  An alternative to the method of finite difference (MFD) which retains most of its attributes but offers advantages for complicated-shaped domains and boundary conditions, is the "finite-element method". Compare MFE with MFD and indicate the essential advantages of MFE over the method of finite differences.

50