Lecture 11Monday, May 9, 2005Supplementary Reading: Osher and Fedkiw, §18.2In the previous lecture we started looking at incompressible flow, where wehave the incompressibility assumption, ∇ ·~V = 0. Under this assumption,we have a relatively smooth flow, with no shocks or rarefactions (although wemay still have contact discontinuities). The equations for conservation of mass,momentum, and energy, assuming incompressibility, are given byρt+~V · ∇ρ = 0 (1)~Vt+~V · ∇~V +∇pρ= 0 (2)et+~V · ∇e = 0 (3)Body forces, e.g. gravity, are added to the RHS of the momentum equation, sothat it becomes~Vt+~V · ∇~V +∇pρ= ~g (4)where ~g = (0, g, 0).1 MAC GridHarlow and Welch [1] proposed the use of a special grid for incompressibleflow computations. This specially defined grid decomposes the computationaldomain into cells with velocities defined on the cell faces and scalars defined atcell centers. That is, in 2D, pi,j, ρi,jare defined at cell centers while ui±12,jandvi,j±12are defined at the appropriate cell faces.Equation (1) is solved by first defining the cell center velocities with simpleaveragingui,j=ui−12,j+ ui+12,j2vi,j=vi,j−12+ vi,j+1221Then the spatial derivatives are evaluated in a straightforward manner, for ex-ample using 3rd order accurate Hamilton-Jacobi ENO. The temporal derivativecan be evaluated with a TVD RK scheme.In order to update the velocity based on equation (2), we first need u andv at all the cell faces. Again, we obtain the values by simple averaging. Forexample,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.2 Projection MethodIn order to update the velocity, we use the projection method due to Chorin [2].The projection method is applied by first computing an intermediate velocityfield~V?ignoring the pressure term,~V?−~Vn4t+~Vn· ∇~Vn= ~gn, (5)and then computing a divergence free velocity field~Vn+1,~Vn+1−~V?4t+∇pn+1ρn+1= 0, (6)using the pressure as a correction. Note that combining equations 5 and 6 toeliminate~V?results in equation 4 exactly.Taking the divergence of equation 6 results in∇ ·∇pn+1ρn+1=∇ ·~V?4t(7)after setting ∇ ·~Vn+1to zero, i.e. after assuming that the new velocity fieldis divergence free. Equation 7 defines the press ure in terms of the value of 4tused in equation 5. Defining a scaled pressure of p?= p4t leads to~Vn+1−~V?+∇p?ρn+1= 0 (8)and∇ ·∇p?ρn+1= ∇ ·~V?(9)in place of equations 6 and 7 where p?does not depend on 4t. When the densityis spatially constant, we can define ˆp = p4t/ρ leading to~Vn+1−~V?+ ∇ˆp = 0 (10)2and4ˆp = ∇ ·~V?(11)where ˆp does not depend on 4t or ρ.This method utilizes the Helmholtz-Hodge decomposition of the vector field~V?,~V?=~Vn+1+ ∇ˆp.In general, the Helmholtz-Hodge decomposition of a vector field expresses thevector field as a divergence free vector field plus the gradient of a scalar field.3 Laplace EquationIn order to discuss the solution of (11), let us first consider solving the Laplaceequation.In 1D, the equation is give bypxx= 0The solution is simply a linep = ax + b (12)The values of the constant a and b are determined by boundary conditions.Assume that the domain is the interval [0, 1]. We may have Dirichlet boundaryconditions, where the value of the function p is given at the boundary. Forexample,p(0) = p0p(1) = p1Plugging the b oundary conditions in the the equation(12), we getp(0) = b = p0p(1) = a + b = p1⇒ a = p1− p0so the coefficients a and b are uniquely determined. Alternatively, Neumannboundary conditions specify the value of pxat the boundary. For example,px(0) = 0 ⇒ a = 0.This gives us a family of lines with slope 0. To find b, we would need anotherpiece of information. A Dirichlet boundary condition would pick out one of thelines with slope 0, thus determining the solution. But observe that specifyingtwo Neumann conditions could lead to no solution. For example,px(0) = 0px(1) = 13These two boundary conditions are inconsistent, hence there is no solution.Another example ispx(0) = 0px(1) = 0In this case, the given boundary conditions are consistent, but incomplete. Westill do not have enough information to identify a unique solution. The aboveexamples illustrate the fact that in 1D, for the Laplace equation, we can deter-mine the solution if we have two Dirichlet boundary conditions or one Neumannand one Dirichlet boundary condition, but will have either no solution or an un-derdetermined solution in the case of two Neumann boundary conditions.For a discretization in 2D, we can determine the solution if we have ei-ther Dirichlet boundary conditions at each point on the boundary, or Neumannboundary conditions at each point except for at least one where we have aDirichlet boundary condition.References[1] Harlow, F. and Welch, J., Numerical Calculation of Time-Dependent Vis-cous Incompressible Flow of Fluid with a Free Surface, The Physics ofFluids 8, 2182-2189 (1965).[2] Chorin, A., A Numerical Method for Solving Incompressible Viscous FlowProblems, J. Comput. Phys. 2, 12-26