Aman Raj

COMPUTATIONAL FLUID DYNAMICS May 11, 2025 QUESTION A: Start with 2D incompressible Navier-Stokes equations: ∂2u ∂2u + 2 ∂x2 ∂y 2 ∂p 1 ∂2v ∂v ∂(uv) ∂(vv) ∂ v + + =− + + ∂t ∂x ∂y ∂y Re ∂x2 ∂y 2 ∂u ∂(uu) ∂(vu) ∂p 1 + + =− + ∂t ∂x ∂y ∂x Re ∂u ∂v + =0 ∂x ∂y Figure 1: Description of the image. 1 SOLUTION: In this question we have to solve the 2D incompressible Navier-Stokes equations which is given in the question using the condition given. We know that N-S equation for x-momentumn is: ∂u ∂(uu) ∂(vu) ∂p 1 + + =− + ∂t ∂x ∂y ∂x Re ∂2u ∂2u + 2 ∂x2 ∂y ∂v ∂(uv) ∂(vv) ∂p 1 + + =− + ∂t ∂x ∂y ∂y Re ∂2v ∂2v + ∂x2 ∂y 2 y-momentum is: Continuity is: ∂u ∂v + =0 ∂x ∂y Now, from the given figure we come to know that u(x, 1) = 1, v(x, 1) = 0 u(0, y) = 0, v(0, y) = 0 u(1, y) = 0, v(1, y) = 0 u(x, 0) = 0, v(x, 0) = 0 for Top wall, Left wall, Right wall, Bottom wall respectively. Also for For a grid with ghost cells at i = 0 and i = Nx + 1 (x-direction), j = 0 and j = Ny + 1 (y-direction): - Left Wall (x = 0): u0,j = −u1,j , v0,j = v1,j , p0,j = p1,j - Right Wall (x = 1): uNx +1,j = −uNx ,j , vNx +1,j = vNx ,j , pNx +1,j = pNx ,j - Bottom Wall (y = 0): ui,0 = −ui,1 , vi,0 = −vi,1 , pi,0 = pi,1 - Top Wall (y = 1): ui,Ny +1 = 2 − ui,Ny , vi,Ny +1 = vi,Ny , pi,Ny +1 = pi,Ny From the past knowledge we know that Convection Terms for Explicit Forward Euler is given by ∂(uu) ∂x ∂(vu) ∂y n ≈ x=0 n ≈ x=0 (uni+1,j )2 − (uni−1,j )2 2∆x n n n vi,j ui,j − vi,j−1 uni,j−1 ∆y and Viscous Terms for Implicit Backward Euler is given by 2 ∂2u ∂x2 ∂2u ∂y 2 n+1 ≈ n+1 n+1 un+1 i+1,j − 2ui,j + ui−1,j ∆x2 ≈ n+1 n+1 un+1 i,j+1 − 2ui,j + ui,j−1 ∆y 2 x=0 n+1 x=0 Now, using the above data we will get the intermediate velocity u∗ and v ∗ We have x-Momentum: ! n n n uni,j−1 ui,j − vi,j−1 vi,j (uni+1,j )2 − (uni−1,j )2 + 2∆x ∆y ∗ ∗ ∗ ∗ ui+1,j − 2ui,j + ui−1,j ui,j+1 − 2u∗i,j + u∗i,j−1 1 + + Re ∆x2 ∆y 2 u∗i,j − uni,j =− ∆t y-Momentum: ! n n n n uni,j vi,j − uni−1,j vi−1,j (vi,j+1 )2 − (vi,j−1 )2 + ∆x 2∆y ∗ ∗ ∗ ∗ ∗ ∗ vi+1,j − 2vi,j + vi−1,j vi,j+1 − 2vi,j + vi,j−1 1 + + Re ∆x2 ∆y 2 ∗ n vi,j − vi,j =− ∆t After Rearranging we get For u∗ : 2∆t 2∆t + 1+ 2 Re∆x Re∆y 2 ∆t u∗i+1,j + u∗i−1,j 2 Re∆x ∆t − u∗ + u∗i,j−1 Re∆y 2 i,j+1 " # n n n (uni+1,j )2 − (uni−1,j )2 vi,j ui,j − vi,j−1 uni,j−1 n = ui,j + ∆t − − 2∆x ∆y u∗i,j − For v ∗ : 1+ 2∆t 2∆t + 2 Re∆x Re∆y 2 ∆t ∗ ∗ vi+1,j + vi−1,j 2 Re∆x ∆t ∗ ∗ − vi,j+1 + vi,j−1 2 Re∆y " # n n n n uni,j vi,j − uni−1,j vi−1,j (vi,j+1 )2 − (vi,j−1 )2 n = vi,j + ∆t − − ∆x 2∆y ∗ vi,j − Applying boundary condition for u∗ and v ∗ - Top Lid (j = Ny ): u∗i,Ny = 1, ∗ vi,N =0 y - Other Walls: u∗0,j = −u∗1,j , u∗Nx +1,j = −u∗Nx ,j 3 (No-slip) ∗ ∗ v0,j = v1,j , u∗i,0 = ∗ ∗ vN = vN x +1,j x ,j −u∗i,1 , ∗ vi,0 (Zero normal gradient) ∗ −vi,1 = (Bottom wall) After solving we get the intermediate velocities as 1 2 ∗ u∗i,j = uni,j + ∆t Convu + ∇ u , Re 1 2 ∗ ∗ n vi,j = vi,j + ∆t Convv + ∇ v . Re where Convu = − Convv = − uni+1,j n uni,j vi,j 2 2 ! n n n vi,j ui,j − vi,j−1 uni,j−1 + ∆y ! n n n (vi,j+1 )2 − (vi,j−1 )2 − uni−1,j vi−1,j + ∆x 2∆y − uni−1,j 2∆x Now, we will use u∗ and v ∗ to solve the Pressure Poisson equation The Pressure Poisson Equation is derived from the incompressibility condition: 1 ∇ · u∗ ∆t where u∗ is the intermediate velocity field from the predictor step. ∇2 p = (1) Discretization Using second-order central differences on a uniform grid (∆x = ∆y = h): ∗ ∗ vi,j+1 − vi,j−1 pi,j+1 − 2pi,j + pi,j−1 1 u∗i+1,j − u∗i−1,j pi+1,j − 2pi,j + pi−1,j + = + h2 h2 ∆t 2h 2h (2) Simplified Form Rearranging terms gives the final discretized form: pi+1,j + pi−1,j + pi,j+1 + pi,j−1 − 4pi,j = h ∗ ∗ u∗i+1,j − u∗i−1,j + vi,j+1 − vi,j−1 2∆t Boundary Conditions • Neumann BCs (∂p/∂n = 0 on walls): p0,j = p1,j pNx +1,j = pNx ,j pi,0 = pi,1 • Reference pressure: pNx /2,Ny = 0 at top center 4 (Left wall) (Right wall) (Bottom wall) (3) Final Expression The pressure at each grid point is obtained by solving: 1 h ∗ ∗ u∗i+1,j − u∗i−1,j + vi,j+1 − vi,j−1 pi,j = pi+1,j + pi−1,j + pi,j+1 + pi,j−1 − 4 2∆t (4) This is solved iteratively (e.g., Gauss-Seidel) until convergence. After the we can easily update velocity with pressure gradient Corrector Step The final velocity field un+1 is obtained by correcting the intermediate velocity u∗ with the pressure gradient: un+1 = u∗ − ∆t∇p (5) Discretized Formulation Using second-order central differences on a collocated grid (∆x = ∆y = h): un+1 i,j = u∗i,j n+1 ∗ vi,j = vi,j pi+1,j − pi−1,j − ∆t 2h pi,j+1 − pi,j−1 − ∆t 2h (6) (7) Boundary Conditions • No-slip walls: un+1 i,1 = 0 (Bottom wall) un+1 1,j n+1 uNx ,j n+1 vi,1 n+1 v1,j n+1 vN x ,j =0 (Left wall) =0 (Right wall) =0 (Bottom wall) =0 (Left wall) =0 (Right wall) • Moving lid (Top wall): un+1 i,Ny = 1 n+1 vi,N =0 y • Ghost cells (for pressure gradient at boundaries): p0,j = p1,j pNx +1,j = pNx ,j pi,0 = pi,1 pi,Ny +1 = pi,Ny (Left) (Right) (Bottom) (Top) The updated velocity field satisfies: ∇ · un+1 ≈ n+1 n+1 n+1 un+1 vi,j+1 − vi,j−1 i+1,j − ui−1,j + =0 2h 2h 5 (8) QUESTION B: Derive the finite difference expressions for: • Forward Euler for time derivatives • Explicit first order treatment of the non-linear convection terms • Implicit treatment of the viscous terms using central difference • Central difference for pressure Laplacian • RHS of the pressure Poisson equation Derive these expressions for (i) an interior cell, (ii) a cell on the four boundary walls, and (iii) four corner cells. The expressions should be in terms of values of the variables at the neighboring cells and/or the boundary values. Note that the boundary conditions are defined on the surface of the cell; you can use a linear interpolation scheme to find the node values of the boundary cells. SOLUTION: For Interior Cell: • Forward Euler for time derivative We know that general form of the unsteady term in the Navier-Stokes equation for velocity u(x, y, t): ∂u ∂t The Forward Euler method approximates the time derivative as: ∂u ∂t n ≈ i,j n un+1 i,j − ui,j ∆t This expression estimates the rate of change of velocity at the point (i, j) between time levels tn and tn+1 using a forward difference. – uni,j : known velocity at the current time step – un+1 i,j : unknown velocity at the next time step – ∆t: time step Thus, we isolate un+1 i,j to update the value explicitly: n un+1 i,j = ui,j + ∆t · (RHS of Navier-Stokes) Where ”RHS” includes convection, diffusion, and pressure gradient terms evaluated at time n (for explicit method). Just the time derivative term written explicitly for any scalar or velocity component ϕ at an interior cell (i, j): n ϕn+1 ∂ϕ i,j − ϕi,j n n ≈ ⇒ ϕn+1 i,j = ϕi,j + ∆t · RHSi,j ∂t i,j ∆t This is the explicit time advancement for Forward Euler. 6 • Explicit 1st order convection terms We know that 2D Convection Term in the x-Momentum Equation We want to discretize the convective term in the x-momentum equation: ∂u u ∂x |{z} ∂u v ∂y |{z} + x-convection y-convection Assume positive flow direction for both u and v, i.e., flow moves left to right and bottom to top. – First-order Upwind in x-direction: ui,j − ui−1,j ∂u ≈ ∂x ∆x if ui,j > 0 – First-order Upwind in y-direction: ∂u ui,j − ui,j−1 ≈ ∂y ∆y ∂u ∂u +v u ∂x ∂y ≈ ui,j · i,j if vi,j > 0 ui,j − ui−1,j ui,j − ui,j−1 + vi,j · ∆x ∆y This is a first-order upwind discretization and is explicit because all terms are evaluated at time level n. We get the Full Expression for x-Momentum Convection Term (Upwind) Convectionu,i,j = uni,j · uni,j − uni−1,j uni,j − uni,j−1 n + vi,j · ∆x ∆y This would be plugged into the RHS of the time advancement equation from the Forward Euler step. For any component like ϕ, the upwind formula is:  ϕ − ϕi−1,j   i,j  ∂ϕ ∆x ≈ ∂x  ϕ   i+1,j − ϕi,j ∆x if ui,j > 0 if ui,j < 0 For v momentum convection term We know that 2D Convection Term in the y-Momentum Equation For the v-component of velocity, the convective term is: ∂v ∂v u + v ∂x ∂y |{z} |{z} x-convection y-convection Again, assume positive flow: u > 0, v > 0 – x-direction (∂v/∂x): ∂v vi,j − vi−1,j ≈ ∂x ∆x if ui,j > 0 ∂v vi,j − vi,j−1 ≈ ∂y ∆y if vi,j > 0 – y-direction (∂v/∂y): 7 So, the convection term becomes: ∂v vi,j − vi−1,j vi,j − vi,j−1 ∂v +v ≈ ui,j · + vi,j · u ∂x ∂y i,j ∆x ∆y Hence, We get the Full Discretized Expression as: Convectionv,i,j = uni,j · n n n n vi,j − vi−1,j vi,j − vi,j−1 n + vi,j · ∆x ∆y This is evaluated explicitly at time level n, matching the Forward Euler time scheme. • Implicit viscous terms for y-Momentum Equation For the v-component, the viscous diffusion term is: 2 ∂ v ∂2v ν + 2 ∂x2 ∂y We will now discretize both second derivatives using second-order central differences, and treat them implicitly (at time level n + 1). 1. ∂2v : ∂x2 ∂2v ∂x2 ≈ n+1 n+1 n+1 vi+1,j − 2vi,j + vi−1,j (∆x)2 ≈ n+1 n+1 n+1 vi,j+1 − 2vi,j + vi,j−1 (∆y)2 i,j 2 2. ∂ v : ∂y 2 ∂2v ∂y 2 i,j After Combining Implicit Viscous Term (∇2 v)n+1 i,j ≈ n+1 n+1 n+1 n+1 n+1 n+1 vi,j+1 − 2vi,j + vi,j−1 vi+1,j − 2vi,j + vi−1,j + (∆x)2 (∆y)2 Hence,we get the Full Expression for Viscous Term as " n+1 # n+1 n+1 n+1 n+1 n+1 vi+1,j − 2vi,j + vi−1,j vi,j+1 − 2vi,j + vi,j−1 n+1 Viscousv,i,j = ν + (∆x)2 (∆y)2 • Central difference for pressure Laplacian In the pressure Poisson equation (PPE), the Laplacian of pressure is: ∇2 p = ∂2p ∂2p + 2 2 ∂x ∂y We know that Central Difference Approximation is as follows 1. ∂2p discretization: ∂x2 ∂2p ∂x2 ≈ i,j pi+1,j − 2pi,j + pi−1,j (∆x)2 8 2. ∂2p discretization: ∂y 2 ∂2p ∂y 2 pi,j+1 − 2pi,j + pi,j−1 (∆y)2 ≈ i,j After Combining Pressure Laplacian (for PPE) we get: Combining both gives the second-order central difference approximation: ∇2 p i,j ≈ pi+1,j − 2pi,j + pi−1,j pi,j+1 − 2pi,j + pi,j−1 + 2 (∆x) (∆y)2 This expression is used on the left-hand side of the pressure Poisson equation. • RHS of pressure Poisson equation We know that General Form for incompressible flow, the PPE is derived from continuity and momentum equations: ρ ∂u∗ ∂v ∗ 2 ∇ p= + ∆t ∂x ∂y where: – u∗ , v ∗ are intermediate velocities (without pressure correction) – ρ is fluid density – ∆t is time step – ∇ · ⃗u∗ is the divergence of intermediate velocity 1. ∂u∗ discretization: ∂x ∂u∗ ∂x ≈ u∗i+1,j − u∗i−1,j 2∆x ≈ ∗ ∗ vi,j+1 − vi,j−1 2∆y i,j ∗ 2. ∂v discretization: ∂y ∂v ∗ ∂y i,j We get the RHS Expression as: ρ ∆t ∗ ∗ u∗i+1,j − u∗i−1,j vi,j+1 − vi,j−1 + 2∆x 2∆y RHSi,j = ρ ∆t ∗ ∗ u∗i+1,j − u∗i−1,j vi,j+1 − vi,j−1 + 2∆x 2∆y RHSi,j = Hence, at the end we get This term ensures the updated velocity field satisfies the continuity equation (∇ · ⃗u = 0). 9 For Boundary wall cells: • Forward Euler for time derivative We know that the General Form of Forward Euler and explicit time update for any velocity component (say u) using Forward Euler is: n n un+1 ∂u i,j − ui,j ≈ ∂t ∆t So, n n un+1 i,j = ui,j + ∆t · RHSu Similarly, n+1 n vi,j = vi,j + ∆t · RHSnv where RHSnu and RHSnv are the convection and diffusion terms evaluated at time level n. For boundary cells, velocity values are not solved directly from Navier-Stokes equations but are instead overwritten by boundary conditions at each time step. At the top boundary: un+1 i,N = 1 n+1 and vi,N =0 The boundary condition is enforced directly regardless of internal updates. – Bottom wall: un+1 i,0 = 0, n+1 vi,0 =0 un+1 0,j = 0, n+1 v0,j =0 un+1 N,j = 0, n+1 vN,j =0 – Left wall: – Right wall: Hence, at the end we get While the Forward Euler time derivative mathematically is: n n un+1 i,j = ui,j + ∆t · RHSu for boundary cells this is not used directly. Instead: n+1 – Top lid: un+1 i,N = 1, vi,N = 0 – Other walls: u = 0, v = 0 • Explicit 1st order convection terms We know that General Form of Convection Terms is as follows The convective terms in the momentum equations are: For u-momentum: For v-momentum: ∂(uv) ∂u2 + ∂x ∂y ∂(uv) ∂v 2 + ∂x ∂y In the upwind scheme: – If u > 0: use backward difference (upstream value) – If u < 0: use forward difference (downstream value) Now, Application to Left Wall Boundary (i = 0) 10 0.1 ∂u2 at i = 0 ∂x Given u0,j = 0 (wall) and assuming u1,j > 0 (flow into domain): ∂u2 ∂x 0.2 ≈ 0,j u21,j − u20,j u21,j = ∆x ∆x (since u0,j = 0) ∂(uv) at i = 0 ∂y Assuming v0,j > 0 (upward flow): ∂(uv) ∂y ≈ 0,j u0,j v0,j − u0,j−1 v0,j−1 u0,j−1 v0,j−1 =− ∆y ∆y (since u0,j = 0 due to no-slip wall) After Combining u-momentum convection term we get 2 u21,j ∂u ∂(uv) u0,j−1 v0,j−1 + = − ∂x ∂y ∆x ∆y i=0,j – For u < 0, use forward difference instead – At top lid boundary (u = 1, v = 0), upwind based on these fixed values – Boundaries with zero velocity lead to significant simplifications • Implicit viscous terms for y-Momentum Equation We know that Viscous Terms in Navier-Stokes (xmomentum) is as follows 2 ∂ u ∂2u ν + ∂x2 ∂y 2 Left Wall Boundary (i = 0) Discretization 0.3 ∂ 2u at i = 0 ∂x2 Using central difference with ghost node: ∂2u ∂x2 ≈ 0,j u1,j − 2u0,j + u−1,j (∆x)2 With no-slip condition (u0,j = 0) and ghost node reflection (u−1,j = −u1,j ): u1,j − 0 + (−u1,j ) =0 (∆x)2 11 0.4 ∂ 2u at i = 0 ∂y 2 Using interior points: ∂2u ∂y 2 u0,j+1 − 2u0,j + u0,j−1 u0,j+1 + u0,j−1 = (∆y)2 (∆y)2 ≈ 0,j Hence, we get the Final Viscous Term (x-momentum) at left wall 2 ∂ u ∂2u u0,j+1 + u0,j−1 ν + = ν 0 + ∂x2 ∂y 2 (∆y)2 Similarly, Bottom Wall Boundary (j = 0) for y-momentum 2 ∂ v ∂2v ν + ∂x2 ∂y 2 Now, Using ghost node reflection (vi,−1 = −vi,1 ): ∂2v ∂x2 ≈ i,0 2 ∂ v ∂y 2 ≈ i,0 vi+1,0 − 2vi,0 + vi−1,0 (∆x)2 vi,1 − 2vi,0 + vi,−1 −2vi,0 = 2 (∆y) (∆y)2 We, get the Final Viscous Term (y-momentum) at bottom wall as the below: 2 ∂ v ∂2v 2vi,0 vi+1,0 − 2vi,0 + vi−1,0 ν + − = ν ∂x2 ∂y 2 (∆x)2 (∆y)2 • Central difference for pressure Laplacian We know that Poisson Equation for Pressure is: ∇2 p = ∂2p ∂2p + ∂x2 ∂y 2 At Left Boundary (i = 0) Discretization 0.5 ∂ 2p at i = 0 ∂x2 Using central difference with ghost node: ∂2p ∂x2 With Neumann boundary condition ≈ 0,j p1,j − 2p0,j + p−1,j (∆x)2 ∂p = 0 at the wall: ∂x p1,j − p−1,j = 0 =⇒ p−1,j = p1,j 2∆x Thus: ∂2p 2(p1,j − p0,j ) ≈ ∂x2 (∆x)2 12 0.6 ∂ 2p at i = 0 ∂y 2 Using standard central difference: ∂2p ∂y 2 ≈ 0,j p0,j+1 − 2p0,j + p0,j−1 (∆y)2 Hence, We get the Final Pressure Laplacian at i = 0 as ∇2 p 0,j = 2(p1,j − p0,j ) p0,j+1 − 2p0,j + p0,j−1 + (∆x)2 (∆y)2 This represents the central difference approximation of the pressure Laplacian at a left boundary cell using: – Ghost nodes with Neumann boundary condition – Standard central difference for interior y-direction – Uniform grid spacing (∆x = ∆y) • RHS of pressure Poisson equation We know that General Form for the pressure Poisson equation in incompressible flow: ρ ∂u∗ ∂v ∗ 2 n+1 ∇ p = + ∆t ∂x ∂y where u∗ , v ∗ are intermediate velocities and ρ is fluid density (often ρ = 1). Left Wall Boundary (i = 0, j) Using Divergence Approximation we get – – ∂u∗ at (0, j): ∂x ∂u∗ ∂x ≈ 0,j ∂v ∗ at (0, j): ∂y ∂v ∗ ∂y ≈ 0,j u∗1,j − u∗0,j ∆x ∗ ∗ v0,j+1 − v0,j−1 2∆y (forward difference) (central difference) At RHS Expression Assuming ρ = 1: RHS0,j 1 = ∆t ∗ ∗ u∗1,j − u∗0,j v0,j+1 − v0,j−1 + ∆x 2∆y Hence, we get the Final Form for Left Boundary as: ∗ ∗ v0,j+1 − v0,j−1 1 u∗1,j − u∗0,j RHS0,j = + ∆t ∆x 2∆y 13 For Corner cells: • Forward Euler for time derivative For the bottom-left corner cell (i, j): – Velocity: U = 0 (left wall), V = 0 (bottom wall) ∂p ∂p – Pressure: = 0 (left wall), = 0 (bottom wall) ∂x ∂y For any variable ϕ (U, V, or p): n n ϕn+1 i,j = ϕi,j + ∆t · (RHS ) where: – ϕni,j = current time step value – ϕn+1 i,j = next time step value – ∆t = time step – RHSn = right-hand side terms evaluated at time n We know that U-Momentum Equation for Corner Cell is n 2 n n ∂U ∂U ∂ U ∂2U ∂p n+1 n Ui,j = Ui,j + ∆t − U +V +ν + − 2 2 ∂x ∂y ∂x ∂y ∂x We have Convection Terms (Explicit, 1st-order Upwind) as n – If Ui,j ≥ 0: n Ui,j −0 ∂U ≈ ∂x ∆x n – If Vi,j ≥ 0: n Ui,j −0 ∂U ≈ ∂y ∆y We have Viscous Terms (Central Difference) as n n Ui+1,j − 2Ui,j ∂2U ≈ , ∂x2 ∆x2 n n Ui,j+1 − 2Ui,j ∂2U ≈ ∂y 2 ∆y 2 We have Pressure Gradient (Neumann BC) as ∂p ≈0 ∂x (using pni−1,j = pni+1,j ) Hence, we get the Final Explicit Update for U as follows n n n n n n Ui,j Vi,j Ui+1,j − 2Ui,j Ui,j+1 − 2Ui,j n+1 n n Ui,j = Ui,j + ∆t −Ui,j + +ν + − 0 ∆x ∆y ∆x2 ∆y 2 • Explicit 1st order convection terms For the bottom-left corner cell (i, j): – U = 0 (left wall) – V = 0 (bottom wall) 14 We know that U-Momentum Convection Terms is U ∂U ∂U +V ∂x ∂y Writing Convection in x-Direction: U ∂U ∂x n n 2 Ui,j (Ui,j −0 ) ∂U ⇒ U ∂U ∂x ≈ ∆x ∂x ≈ ∆x n n Ui+1,j −Ui,j ∂U n ⇒ U ∂U ∂x ≈ ∆x ∂x ≈ Ui,j Uni,j ≥ 0 (flow to right) : n Ui,j < 0 (flow to left) : Writing Convection in y-Direction: V Vni,j ≥ 0 (flow upward) : ∂U ∂y n < 0 (flow downward) : Vi,j n n Ui+1,j −Ui,j ∆x ∂U ∂y n n n Ui,j Vi,j −0 Ui,j ⇒ V ∂U ∆y ∂y ≈ ∆y n n Ui,j+1 −Ui,j ∂U n ⇒ V ∂U ∂y ≈ ∆y ∂y ≈ Vi,j ≈ n n Ui,j+1 −Ui,j ∆y Hence, the Combined Convection final Term is given The total convection term for U is: n Uni,j ≥ 0, Vi,j ≥0: n n Ui,j < 0, Vi,j <0: n 2 n n (Ui,j ) Vi,j Ui,j ∆x + ∆y n n Ui+1,j −Ui,j n Ui,j ∆x n + Vi,j n n Ui,j+1 −Ui,j ∆y n n Mixed cases (e.g., Ui,j ≥ 0, Vi,j < 0) are handled similarly by combining the appropriate cases. • Implicit viscous terms for y-Momentum Equation For the bottom-left corner cell (i, j): – U = 0 (left wall) – V = 0 (bottom wall) We know that Viscous Terms in U-Momentum Equation is given by 2 ∂2U ∂ U + ν ∂x2 ∂y 2 Discretized implicitly (evaluated at time n + 1) using central differences. We can Second Derivative in x-Direction as below: ∂2U ∂x2 n+1 ≈ n+1 n+1 n+1 Ui+1,j − 2Ui,j + Ui−1,j ∆x2 = n+1 n+1 Ui+1,j − 2Ui,j ∆x2 i,j n+1 (since Ui−1,j = 0) Similarly Second Derivative in y-Direction can be written as: ∂2U ∂y 2 n+1 ≈ n+1 n+1 n+1 Ui,j+1 − 2Ui,j + Ui,j−1 ∆y 2 = n+1 n+1 Ui,j+1 − 2Ui,j ∆y 2 i,j 15 n+1 (since Ui,j−1 = 0) Hence after Combining all Implicit Viscous Term we get the final expression as: ! 2 n+1 n+1 n+1 n+1 Ui+1,j − 2Ui,j Ui,j+1 − 2Ui,j ∂2U ∂ U + ≈ν + ν ∂x2 ∂y 2 ∆x2 ∆y 2 We can Extent our methodology to Other Corners too n+1 n+1 – Top-left: Ui−1,j = 0 (left wall), Ui,j+1 is interior n+1 – Top-right/Bottom-right: Adjust for right wall (Ui+1,j = 0) or top/bottom walls – Similar approach applies to all corner cells with appropriate boundary conditions • Central difference for pressure Laplacian For the bottom-left corner cell (i, j): ∂p =0 ∂x ∂p =0 – Bottom wall: ∂y – Left wall: We know that Discretization of Pressure Poisson Equation is given as ∂2p ∂2p + = RHS ∂x2 ∂y 2 Now, solving for x-Direction Laplacian Using central difference with ghost cell we can write: ∂2p ∂x2 From ∂p = 0 at left wall: ∂x ≈ i,j pi+1,j − 2pi,j + pi−1,j ∆x2 pi+1,j − pi−1,j = 0 =⇒ pi−1,j = pi+1,j 2∆x Thus: ∂2p 2(pi+1,j − pi,j ) ≈ 2 ∂x ∆x2 Similarly for y-Direction Laplacian Using central difference with ghost cell we can also write: ∂2p ∂y 2 From ≈ i,j pi,j+1 − 2pi,j + pi,j−1 ∆y 2 ∂p = 0 at bottom wall: ∂y pi,j+1 − pi,j−1 = 0 =⇒ pi,j−1 = pi,j+1 2∆y Thus: ∂2p 2(pi,j+1 − pi,j ) ≈ ∂y 2 ∆y 2 Hence after Combing Pressure Laplacian, we get the final expression as pi+1,j pi,j+1 1 1 ∇2 p ≈ 2 + − p + i,j ∆x2 ∆y 2 ∆x2 ∆y 2 We can also write our results to Other Corners too 16 – Top-left: ∇2 p ≈ 2pi,j−1 − 2pi,j 2pi+1,j − 2pi,j + 2 ∆x ∆y 2 – Top-right/Bottom-right: Similar treatment with appropriate ghost cells • RHS of pressure Poisson equation For the bottom-left corner cell (i, j): ∂p =0 ∂x ∂p – Bottom wall: =0 ∂y – Left wall: The RHS (from continuity enforcement) is written as: ∗ ∗ ∗ ∗ − Ui−1,j − Vi,j−1 Vi,j+1 1 Ui+1,j + RHSi,j = ∆t 2∆x 2∆y On Applying Neumann Boundary Conditions 0.7 From U ∗ Terms ∂p = 0: ∂x Thus: 0.8 From Thus: ∗ ∗ Ui−1,j = −Ui+1,j ∗ ∗ ∗ ∗ ∗ Ui+1,j − (−Ui+1,j ) Ui+1,j Ui+1,j − Ui−1,j = = 2∆x 2∆x ∆x V ∗ Terms ∂p = 0: ∂y ∗ ∗ Vi,j−1 = −Vi,j+1 ∗ ∗ ∗ ∗ ∗ Vi,j+1 − (−Vi,j+1 ) Vi,j+1 Vi,j+1 − Vi,j−1 = = 2∆y 2∆y ∆y Hence, we can write the Final RHS Expression as RHSi,j = 1 ∆t 17 ∗ ∗ Ui+1,j Vi,j+1 + ∆x ∆y QUESTION C: Construct a detailed flow-chart of your fractional step solver which should show all the major steps in the solution procedure. You can use PowerPoint, Illustrator, Inkscape, CorelDraw, etc to draw the flow charts. Include a snapshot to the project report. SOLUTION: Figure 2: Description of the image. 18 QUESTION D: Show that the solver produces results without crashing at very low Re (Re = 1 10). Check if the values of u and v on the boundaries are as expected. Check if the vectors are pointing in the expected direction. SOLUTION: Below is the python code for Re=1 The output for Re=1 is given below which clearly show that code run without crashing at Re=1 19 Now, Python code for Re=10 is given below The output for Re=10 is given as which clearly show that code run without crashing at Re=10 Hence, the solver produces results without crashing at very low Re (Re = 1 10) Now, we can check whether the values of u and v on the boundaries are as expected or not by using the below python code. 20 Below is the python to check for u and v The respective output for the code is which shows that boundary condition fail, u and v on boundaries are not as expected. Velocity Vector Verification The velocity vectors are pointing in the expected direction: • Top boundary (lid): – Vectors point rightward (positive x-direction) – Indicates U = 1, V = 0, matching the boundary condition • Side and bottom walls: – Vectors are zero – Satisfies no-slip condition: U = 0, V = 0 • Interior: – All velocities initially zero – Standard initial condition before flow develops Below is the graph for it 21 QUESTION E: Simulate the lid driven cavity flow at Re = 100 on a uniform 32 × 32 grid. The time step chosen should be large but satisfies the CFL constraint. Simulate the flow until a steady state is reached. • Plot the streamlines of the flow at the steady state. You can use the streamslice function in Matlab and streamplot in Python. • Compare the results with Ghia 1982. Do you see good agreement? Why or Why not? Please discuss. SOLUTION: The python code for Re=100 is given below 22 The output for the code and the plot is given below Plot the streamlines of the flow at the steady state. QUESTION F: Simulate the lid driven cavity flow at Re = 200 on a uniform 32 × 32 grid. The time step chosen should be large but satisfies the CFL constraint. Simulate the flow until a steady state is reached. • Plot the streamlines of the flow at the steady state. • Compare the results with Re = 100. Do you see any difference? Please discuss. 23 SOLUTION: The python code for Re=200 is given below The output for the code is given below 24 The corresponding plot for Re=200 is below Expected Flow Differences (Re = 100 vs Re = 200) • At Re = 100: – Flow is steady and laminar – Single central vortex forms – Streamlines are symmetric – Weak secondary vortices (if any) in corners • At Re = 200: – Still steady (transient at early stages, but eventually stabilizes) – The primary vortex shifts downward slightly – Stronger secondary vortices start to appear in the bottom corners – Velocity gradients near the wall become steeper (thinner boundary layers) Physical Insight / Discussion Points • Inertia increases at higher Re → more pronounced circulation • Viscous diffusion weakens → sharper flow features • Secondary recirculations emerge near corners (especially bottom corners) • The flow field shows higher velocity magnitudes near the lid and stronger shear layers 25 Table 1: Comparison of flow quantities between Re = 100 and Re = 200 Quantity Description Difference to Look For Streamlines Vorticity U-velocity along vertical centerline V-velocity along horizontal centerline Contour of flow direction Curl of velocity field Compare u(x = 0.5, y) Compare v(x, y = 0.5) More complex structure at Re = 200 Sharper and more spread at Re = 200 Peak moves and steepens Deeper trough at Re = 200 QUESTION G: Simulate the lid driven cavity flow at Re = 400 on a uniform 32 × 32 grid. The time step chosen should be large but satisfies the CFL constraint. Simulate the flow until a steady state is reached. • Can you reach a steady state solution? Why or Why not? • What artifacts do you observe during the simulation? Do you have an explanation for that?. SOLUTION: The python code for Re=400 is given below 26 We get the output of this code till Step 9900 27 The corresponding plot is: Steady-State Solution at Re = 400 In principle, yes, a steady-state solution can be reached for lid-driven cavity flow at Re = 400. At this moderate Reynolds number, the flow remains laminar and steady. However, the ability to actually reach that state numerically depends on a few critical aspects: • Time step size: While you mentioned using a large time step that satisfies the CFL condition, even CFL-compliant large steps can introduce numerical instability or inaccuracy, especially in a semiimplicit scheme. • Solver accuracy and convergence criteria: If your pressure-velocity coupling or Poisson solver isn’t tight enough (e.g., loose residual tolerances), the solver might stagnate before reaching the true steady state. Observed Artifacts During Simulation From the plot: • The velocity field looks incorrect. The only non-zero velocities appear to be along the top lid, with no apparent recirculating flow inside the cavity. • The interior flow is practically zero, suggesting that the lid motion is not properly propagating momentum into the domain. Possible Explanations for These Artifacts Here are several likely causes: 28 (a) Incorrect Boundary Conditions or Initialization • If you only set velocity at the top boundary and didn’t let momentum diffuse into the domain properly (via convection + diffusion), the flow will remain stagnant. • This could also happen if you don’t apply no-slip conditions on the other three walls correctly. (b) Projection Step (Pressure Correction) Issue • If you’re using a fractional step method, any error in solving the pressure Poisson equation (e.g., wrong Neumann/Dirichlet BCs or loose tolerance) can prevent correct enforcement of incompressibility. • This would suppress any vorticity or recirculation from forming in the interior. (c) Discretization Error / Numerical Damping • On a 32×32 grid, if the discretization is too diffusive (especially with large time steps), it could smear out the vortex development. • This is especially true if you use first-order upwinding for convection terms. (d) Time Step Still Too Large (Borderline CFL) • Large ∆t can cause numerical damping of transient features before the vortex structure develops. • A smaller ∆t might be required to capture transient evolution properly and allow convergence. QUESTION H: Extra Credit (10% of the total grade): Try to simulate the flow at higher Reynolds number. Do you see the unexpected checkerboard pattern? Explain what happened. To receive this extra credit, you can use any of the methods discussed in class to fix this. A good suggestion is to use a separate set of velocities (U, V ) on the boundaries that are designed to satisfy the divergence free condition. See Zang, Yan, Robert L. Street, and Jeffrey R. Kosoff. “A non-staggered grid, fractional step method for time-dependent incompressible Navier-Stokes equations in curvilinear coordinates.” Journal of Computational Physics 114.1 (1994): 18-33. • Demonstrate that your code can solve the flow at Re = 1, 000. Compare the results with Ghia 1982. • You have to show very good agreement to receive the credit. SOLUTION: When we run the simulation at higher Reynold number such as 800 , 100, etc, we see the following observation • The flow becomes more complex and faster. • Stronger main vortices and corner eddies develop. • Numerical errors grow more easily, especially near boundaries. 29 Yes, We can see the unexpected checkerboard pattern, The reasons for this is written below • It often occurs when using a collocated grid (pressure and velocity stored at the same location). • The solver might produce a mathematically valid but physically incorrect alternating pattern. • The discretization decouples pressure and velocity. • The equations allow for alternating values (e.g., +1, −1, +1, −1, . . . ), which don’t represent real fluid flow. • This satisfies the discrete divergence condition but creates non-physical pressure and velocity fields. We can use the following method to fix it 1. Use a Smooth, Divergence-Free Lid Velocity • Instead of a flat profile (u = 1), use: u(x, y = 1) = 16x2 (1 − x)2 , v(x, y = 1) = 0 • Ensures zero velocity at the corners. • Keeps the flow divergence-free at the boundary. 2. Use Momentum Interpolation (For Collocated Grids) • Interpolate velocities using both velocity and pressure: ui+1/2 = ∆t pi+1 − pi ui + ui+1 − 2 2ρ ∆x • This properly couples pressure and velocity, avoiding the checkerboard issue. 3. Consider Using a Staggered Grid • Pressure is stored at cell centers, while velocity components are stored at cell faces. • This natural separation helps maintain coupling and prevents checkerboard artifacts. Demonstration of the code at Re=1,000 At Re = 1000? • A large central vortex forms. • Two smaller corner vortices appear at the bottom left and right. • Flow looks smooth and realistic — not noisy or unstable. Set the Reynolds number and calculate viscosity like this: Re = 1000, Ulid = 1.0, L = 1.0, 30 ν= Ulid · L = 0.001 Re We can see the following observation and see that, it is working at that Reynold number • The main vortex is stable and smooth. • No unphysical patterns like checkerboarding. • The simulation doesn’t crash or blow up. • Your velocity profile is close to benchmark data (e.g., Ghia et al.). 31