Numerical Solution Methods for PDEs in Advanced Transport Phenomena

1. Advanced Transport Phenomena Module 8 Lecture 36 Outline of Numerical Solution Methods Dr. R. Nagarajan Professor Dept of Chemical Engineering IIT Madras

2. Outline of Numerical Solution Methods

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

4. 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

5. 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 electronic digital computers, and • development of efficient numerical methods based on converting the govern­ing 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)

6. 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: ui,j,k,n,where the indices i,j,k (= 0, 1, 2, ...) denote discrete spatial locations in D, and n (=0,1,2,...) denotes discrete locations in time domain

7. MFD • Question:How can the numerical values of ui,j,k,nbe found such that the original PDE + BC + IC can be acceptably approximated? • Steps: • 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:

8. 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

9. MFD where v is a length coordinate normal to the boundary B and a, b, c are coefficients independent of u itself. • BC includes the commonly encountered special cases a = 0 (uBspecified) and b = 0 (specified "flux”). • Initial conditions (usually u(x, 0)) must also be specified.

10. 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 ui,j,k,n, of the field density at the "node" points (provided the conditions implied by Steps 3, 4, and 6 below are satisfied).

11. MFD 3. Replace the PDE by an appropriate "difference analog" which will provide algebraic interrelations between the ui,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").

12. 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.

13. MFD 4. Replace the BC by a difference analog; in all, sufficient algebraic equations must now be available to calculate all of the required ui,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).

14. MFD 5. Solve the resulting algebraic system of equations for the ui,j,k,nusing matrix methods/notation to organize the calculations. • Solution procedure simplest if the equations are linear and the analog expressions for N(u)at most link variables at three adjacent mesh points.

15. MFD • Solution procedure trivial if an "explicit method" is used since each ui,j,k,n+1, is calculated directly from the "known" (previously determined) values ui.j.k,n; • however, such methods are usually unstable unless the mesh is sufficiently fine. • "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).

16. 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 insensi­tive to mesh size. • Mesh size cannot be so small that inevitable roundoff errors cause large relative errors in computing first and second dif­ferences 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.

17. 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).

18. MFD • For many commonly encountered PDE-BVPs, pre-programmed "subroutines" are available on call. An example is Harwell Subroutine DPFLA (Crank­-Nicolson method) for quasi-linear "parabolic" problems of the form: ut =a(u,…)uxx+ bux + cu + d, with linear two-point boundary data. (Here, the subscripts t, xx, and x denote derivatives with respect to the indicated variable).

19. 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 et al. (1982)).

20. MFD • More generally, one can use MFD to “simultaneously" solve for the fields of interest and 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)).

21. 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.

22. 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.

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

24. 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 uijkat the nodes as the basic unknowns. • Use some interpolation procedure to define u everywhere else.

25. 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 algebraic equations interrelating the nodal values uijk. • This leads to sparse matrix equations for the uijk 4. The solution to these algebraic equations defines the MFE-discretized solution to the original continuum boundary-value problem

26. 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.

27. 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 u is defined via the nodal values uijk and interpolation functions of position ("local coordinates"); • the domain residual within each finite element is made orthogonal to each interpolation formula.

28. 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. 

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

30. 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.

31. MWR • Strategy: • Consider the problem of finding u(x,t) satisfying the PDE: ut = N(u) (where 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, t) which satisfies the BCs, and contains N specified "basis-," trial-, or shape-functions of the spatial coordinates x (defined over the entire domain V) and N as yet undetermined functions of time.

32. MWR • By satisfying the PDE (expressing local conservation) as accurately as possible (globally, not necessarily locally), it is possible to obtain N ODEs (and their ICs) allowing the N previously undetermined functions of time to be found, where N may be as small as 1.

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

34. MWR [εvwould 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.

35. MWR • More generally, one can force any N "weighted-domain residuals" to vanish; that is, where the weighting functions wi(x) can be smooth functions of position, unrelated to the trial functions, chosen to emphasize locations where accuracy is felt to be especially important.

36. 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 xi(i = 1, 2, ... , N), thereby avoiding time-consuming quadratures].

37. 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 infor­mation 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:

38. 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.

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

40. MWR • (involving only one undetermined function of time; viz., the thermal penetration depth, δh(t), which vanishes at t = 0), • combined with a macroscopic energy con­servation condition over the "entire" domain 0 ≤x ≤ δh(t), • leads to an ODE for δh(t)with the analytic solution: δh = 2(αt)1/2. • 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)!

41. MWR Comments • Theoretical basis for systematically choosing the weighting functions wi{x} is as follows: • If the weighting functions wi(x) are members of a "complete" set of functions on the domain V, then the only function εv(u) for which every weighted integral of wiεv(u) vanishes is, in fact, εv(u) = 0. • 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.

42. 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)])

43. 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.

44. 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 ui.j,k,n of the dependent variable.

45. 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?

46. 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?

47. 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?

48. 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?

49. MFD/ MFE/ MWR: Exercises • In a one-dimensional transient problem (say y, t) is it possible to discretize only one variable (say y) 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?

50. MFD/ MFE/ MWR: Exercises • Suppose one obtains a discretized solution (so that all relevant dependent variables are known at the "mesh points"). How could these variables be estimated elsewhere? (Optional: 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.