Equation Sheet for Quiz 2

MASSACHUSETTS INSTITUTE OF TECHNOLOGY DEPARTMENT OF MECHANICAL ENGINEERING CAMBRIDGE, MASSACHUSETTS 02139 2.29 NUMERICAL FLUID MECHANICS — SPRING 2015 EQUATION SHEET – Quiz 2 Number Representation - Floating Number Representation: x = m be , b-1 ≤ m < b0 Truncation Errors and Error Analysis y  f ( x1, x2 , x3 ,..., xn ) - Taylor Series: f ( xi 1 )  f ( xi )  x f '( xi )  Rn  x 2 x 3 x n n f ''( xi )  f '''( xi )  ...  f ( xi )  Rn 2! 3! n! x n 1 ( n 1) f ( ) n  1! n - The Differential Error (general error propagation) Formula:  y   i 1 - The Standard Error (statistical formula): E (  s y ) - Condition Number of f(x): K p  n  f    x  i 1  i  f ( x1 ,..., xn ) i xi 2 i2 x f '( x ) f (x) Roots of nonlinear equations ( xn1  xn  h( xn ) f ( xn ) ) - Bisection: successive division of bracket in half, next bracket based on sign of n 1 f ( x1n1 ) f ( xmid-point ) f ( xU ) ( xL  xU ) f ( xL )  f ( xU ) - Fixed Point Iteration (General Method or Picard Iteration): - False-Position (Regula Falsi): xr  xU  xn1  g ( xn ) or xn 1  xn  h( xn ) f ( xn ) 1 f ( xn ) f '( xn ) ( xn  xn 1 ) - Secant Method: xn 1  xn  f ( xn ) f ( xn )  f ( xn 1 ) - Newton Raphson: xn 1  xn  - Order of convergence p: Defining en  xn  x e , the order of convergence p exists if there exist a constant C≠0 such that: lim n  en 1 en p C Conservation Law for a scalar , in integral and differential forms: d  d -    dV   dV  CS  (v.n )dA   CS q .n dA  CM  dt  dt CVfixed Advective fluxes (Adv.& diff. ="convection" fluxes) - Other transports (diffusion, etc)  CVfixed s dV Sum of sources and sinks terms (reactions, etc)   .( v )  . q  s t Linear Algebraic Systems: aik( k ) ( k 1) (k ) (k ) ( k 1)  bi( k )  mik bk( k ) , - Gauss Elimination: reduction, mik  ( k ) , aij  aij  mik akj , bi akk  followed by a back-substitution. xk   bk   - LU decomposition: A=LU, aij  min( i , j )  k 1 n a j  k 1 (k ) kj  x j  akk( k )  mik ak( kj ) - Choleski Factorization: A=R*R, where R is upper triangular and R* its conjugate transpose. 1 - Condition number of a linear algebraic system: K ( A)  A A - A banded matrix of p super-diagonals and q sub-diagonals has a bandwidth w = p + q + 1 n n n(n  1) n(n  1)(2n  1) - i  and  i 2  2 6 i 1 i 1 - Eigendecomposition: Ax  x and det(A  I)  0 - Norms: m A 1  max  aij 1 j  n A A n  F  max  aij 1 i  m  “Maximum Column Sum” i 1 “Maximum Row Sum” j 1 2  m n  a   ij   i 1 j 1  A 2  max A*A “Frobenius norm”  or “Euclidean norm”  “ L  2 norm”  or “spectral norm”  Iterative Methods for solving linear algebraic systems: xk 1  B x k  c - Necessary and sufficient condition for convergence:  (B)  max i  1, where i  eigenvalue(Bnn ) k  0,1, 2,... i 1... n - Jacobi’s method: xk 1  -D-1 (L  U) xk  D-1b - Gauss-Seidel method: xk 1  (D  L)-1 U xk  (D  L)-1 b - SOR Method: xk 1  (D  L)1[U  (1  )D]x k  (D  L)1 b ri T ri ri , ri  b  Axi - Steepest Descent Gradient Method: xi 1  xi  T ri Ari - Conjugate Gradient: xi 1  xi  i vi (αi such that each vi are generated by orthogonalization of residuum vectors and such that search directions are A-conjugate). Finite Differences – PDE types (2nd order, 2D): Axx  Bxy  Cyy  F ( x, y,  , x , y ) B2  AC > 0: hyperbolic; B2  AC = 0: parabolic; B2  AC < 0: elliptic ( )  0 , Finite Differences – Error Types and Discretization Properties ( - Consistency: ( )  ˆx ( )  0 when x  0 - Truncation error:  x  - Error equation:  x  - Stability: ˆ (ˆ)  0 ) x ( )  ˆx ( )  O(x p ) for x  0 ( )  ˆx (ˆ   )   ˆx ( ) (for linear systems) ˆ 1  Const. x - Convergence:   ˆ 1 x (for linear systems)  x   O(x p ) Finite Differences – General schemes and Higher Accuracy s   mu  Higher Order Accuracy Finite-difference based on Taylor Series:  m    ai u j i   x  x  j i  r Newton’s interpolating polynomial formulas, equidistant sampling: Lagrange polynomial: n f ( x )   Lk ( x ) f ( xk ) with Lk ( x )  k 0 Hermite Polynomials and Compact/Pade’s Difference schemes: n x  xj j 0, j k xk  x j  s q   mu   ai u j i   x  m    j i i  p  b  x i  r i Finite Differences – Non-Uniform Grids, Grid Refinement and Error Estimation For a centered-difference approximation of f’(x) over a 1D grid, contracting/expanding with a constant factor re, xi 1  re xi , the: (1  re ) xi f ''( xi ) 2 (1  re,h )2 - Ratio of the two truncation errors at a common point is: R  re,h Grid-Refinement and Error estimation:  u  u4 x  - Estimate of the order of accuracy: p  log  2 x  log 2  ux  u2 x  - Leading term of the truncation error is:  rex  - Discretization error on the grid Δx:  x  u x  u2 x 2p 1 Richardson Extrapolation for the Trapezoidal Rule: “Romberg” Differentiation Algorithm: D j ,k  2 4k 1 D j 1,k 1  D j ,k 1 4k 1  1 Finite Differences – Fourier Analysis and Error Analysis Fourier transform of a generic PDE, f t  n f : With f (x,t)  f k (t ) ei kx , one obtains:  n x k  d f k (t ) n   ik  f k (t )   f k (t ) dt for    ik  n Finite Difference Methods – Effective wave number and speed   ei kx  sin(k x ) i kx j (for CDS, 2nd order, keff  ) Effective Wave Number:    i keff e x  x  j Effective Wave Speed (for linear convection eq., f f c  0 ): t x   d f knum.   f knum. (t ) c i keff  f numerical ( x, t )   f k (0) ei kx i keff t   f k (0) ei k (x ceff t ) dt k  k   ceff c   eff keff   k Finite Difference Methods – Stability Von Neumann:  ( x, t )     ( t ) e  ,   (t ) e    i x i x  e t ei  x (γ in general complex, function of β)  Strict condition for stability: e t  1 or for   e t ,   1  (for the error not to grow in time) Useful trigonometric relations: ei x  ei x  2cos(x), ei x  ei x  2sin(x) and 1  cos(x)  2sin 2 (x / 2) CFL condition: “Numerical domain of dependence of FD scheme must include the mathematical domain of dependence of the corresponding PDE” Figure 23.1 Chapra and Canale Forward Differences Backward Differences Centered Differences Finite Difference Methods – Schemes for specific PDE types Hyperbolic, 1D: ux + b uy = 0 Elliptic PDEs: 2D Laplace/Poisson Eq. on a Cartesian-orthogonal uniform grid u k  uik1, j  uik, j 1  uik, j 1  h 2 gi , j SOR, Jacobi: uik, j 1  (1   ) uik, j   i 1, j 4 k k 1 ui 1, j  ui 1, j  uik, j 1  uik, j 11  h 2 gi , j k 1 k SOR, Gauss-Seidel: ui , j  (1   ) ui , j   4 Parabolic PDEs: 2D Heat Conduction Eq. on a Cartesian-orthogonal grid n n n n n n Ti ,nj1  Ti ,nj 2 Ti 1, j  2Ti , j  Ti 1, j 2 Ti , j 1  2Ti , j  Ti , j 1 c c • Explicit: t x 2 y 2 • Crank-Nicolson Implicit (for Δx=Δy, with r  (1  2r )Ti ,nj1  (1  2r )Ti ,nj  Ti ,nj1/ 2  Ti ,nj t c 2 ): x 2 r n 1 r Ti 1, j  Ti n1,1j  Ti ,nj11  Ti ,nj11   Ti n1, j  Ti n1, j  Ti ,nj 1  Ti ,nj 1   2 2 2 Ti n1, j  2Ti ,nj  Ti n1, j 2 Ti ,nj1/1 2  2Ti ,nj1/ 2  Ti ,nj1/1 2 c c t / 2 x 2 y 2 n 1 n 1 n 1 n 1/ 2 n 1/ 2 Ti ,nj1  Ti ,nj1/ 2  Ti ,nj1/1 2 2 Ti 1, j  2Ti , j  Ti 1, j 2 Ti , j 1  2Ti , j c c t / 2 x 2 y 2 • ADI: n 1/2 n 1/2 n 1/2 n n n (for Δx=Δy): rTi , j 1  2(1  r )Ti , j  rTi , j 1  rTi 1, j  2(1  r )Ti , j  rTi 1, j n 1/2 rTi n1,1j  2(1  r)Ti ,nj1  rTi n1,1j  rTi ,nj1/2  rTi ,nj1/2 1  2(1  r )Ti , j 1 Finite Volume Methods: V d 1   F  .n dA  S , where     dV and S   s dV S V dt V V Cartesian grids  Surface Integrals: Fe   f dA Se - 2D problems (1D surface integrals)  Midpoint rule (2nd order): Fe   f dA  f e Se  f e Se  O( y 2 )  f e Se  Trapezoid rule (2nd order): Fe   f dA  Se  Simpson’s rule (4th order): Fe   Se Se Se - 3D problems (2D surface integrals)  ( f ne  f se )  O( y 2 ) 2 ( f ne  4 f e  f se ) f dA  Se  O( y 4 ) 6 Midpoint rule (2nd order): Fe   f dA  Se f e Se  Volume Integrals: S   s dV ,   V  O( y 2 , z 2 ) 1 dV V V - 2D/3D problems, Midpoint rule (2nd order): SP   s dV  sP V  sP V V - 2D, bi-quadratic (4th order, Cart.): SP  x y 16sP  4ss  4sn  4sw  4se  sse  ssw  sne  snw  36  Interpolations / Differentiations (obtain fluxes “Fe= f (e)” as a function of cell-average values) P if  v . n e  0 - Upwind Interpolation (UDS): e   E if  v . n e  0 x  xP - Linear Interpolation (CDS): e  E e  P (1  e ) where e  e xE  xP   E   P (1   ), with      P x  xP  E  x e xE  xP xE  xP - Quadratic Upwind interpolation (QUICK): e  U  g1 (D  U )  g2 (U  UU ) 6 3 1 3x 3  3 For uniform grids, e  U  D  UU  8 8 8 48 x 3 - Higher order schemes: For example, for  ( x)  a0  a1 x  a2 x 2  a3 x 3 , 27P  27E  3W  3EE Convective fluxes e 48 Diffusive Fluxes, for a uniform Cartesian grid:  x For a compact high order scheme: e  P  E Solution of the Navier-Stokes Equations Newtonian fluid + incompressible + constant: 2  e  R3 D 27E  27P  W  EE . 24 x x       O( x 4 )  8  x P x E   v  .(  v v )  p  2v   g t .v  0 Strong conservative form, general Newtonian fluid:    u u   vi 2 u  .(  vi v )  .   p ei    i  j  e j   j ei   gi xi ei   x    t 3 x j  j xi    Kinetic energy equation, CV form: 2 2 v v   dV     (v .n )dA   p v .n dA   ( .v ). n dA    : v  p .v   g.v dV CS CS CS CV t CV 2 2   v  .  .(  v v )   .  2v   .   g  t For constant μ and ρ: .p  .  .(  v v )  Pressure equation: .p  2 p  .  Pressure-correction Methods Hi    (  ui u j )   ij +  xj  xj Forward-Euler Explicit in Time:   n  pn n 1 n  u   u   t       Hi   i i  xi     n n     p   Hi   x   x   x i  i i  Backward-Euler Implicit in Time:    (  ui u j ) n 1   ij n 1  p n 1  n 1 n +    ui     ui   t      x  x  x  j j i    n 1 n 1 n 1   ij      p     (  ui u j )         x  x   xj  x j   i  i   xi  Backward-Euler Implicit in Time, linearized momentum update:   (  ui u j ) n  (  uin u j )  (  ui u nj )   ij n   ij  p n p  n 1 n         ui     ui    ui  t     x  x  x  x  x  x  xi  j j j j j i  Steady state solver, matrix notation: uim* Outer iteration, nonlinear solve: A u Outer iteration, pressure update: δ uim* δxi Inner iteration, linear solve: uim* m* i b m 1 uim* δp  δxi δ  uim*  A δxi    m i A u b m uim* δp  δxi m 1 1 m δp  , δxi    m* uim*  A ui 1 bmum*1 m Steady state solver, matrix notation, pressure-correction schemes: Based on the above, but introduce uim  uim*  u' p m  p m1  p ' and further simplify to get varied schemes (SIMPLE, SIMPLER, SIMPLEC, PISO, etc.) i Projection Methods, Pressure-Correction Form Non-Incremental: u  i * n 1    ui  n   (  ui u j ) n 1   ij n 1   t   +  ;   x  x j j      p n 1  1     xi   xi  t  xi   ui  n 1    ui *  n 1   u   i  t * n 1 ; u  i * n 1 D  (bc)  p n 1 0  n D  p n 1  xi Incremental: u  i * n 1    ui  n   (  ui u j ) n 1   ij n 1  p n   t   +   ;   x  x  x j j i   n 1 n   p  p  1     t  xi  xi   xi    ui  n 1    ui *  n 1  t   u   ; i  p n 1 * n 1 p n u  i   p n 1  p n  n * n 1 D  (bc) 0 D   xi Rotational Incremental: u  i * n 1    ui      p   xi   xi   ui  n 1 n 1 n   (  ui u j ) n 1   ij n 1  p n   t   +  ;   xj  xj  xi           ui *  n 1 1  t  xi  t p n 1  p n   p n 1     u   ; i   p   xi * n 1 n 1   xi  u   i * n 1   p n 1  n 0 D u  i * n 1 D  (bc)