CS205b/CME306Lecture 171 Incompressible FlowSupplementary Reading: Osher and Fedkiw, §18.1, §18.2Recall that the system of equations we must solve for incompr essib le flow is∇ · u = 0 (1)ρt+ u · ∇ρ = 0 (2)ut+ u · ∇u +∇pρ= g. (3)(4)1.1 Projection MethodIn order to u pdate the velocity, we use the projection method due to Chorin [1]. The projectionmethod is applied by first computing an intermediate velocity field u⋆ignoring the pressure term,u⋆− un∆t+ (un· ∇)un= gn, (5)and then computing a divergence f ree velocity field un+1,un+1− u⋆∆t+∇pn+1ρn+1= 0, (6)using th e pressure as a correction. Note that combining equations 5 and 1.4 to eliminate u⋆resultsin equation 1.4.Taking the divergence of equation 1.4 results in∇ ·∇pn+1ρn+1=∇ · u⋆∆t(7)after setting ∇·un+1= 0, i.e. after assuming that the new velocity field is divergence free. Equation7 defines th e p ressure in terms of the value of ∆t used in equation 5. Defining a scaled pressu re ofp⋆= p∆t leads toun+1− u⋆+∇p⋆ρn+1= 0and∇ ·∇p⋆ρn+1= ∇ · u⋆1in place of equations 1.4 and 7 where p⋆does not d epend on ∆t. When the density is spatiallyconstant, we can define ˆp = p∆t/ρ leading toun+1− u⋆+ ∇ˆp = 0and∆ˆp = ∇ · u⋆(8)where ˆp does not depend on ∆t or ρ.This method u tilizes the Helmholtz-Hodge decomposition of the vector field u⋆,u⋆= un+1+ ∇ˆp.In general, the Helmholtz-Hodge decomposition of a vector field expresses the vector field as adivergence free vector field plus th e gradient of a scalar field.1.2 MAC GridHarlow and Welch [2] proposed the use of a special grid for incompressible fl ow computations. Thisspecially defined grid decomposes the computational domain into cells with velocities defined onthe cell f aces and scalars defined at cell centers. That is, in 2D, pi,j, ρi,jare defined at cell centerswhile ui±12,jand vi,j±12are defined at the appropriate cell f aces.Equation (2) is solved by first defining the cell center velocities with s imple averagingui,j=ui−12,j+ ui+12,j2vi,j=vi,j−12+ vi,j+122Then the spatial derivatives are evaluated in a straightforward manner, for example us ing 3rd orderaccurate Hamilton-Jacobi ENO. The temporal derivative can be evaluated with a TVD RK scheme.In order to update the velocity based on equation (1.4), we first need u and v at all the cellfaces. Again, we obtain th e values by simple averaging. For example,vi−12,j=14vi−1,j−12+ vi−1,j+12+ vi,j−12+ vi,j+12.Similarly, to get u values on the v faces, we compute the averageui,j−12=14ui−12,j−1+ ui+12,j−1+ ui−12,j+ ui+12,j.1.3 Discretizing Divergence of VelocitySupplementary Reading: Osher and Fedkiw, §18.3, §23.12Here we discuss the discretization in 8 of the term ∇ · u⋆. Since we are solving for the pres-sure, which on the MAC grid lives at the cell centers, we need to discretize the term at the cellcenters. We have that(∇ · u⋆)i,j=u⋆x+ v⋆yi,j=u⋆i+12,j− u⋆i−12,j∆x+v⋆i,j+12− v⋆i,j−12∆y+ O(∆x2) + O(∆y2)So we have used the intermediate velocity stored at the cell faces to get a second order accurateapproximation to ∇ · u⋆at the cell centers.1.4 Semi-Lagrangian Velocity AdvectionEvolving the momentum equation is done by first advecting velocities and applying body forceswith to obtain u⋆. Then, equation is evaluated by solving the elliptic Poisson’s equation, followedby applying the resulting pressure to obtain a final velocity un+1. There is no time step restrictionfor second step, so the only CFL r estriction is on the velocity advection. Therefore, if we usethe semi-Lagrangian method for the velocity advection we can eliminate the remaining time steprestriction. For u⋆, the method isu⋆j= unxj− unj∆tFor v⋆we must also account for gravity, so we havev⋆j= vn(xj− unj∆t) − ∆t gwhere we are computingvt+ u · ∇v = −g.This is a Godunov splitting, which is firs t order accurate.1.5 Algorithm OverviewWe now have the following steps in updating the velocity field for incompressible flow using theprojection method:1. Compute the intermediate velocity field u⋆(at cell faces)u⋆− un∆t+ u · ∇u = g (9)2. Solve an elliptic equation for the pressure (at cell centers)∇ ·1ρ∇p=∇ · u⋆∆t(10)3. Compute the divergence free velocity field un+1(at cell faces)un+1− u⋆∆t+∇pρ= 0 (11)3We can multiply p by ∆t, to get p⋆= p∆t, and rewrite steps 2 and 3 as2a.∇ ·1ρ∇p⋆= ∇ · u⋆(12)3a.un+1− u⋆+∇p⋆ρ= 0 (13)If ρ is constant, then we can move it under the ∇ operator and defining ˆp =p∆tρ, we can rewritesteps 2 and 3 as2b.∆ˆp = ∇ · u⋆(14)3b.un+1− u⋆+ ∇ˆp = 0 (15)References[1] Chorin, A., A Numerical Method for Solving Incompressible Viscous Flow Problems, J. Com-put. Phys. 2, 12-26 (1967).[2] Harlow, F. and Welch, J., Numerical Calculation of Time-Dependent Viscous IncompressibleFlow of Fluid with a Free Surface, The Physics of Fluids 8, 2182-2189