CS205b/CME306Lecture 191 ViscosityWe now focus on the discretization of the viscosity term in the Navier-Stokes equations. Typicallythe inviscid equations are called the Euler equations while the viscous equations are called theNavier-Stokes equations.For incompressible flow with nonzero viscosity we still have the same equation for conservationof mass. It is given byρt+ u · ∇ρ = 0.However, the momentum equation (in 2D) becomes(ut+ u · ∇u +pxρ=(2µux)x+(µ(uy+vx))yρvt+ u · ∇v +pyρ=(µ(uy+vx))x+(2µvy)yρ− g(1)where we have added the viscosity terms to the RHS of the equation. In vector form, this is canbe written asut+ u · ∇u +∇pρ= g +(∇ · τ )Tρ.Now consider the special case where µ = constant in (1). In that case we can simplify theviscosity term on the RHS as follows.(2µux)x+ (µ(uy+ vx))yρ=2µuxx+ µuy y+ µvxyρ=µ(uy y+ uxx)ρ+µ(uxx+ vxy)ρ=µ(uy y+ uxx)ρ+µ(ux+ vy)xρ=µ(uy y+ uxx)ρ+ 0=µρ∆u1(µ(uy+ vx))x+ (2µvy)yρ=µuy x+ µvxx+ 2µvy yρ=µ(vxx+ vy y)ρ+µ(vy y+ uxy)ρ=µ(vxx+ vy y)ρ+µ(vy+ ux)yρ=µ(vxx+ vy y)ρ+ 0=µρ∆vTherefore for µ = constant, the equations (1) become(ut+ u · ∇u +pxρ=µρ∆uvt+ u · ∇v +pyρ=µρ∆v − g(2)1.1 DiscretizationIn the projection method for incompressib le flow the viscosity term is included in the computationof u⋆, the intermediate velocity field. That is, the steps in the projection method become1. Compute the intermediate velocity field u⋆u⋆− un∆t+ un· ∇un=(∇ · τ )Tρ+ g2. Solve an elliptic equation for the pressure∆ˆp = ∇ · u⋆3. Compute the divergence free velocity field un+1un+1− u⋆+ ∇ˆp = 0where we have again assume that ρ = constant, and set ˆp =p∆tρ.Next we will discretize the viscous terms in (2). Since we are using a MAC grid and u⋆isdefined at the cell walls, we need the viscous term discretized at the cell walls. We approximatethe Laplacian of u at the grid point i +12, j as(∆un)i+12,j≈uni−12,j− 2uni+12,j+ uni+32,j∆x2+uni+12,j−1− 2uni+12,j+ uni+12,j+1∆y2This is a second order central difference approximation. The problem with this approximationis that it requires that ∆t ∼ ∆x2for stability. This is a severe restriction on the time step and wewou ld like to avoid it. One solution, due to Kim and Moin, is to treat the viscosity implicitly. Sofor step 1 in the projection method, we solve the equationu⋆− un∆t+ un· ∇un=(∇ · τ⋆)Tρ+ g2The term un· ∇unis still treated the same as before. Th en the terms at time step n will be onthe RHS, while the ⋆ terms are on the LHS. In the case of constant µ, we get a decoupled linearsystem of the formA1u = b1A2v = b2Another possibility is to use trapezoidal ru leu⋆− un∆t+ un· ∇un=(∇ · τ⋆)T+ (∇ · τn)T2ρ+ gOne problem in incompressible flow is that the numerical viscosity may be larger than thephysical viscosity. We want the numerical viscosity arising from the d iscretization of the u · ∇uterm to be smaller than the physical viscosity∇·τρ.Recall the first order upwind discretization of the advection equationut+ ux= 0.The discretization isut+ui− ui−1∆x= 0.⇒ ut+ui−ui− ∆x(ux)i+∆x22(uxx)i+ O(∆x3)∆x= 0⇒ ut+ (ux)i−∆x2(uxx)i= O(∆x2)⇒ ut+ (ux)i=∆x2(uxx)i+ O(∆x2)Now suppose you want to solveut+ ux= µuxx.From the above, we see that using a first order upwind discretization for uxour modified equationwill beut+ ux=µ +∆x2uxx.µ is the real viscosity and∆x2is the numerical viscosity. One of the big problems with solvingNavier-Stokes is that the numerical viscosity is often larger than the real viscosity.2 VorticityHere we describe a method to counteract the numerical dissipation that damps out many interestingfeatures in the flow.Taking the curl of the momentum equationut+ u · ∇u +∇pρ= g3givesΩt+ u · ∇Ω − Ω · ∇u −1ρ2∇p × ∇ρ = ∇ × gwhereΩ = ∇ × u.In 2D,Ω =i j k∂∂x∂∂y∂∂zu v 0=−∂∂zv∂∂zu∂∂xv −∂∂yuSince∂∂zu =∂∂zv = 0we haveΩ =00vx− uy=00ΩSo this is particularly nice in 2D as we get one scalar equation for Ω (in 3D we still get a 3-vector).Since Ω will be either positive or negative, the vorticity vector Ω is pointing either into or out ofthe x − y plane. Vorticity can be thought of as a p addle wheel which is trying to spin the flow.The direction of the spinning depends on the sign of Ω .Some points of interest regarding vorticity are• Vorticity is conserved.• Vorticity stays confinedin high Reynolds number flows.Here we discuss a simple turbulence model due to Steinhoff. First we compute vorticity locationvectorsn =∇kΩkk∇kΩkk.Then we compute the paddle wheel force asf = n × Ω.Steinhoff’s idea was to add a forcing term to the momentum equationsut+ u · ∇u +∇pρ= g + ǫ∆xfIt is interesting to note that if you linearize the forcing term, it looks like