ME 608
Numerical Methods for Heat, Mass and Momentum Transfer Jayath Jayathii Y. Murthy Murthy Professor, School of Mechanical Engineering Purdue University
[email protected]
Spring 2006
Lecture 1: Introduction to ME 608 Conservation Equations
Outline of Lecture
Course organization
Introduction to CFD
Conservation equations, general scalar transport equation
Conservation form
Motivation
Huge variety of industrial flows: •Rotating machinery •Compressible/incompressible aerodynamics •Manifolds, piping •Extrusion, mixing •Reacting flows, combustion …. Impossible to solve Navier-Stokes equations analytically for these applications!
History
Earliest Earliest “CFD” “CFD” work by L.F. Richard Richardson son (1910) » Use Used hu huma man n com compu pute ters rs » Iter Iterat ativ ive e solu soluti tion ons s of of Lapl Laplac ace’ e’s s eq eqn n using finite-difference methods, flow over cylinder etc. » Erro Errorr esti estima mate tes, s, ext extra rapo pola lati tion on to to zero error
“So far I have paid piece rates for the operation (Laplacian) of about n/18 pence per coordinate point, n being the number number of digits digits … one of the the quickes quickestt boys averaged 2000 operations (Laplacian) per week for numbers numbers of 3 digits, those done wrong being discounted …” Richardson, 1910
Lewis F. Richardson (1881-1953)
Also researched mathematical models for causes of war : Generalized Foreign Politics (1939) Arms and Insecurity(1949) Statistics of Deadly Quarrels (1950)
History
Relaxation methods (1920’s-50’s) Landmark Landmark paper paper by Courant, Courant, Friedric Friedrichs hs and Lewy for hyperbol hyperbolic ic equation equations s (1928) (1928) Von Neumann stability criteria for parabolic problems (1950) Harlow Harlow and Fromm Fromm (1963) (1963) computed computed unsteady unsteady vortex street using a digital computer. They published a Scientific American article (1965) which ignited interest in modern CFD and the idea of computer experiments Boundary-layer codes developed in the 19601970’s (GENMIX (GENMIX by Patank Patankar ar and Spalding Spalding in 1972 for eg.) Solution techniques for incompressible flows published through the 1970’s (SIMPLE family of algorithms algorithms by by Patankar Patankar and Spalding Spalding for for eg.) Jameson computed Euler flow over complete complete aircraft (1981) Unstructured mesh methods developed in 1990’s
John von Neumann (1903-1957)
Richard Courant (1888-1972)
Conservation Equations
Nearly all physical processes of interest to us are governed by conservation equations » Mass, Mass, moment momentum um ene energy rgy conser conservat vation ion
Written Written in terms of specific specific quantities quantities (per unit unit mass basis) basis) » Momentum per unit mass (velocity) » Ener Energy gy per un unit it mass ass e
Consider a specific quantity φ » Cou Could ld be momen momentum tum per per unit unit mass, mass, ener energy gy per per unit unit mass.. mass..
Write conservation statement for φ for control volume of size
∆x x ∆y x ∆z
Conservation Equations (cont’d) Accumulation of φ in control volume over time step ∆t = Net influx of φ into control volume - Ne Nett eff efflu lux x of of φ out of control volume + Net generation of φ inside control volume
Conservation Equations (cont’d) Accumulation:
Generation:
Influx and Efflux:
Diffusion and Convection Fluxes Diffusion Flux Convection Flux Net flux
Velocity Vector Diffusion coefficient
Γ
Combining…
Taking limit as ∆x, ∆y, ∆z -> 0
General Scalar Transport Equation
Or, in vector form:
Conservation Form Consider steady state. The conservation form of the scalar transport equation is:
Non-Conservation Form
Finite volume methods always start with wit h the conservation form
General Scalar Transport Equation
Storage
Recall: say)
Convection
Diffusion
Generation
mass φ is a specific quantity (energy per unit mass
V : velocity vector
Γ: Diffusion coefficient ρ: density S: Source term (Generation per unit volume W/m3)
Continuity Equation ∂ ρ + ∇ ⋅ ( ρ V ) = 0 ∂t Here,
φ= 1 Γ= 0 S = 0
Energy Equation
h = sensible enthalpy per unit mass, J/kg k = thermal conductivity Sh = energy generation W/m3
Note: h in in convection convection and storage storage terms T in diffusion terms How to cast in the form of the general scalar transport equation?
Energy Equation (cont’d) Equation of State
Substitute to Find
Here,
φ= h Γ= k/C p S = S h
Momentum Equation X-Momentum Equation
Here,
φ= u
τ ij
⎛ ∂ui ∂u j ⎞ = µ ⎜⎜ + ⎟⎟ ⎝ ∂ x j ∂xi ⎠
Γ= µ
∂p S = S u - ∂ x
S is good good “dumpin “dumping g ground ground”” for everything that doesn’t fit into the other terms
Species Transport Equation
Yi = kg of specie i /kg of mixture
Γi = diffusion coefficient of i in mixture i Ri = reaction source
Closure
In this lecture we » Develo Developed ped the the proce procedur dure e for devel developi oping ng the gove governi rning ng equation for the transport of a scalar φ » Re Reco cogn gniz ized ed the the com commo mona nalility ty of tran transp spor ortt of of – Mass, Mass, momen momentum tum,, energy, energy, spec species ies
» Casting Casting a allll these these different different equations equations into this single single form form is very useful » Can devi devise se a single single meth method od to solve solve this this class class of of governi governing ng equation
Lecture 2: The General Scalar Transport Equation Overview of Numerical Methods
Last time…
Wrote conservation statement for a control volume
Derived a general scalar transport equation
Discovered that all transport processes commonalities » » » »
Storage Diffusion Convection Generation
This time…
Examine important classes of partial differential equations and understand their behavior See how this knowledge applies to the general scalar transport equation Start a general overview of the main elements of all numerical methods
General Scalar Transport Equation
Storage
Recall: say)
Convection
Diffusion
Generation
mass φ is a specific quantity (energy per unit mass
V : velocity vector
Γ: Diffusion coefficient ρ: density S: Source term (Generation per unit volume W/m3)
Classification of PDEs Consider the second-order second-order partial differential differential equation for φ (x,y):
Coeffi Coef fici cien ents ts a, a,b, b,c, c,d, d,e, e,ff are are linea linearr -- no nott fu func ncti tion ons s of φ, but can be functions of (x,y) Discriminant D<0
Elliptic PDE
D=0
Parabolic PDE
D>0
Hyperbolic PDE
Elliptic PDEs Consider 1-D heat conduction in a plane wall with constant thermal conductivity T o T L
Boundary conditions
Solution:
Elliptic PDE’s
T o • T(x) is influenced by both
boundaries • In the the absence absence of source source terms, terms, T(x) is bounded by the values on both boundaries •Can we devise numerical schemes which preserve these properties?
T L
Parabolic PDEs Consider 1D unsteady conduction in a slab with constant properties:
Boundary and initial conditions
T 0
Solution:
T i
T 0
Para Parabo bolic lic PDEs PDEs (con (cont’ t’d) d) T 0
T i
T 0
• The solution at T(x,t) is influenced by the boundaries, just just as with elliptic elliptic PDEs •We need only only initial condtions T(x,0). We do not need future conditions •Initial conditions only affect future conditions, not past conditions • Initial conditions conditions affect affect all spatial points points in the future • A steady steady state state is reache reached d as t->∞. In this limit we recover the elliptic PDE. •In the absence of source terms, the temperature is bounded by initial and boundary conditions •Marching solutions are possible
Hyperbolic PDEs Consider the convection of a step change in temperature:
Initial and boundary conditions
Solution:
Hype Hyperb rbol olic ic PDE PDEss (con (cont’ t’d) d)
Hype Hyperb rbol olic ic PDE PDEss (con (cont’ t’d) d) • Upstream Upstream conditions conditions can potentially affect the solution at a point x; downstream conditions do not • Inlet Inlet condition conditions s propagate propagate at a finite speed U •Inlet condition is not felt at location x until a time x/U
Relation to Scalar Transport Equation
• Contains Contains all three three canonic canonical al PDE terms terms • If Re is low low and situation situation is steady, steady, we get an elliptic equation • If diffusion coefficient is zero , we get a hyperbolic hyperbolic equation • If Re is low low and situation situation is unsteady, unsteady, we get a parabolic equation • For mixed mixed regimes, regimes, we get mixed mixed behavio behaviorr
Components of CFD Solution
Geometry creation
Domain Domain discre discretizatio tization n (mesh (mesh genera generation) tion)
Discretiza Discretization tion of governing governing equations equations
Solution of discrete equations; accounting for nonlinearities linearities and inter-equa inter-equation tion coupling coupling Visualization and post-processing
Solution Process
Analytical solution gives us φ(x,y,z,t). Numerical solution gives us φ only at discrete grid points. The process of converting the governing partial differential equation into discrete algebraic equations is call discretization. Discretization involves
» Disc Discre reti tiza zatio tion n of spac space e usin using g mesh mesh gene genera ratio tion n » Discre Discretiz tizati ation on of govern governing ing eq equat uation ions s tto o yield yield sets sets of algebraic equations
Mesh Types
Regular and body-fitted meshes
Stair-stepped representation of complex geometry
Mesh types (cont’d)
Blockstructured meshes
Unstructured meshes
Mesh Types Cell shapes
Nonconformal mesh Hybrid mesh
Mesh Terminology
• Node-based
finite volume scheme : φ stored at vertex
• Cell-based finite volume scheme : φ stored at cell centroid
Overview of Finite Difference Method Consider diffusion equation:
Step 1: Discre Discretize tize domain domain using using
a mesh. Unknowns are located at nodes
Step 2: Expand φ in Taylor series about point 2
Subtracting equations yields
Finite Difference Method (cont’d)
Step 3: Adding equations yields
Drop truncated terms:
Step 4: Evaluate source term at point 2:
Second order truncation error
Finite Difference Method (cont’d)
Step 5: Assemble discrete equation
Comments » We can can writ write e one one such such equat equation ion for each each grid grid point point » Bounda Boundary ry cond conditi itions ons give give us bou bounda ndary ry valu values es » Seco Second nd-o -ord rder er accu accura rate te » Nee Need d to find find a way way to solve solve coup couple le algeb algebrai raic c equati equation on set set
Overview of Finite Volume Method Consider the diffusion equation:
Step 1: Integrate over control volume
Finite Volume Scheme (cont’d) Step 2: Make linear profile assumption between cell cent centro roid ids s for for φ. Assume S varies linearly over CV Step 3: Collect terms and cast into algebraic equation:
Comments
Process starts with conservation statement over cell. We find φ such that it satisfies conservation. Thus, regardless of how coarse the mesh is, the finite volume scheme always gives perfect conservation
This does not guarantee accuracy, however. The proces process s of discretiza discretization tion yields yields a flux balance balance involving face values of the diffusion flux, for example:
∂φ ⎞ ⎛ −Γ ⎜ ⎟ ⎝ ∂ x ⎠ e
e
Profile assumptions for φ and S need not be the same.
Comments (cont’d)
As with finite difference method, we need to solve a set of coupled algebraic equations Though finite difference and finite volume schemes use different procedures to obtain discrete equations, we can use the same solution techniques to solve the discrete equations
Closure In this lecture we
Considered Considered different different canonica canonicall PDEs and examined examined their behavior Understood how these model equations relate to our general scalar transport equations Started an overview of the important elements of any numerical method In the next lecture we will complete this overview and start looking more closely at the finite volume method for diffusion problems.
Lecture 3: Overview of Numerical Methods
Last time…
Examined important classes of partial differential equations and understood their behavior Saw how this knowledge would apply to the general scalar transport equation Started an overview of numerical methods including mesh terminology and finite difference methods
This time… We will continue the overview and examine
Finite difference, finite volume and finite element methods Accuracy, consistency, stability and convergence of a numerical scheme
Overview of Finite Element Method
Consider diffusion equation Let φ be an approximation to φ Since nce φ is an app pprroximation ion, it doe oes s not satisfy the diffusion equation, and leaves a residual R: 2
Γ
d
φ
dx
2
+S = R
Galerkin Galerkin finite finite element method method minimize minimizes s R with respect respect to to a weight function:
Finite Element Method (cont’d)
A family of weight functions Wi, I = 1,…N, (N: number of grid points) is used. This generates N discrete equations for the N unknowns:
Weight function is local – i.e. i.e. zero zero ever everyw ywhe here re except close to i i-1
i Element i-1
w i Element i
i+1
Finite Element Method (cont’d)
In addition a local shape function Ni is used used to d disc iscret retize ize R. Under a Galerki Galerkin n formulation formulation,, the weight weight and shape functio functions ns are chosen to be the same. Shape function is non- zero only in the vicinity of node i => “local basis” N i
N i-1 i-1
i+1
i Element i-1
Element i
Finite Element Method (cont’d)
The discreti discretization zation process process again leads leads to a set set of algebraic algebraic equations of the form:
ai ,iφi
= ai ,i +1φi +1 + ai ,i −1φ i ,i −1 + bi
Comments » Note Note how the use use of a local local basis basis rest restric ricts ts th the e relati relations onship hip between a point i and its neighbors to only nearest neighbors » Agai Again, n, we we have have an alg algeb ebra raic ic equa equati tion on set set to solv solve e – can can u use se the same solvers as for finite volume and finite difference methods
Comparison of methods
All three yield discrete algebraic equation sets which must be solved
Local basis – only near-neighbor near-neighbor dependence dependence
Finite volume method is conservative; the others are not
Order or accuracy of scheme depends on » Taylor Taylor seri series es trun truncat cation ion in in finite finite differ differenc ence e schem schemes es » Profil Profile e assump assumptio tions ns in finite finite volume volume scheme schemes s » Order Order of of s shape hape functi functions ons in finit finite e eleme element nt sche schemes mes
Solution of Linear Equations
Linear equation set has two important characteristics » Ma Matr trix ix is is spa spars rse, e, ma may y be be band banded ed » Coe Coeffi fficie cients nts are provisi provisional onal for non non-lin -linear ear problem problems s Two different approaches » Direct methods » Iterative me metho hod ds Approach defines “path to solution” » Fina Finall answ answer er onl only y dete determ rmin ined ed by by disc discre reti tiza zati tion on
Direct Methods
All disc discret retiza izatio tion n scheme schemes s lead lead to
Here φ is solution vector [φ1 , φ2 ,…, φN]T.
Can invert: Inversion is O(N3 ) operation. Other more efficient methods exist. » Take Take advant advantage age of of band band struct structure ure ifif it exist exists s » Take Take ad adva vant ntag age e of of spa spars rsit ity y
Direct Methods (cont’d)
Large storage and operation count » For For N grid grid po poin ints ts,, mus mustt sto store re NxN NxN ma matr trix ix » Only Only stor store e nonnon-zer zero o entr entries ies and fill fill patt pattern ern
For non-linear problems, A is provisional and is usually updated as a part of an outer loop » No Nott wort worth h solv solving ing syst system em too too “ex “exac actl tly” y”
As a result, direct methods not usually usuall y preferred in CFD today
Iterative Methods
Guess and correct philosophy Gauss-Seidel scheme is typical: » Visi Visitt ea each ch grid grid po poin intt Update using » Sweep Sweep repea repeated tedly ly throu through gh grid grid point points s until until conv converge ergence nce criterion is met » In ea each ch sweep sweep,, points points alre already ady visit visited ed ha have ve new new values values;; points points not yet visited have old values
Iterative Methods (cont’d)
Jacobi scheme is similar similar to Gauss-S Gauss-Seidel eidel scheme scheme but does does not use latest available values » All values values are updated updated simult simultaneo aneousl usly y at end of swee sweep. p.
Iterative are not guaranteed to converge to a solution unless Scarborough criterion is criterion is satisfied
Scarborough Criterion
Scarborough criterion states that convergence of an iterative scheme is guaranteed if:
This means that coefficient matrix must be diagonally dominant
Gauss-Seidel Scheme
No need to store coefficient matrix
Operation count per sweep scales as O(N)
However, convergence, even when guaranteed, is slow for large meshes Will examine alternatives later in course
Accuracy
While looking at finite difference methods, we wrote: Second- order truncation error
Halving grid size reduces error by factor of four for second-order scheme Cannot Cannot say what absolut absolute e error is – truncation truncation error error only gives gives rate of decrease
Accuracy
Order Order of discret discretiza izatio tion n scheme scheme is n n if if truncation error is O(∆xn ) When more than one term is involved, the order of the discretizat discretization ion scheme scheme is that of the the lowest order term. order term. Accuracy Accuracy is a property property of the discretiza discretization tion scheme, scheme, not the path to solution
Consistency
A discretizati discretization on scheme scheme is is consist consistent ent if the truncation truncation error vanishes as ∆x ->0
Does not always happen: What if truncation error is O(∆x/ ∆t) ? Consistenc Consistency y is a property of the discreti discretization zation scheme, scheme, not the path to solution
Convergence
Two uses of the term » Conver Convergen gence ce to to a meshmesh-inde independ pendent ent soluti solution on thro through ugh mesh mesh refinement » Con Conver vergen gence ce of an itera iterativ tive e schem scheme e to a final final unch unchang anging ing answer (or one meeting convergence criterion)
We will usually use the latter meaning
Stability
Property of the path to solution
Typically used to characterize iterative schemes
Depending on the characteristics of the coefficient matrix, errors may either be damped or may grow during iteration An iterative scheme is unstable if it fails to produce a solution to the discrete equation set
Stability
Also possible to speak of the stability of unsteady schemes » Unstable : when solving a time-dependent problem, the solution “blows up” Von-Neumann (and other) stability analyses determine whether linear systems stable under various iteration/time-stepping schemes For non-linear/coupled problems, stability analysis is difficult and not much used » Take Take guidanc guidance e from line linear ar analys analysis is in appro appropri priate ate parameter range; intuition
Closure
This time we completed an overview of the numerical discre discretiz tizati ation on and solutio solution n process process » » » »
Dom omai ain n disc discre reti tiza zati tion on Disc Di scre reti tiza zati tion on of go gove vern rnin ing g eq equa uati tion ons s Solu Soluti tion on of of line linear ar alg algeb ebra raic ic set set Prope Propert rtie ies s of of disc discre reti tiza zati tion on an and d pat path h to to solu soluti tion on – Accuracy, Accuracy, consisten consistency, cy, convergenc convergence, e, stabil stability ity
Next time, time, we will start start looking at finite finite volume discreti discretization zation of diffusion equation
Lecture Lecture 4: The The Diffu Diffusion sion Equati Equation on – A First Look
Last Time…
We completed an overview of the numerical discre discretiz tizati ation on and solut solution ion proce process ss » Doma Domain in dis discr cret etiz izat atio ion n » Di Disc scre reti tiza zati tion on of go gove vern rnin ing g eq equa uati tion ons s – fini finite te difference, finite volume, finite element » Solu Soluti tion on of of line linear ar alg algeb ebra raic ic set set » Prop Proper erti ties es of dis discr cret etiz izat atio ion n an and d pat path h to solu soluti tion on – Accuracy, Accuracy, consisten consistency, cy, convergence convergence,, stability
This Time… We will
Apply the finite volume scheme to the steady diffusion equation on Cartesian structured meshes
Examine the properties of the resulting discretization
Describe Describe how to discret discretize ize boundary boundary conditi conditions ons
2D Steady Diffusion • Con Consid sider er stead steady y diffus diffusion ion with with a source term: • Here
• Int Integr egrate ate over over contro controll volum volume e to yield
2D Steady Diffusion
Apply divergence theorem to yield
Discrete Flux Balance • Writing
integral over control volume:
•Compactly:
Discrete Flux Balance (cont’d)
Area vectors given by:
Fluxes given by
Discretization
Assume φ varies linearly between cell centroids
Note: » Sym Symme mettry of (P, E ) an and d (P,W) in flux expression » Opp ppos osit ite e sig signs ns on (P, P,E E) and (P,W) terms
Source Linearization
Source term must be line linear ariized zed as: as:
Assume SP <0
More on this later!
Final Discrete Equation N
W
P
S
E
Comments
Discrete equation reflects balance of flux*area with generation inside control volume As in 1-D case, we need fluxes at cell faces These are written written in terms terms of cell-ce cell-centroid ntroid values values using using profile assumptions.
Comments (cont’d)
Formulation is conservative: Discrete equation was derived by enforcing conservation. Fluxes balance source term regardless of mesh density For a structured mesh, each point P is coupled to its four nearest neighbors. Corner points do not enter the formulation.
Properties of Discretization
aP, anb have same sign: This implies that if neighbor φ goes up, φP also goes up If S=0:
Thus φ is bounded by neighbor values, in keeping with properties of elliptic partial differential equations
Prope Properti rties es of of Discr Discreti etizat zation ion (cont’ (cont’d) d)
What about Scarborough Criterion ? Satisfied in the equality
What about this?
Boundary Conditions
Flux Balance
Different boundary conditions require different representations of Jb
Dirichlet BCs
Dirich Dirichlet let boundar boundary y conditi condition: on:
φb = φgiven
Put in the requisite flux into the near-boundary cell balance
Diri Dirich chle lett BC’s BC’s (con (cont’ t’d) d) For near-boundary cells: aP
> ∑ anb nb
Satisfies Scarborough Criterion !
Also, φP bounded by interior neighbors and boundary value in the absence of source terms
Neumann BC’s
Neumann boundary conditions : qb given
Replace Jb in cell balance with given flux
Neumann BC’s (cont’d) For Neumann boundaries aP
= ∑ anb nb
So inequality constraint in Scarborough criterion is not satisfied Also, φP is not bounded by interior neighbors and boundary value even in the absence of source term te rms s – th this is is is fine ine because of the added flux at the boundary
Boundary Values and Fluxes
Once we solve for the interior values of φ, we can recover the boundary value of the flux for Dirichlet boundary conditions using
Similarly, for Neumann boundary conditions, we can find the boundary value of φ using
Closure
In this lecture we » De Desc scri ribe bed d th the e disc discre reti tiza zati tion on proc proced edur ure e fo forr th the e diffusion equation on Cartesian meshes » Saw Saw tha thatt the the resu result ltin ing g dis discr cret etiz izat atio ion n proc proces ess s preserves the properties of elliptic equations » Sinc Since e we we get get dia diago gona nall domi domina nanc nce e wit with h Dir Diric ichl hlet et bc, bc, the discret discretizati ization on allows allows us to use iterati iterative ve solvers solvers Next time, we will look at one more boundary condition (Robbins or mixed bc), source linearization and conjugate heat transfer
