CFD 3 Computer Fluid Dynamics 2181106 E 181107

CFD 3 Computer Fluid Dynamics 2181106 E 181107 Solvers, Approximation, Stability, Boundedness of Numerical schemes Remark: foils with „black background“ could be skipped, they are aimed to the more advanced courses Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2013

CFD 3 SOLVERS FEM, BEM, FVM, FD transfer PDE into system of algebraic equations for Tj (nodal pressures, velocities, temperature, concentrations…) solved by ØFinite methods (Gauss, SVD, LU decomposition, frontal methods) N 3 operations are required – suitable for smaller systems. ØIterative methods (GS, multigrid, GMRES, conjugated gradients). Prevails at CFD calculations, characterized by number cells (nodes) of several millions and parallel processing (external as well as internal aerodynamics of cars requires up to 108 finite volumes, solved in clusters of e. g. 512 and more processors) Iterative methods are not so sensitive to round-off errors (that’s why they can be applied for such huge systems) Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

CFD 3 Mathematical Requirements Three Mathematical requirements Øconsistency (discretized equation for must be identical with PDE) order of accuracy (m-with respect time, n-with respect to spatial approximation) Østability (attenuation of round-off errors or glitches of initial conditions) Øconvergency. Lax theorem: consistent and stable numerical scheme converges to exact solution (but it holds only for linear systems) Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

CFD 3 Physical Requirements Three Physical requirements ØConservativeness. Balance of mass should hold exactly at an element level and globally. Fulfilled by FVM (Finite Volume Method). Not exactly satisfied fy FEM. ØBoundedness. Solution should not exhibit local min/max in the absence of internal sources (of mass, momentum or heat). Solution should be bounded by boundary values. Min/max principle. ØTransportivness. Numerical scheme should reflect directionality of information transfer (convection along streamlines) Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

CFD 3 Numerical method-analysis Few examples, how to analyze order of accuracy and stability of suggested numerical schemes (FD methods) Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2013 Benton

CFD 3 Order of accuracy Ti-1 x Ti Ti+1 Taylor expansion Approximation of first derivative = Higher Order Terms Accurate for 2 nd order polynomials = T=1, x, x 2 Accurate for 1 nd order polynomials T=1, x Approximation of second derivative = Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010 Accurate for 3 rd order polynomials T=1, x, x 2, x 3

CFD 3 Order of accuracy Ti-1 x Ti Ti+1 Therefore finite differences substituting derivatives at node xi are ØFirst order ØSecond order ØThird order Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

CFD 3 Stability example (explicit method) Unsteady heat transfer (Fourier equation – parabolic) T-temperature, a-temperature diffusivity Finite difference method EXPLICIT (explicit means that unknown temperatures at a new time level can be expressed explicitly, without necessity to solve a system of algebraic equations). ∆t j-1 What is the order of accuracy? Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010 j ∆x j+1 n

CFD 3 Stability example (explicit method) ∆t j-1 j ∆x j+1 Residual of this PDE is therefore Scheme is consistent, linear with respect time, cubic with respect space. Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010 n+1 n

CFD 3 Stability example (explicit method) Rewrite the explicit formula to the following (explicit) form Unknown temperature at a new time level Known temperatures at “old” time level Rules: ØSum of coefficients must be the same on the left and on the right side (1=A+A+1 -2 A). Why? A constant solution must be fulfilled exactly! ØAll coefficients must be positive for bounded solution. Why? ØScarborough criterion (sum of absolute values of off-diagonal elements diagonal element, criterion for convergency of iterative methods) Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

CFD 3 Stability example (explicit method) So why all coefficients must be positive for bounded solution? Resulting value T is calculated as a weighted average of values (sum of weighting coefficients must be 1). Let us assume only two values for simplicity and T 1<T 2. The solution is bounded if T 1<T<T 2. Let us assume, that the result is not bounded and T<T 1. Then For positive value (1 -A)>0 it follows that T 1>T 2 and this is contradiction. Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

CFD 3 Stability example of unbounded solution Let us consider what would happen if A=1 (negative value 1 -2 A) 19 -7 3 -1 1 Initial condition is 0 in all nodes and only in one node is 1. Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010 Evolution of initial condition in node

CFD 3 Stability example (explicit method) Stability condition can be expressed as a restriction of time step Interpretation in terms of penetration theory. Effective velocity of a thermal pulse Effective domain of dependency Δx-distance of penetrated thermal pulse at a time Δt Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

CFD 3 Stability (hyperbolic/parabolic eqs. ) Stability criterion CFL (Courant Friedrichs Levi) for hyperbolic equations was presented in the previous lecture as where c is the speed of sound or a transport velocity. This CFL criterion represents a linear restriction of the time step with the spatial step and seems to be quite different than the stability criterion for the diffusion driven phenomena. However, these criteria are almost the same, taking into account the penetration depth theory Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010 this is qualitatively the same result (up to a scale constant)

CFD 3 Stability example of wrong scheme Richardson’s scheme for the solution of previous equation n+1 ∆t operates at 3 time levels, n-1, n, n+1 and has higher orders of accuracy j-1 However, the scheme is practically useless. WHY? Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010 j ∆x j+1 n n-1

CFD 3 Stability example of wrong scheme Because this coefficient is always negative n+1 ∆t j-1 j ∆x j+1 n n-1 Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

CFD 3 Stability how to improve Richardson? Richardson’s scheme du. Fort Frankel scheme and this solution will be bounded for A<1/2. Order of accuracy remains high. However: . ∆t j-1 Consistency with Fourier equation is assured only if otherwise the hyperbolic equation of heat transfer would be solved where Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010 n+1 j ∆x j+1 n n-1

CFD 3 Stability hyperbolic equations Leap frog (conditionally stable - von Neumann analysis) Always unstable (negative coefficient at Tj+1) n+1 ∆t. Lax Fridrichs (conditionally stable) Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2013 j-1 j ∆x j+1 n n-1

CFD 3 Stability Neumann More precise (and more complicated) is the stability analysis suggested by von Neumann. It is based upon spectral decomposition of solution, i. e. at a time level n the spatial profile is substituted by Fourier expansion This Fourier component is substituted into differential equation and amplification factor G is evaluated. Numerical scheme is stable, as soon as the magnitude of identified amplification factor decreases. Gain factor G is generally a complex number. Real part is a measure of dumping error (artificial viscosity) and imaginary part is a measure of phase error (dispersion error).

CFD 3 PDE stability analysis (von Neuman) A general Fourier term of solution of a linear partial differential equation km= m/L is wave number (discrete frequency). Arbitrary initial condition T(0, x) can be expressed by Fourier expansion and evolution of individual Fourier terms can be analysed (|exp(at)|>1 indicate growth - instability, |exp(at)|<1 stability). Why is Euler formula unstable? Courant number CFL criterion (Courant-Fridrichs-Levi) Gain G –absolute value greater than 1 for any wave number