CS205b/CME306Lecture 161 Incompressible Flow1.1 Laplace Equation in 1DSupplementary Reading: Osher and Fedkiw, §18.1, §18.2Recall that the system of equations we must solve for incompr essib le flow is∇ · u = 0ρt+ u · ∇ρ = 0ut+ u · ∇u +∇pρ= g.The Laplace equation in 1D is given bypxx= 0.The solution is simply a linep = ax + b. (1)The values of the constant a and b are determined by boundary conditions. Assume that thedomain is the interval [0, 1]. We may have Dirichlet boundary conditions, where the value of thefunction p is given at the boundary. For example,p(0) = p0p(1) = p1.Plugging the boundary conditions in the the equation(1), we getp(0) = b = p0p(1) = a + b = p1⇒ a = p1− p0so the coefficients a and b are uniquely determined. Alternatively, Neumann boundary conditionsspecify the value of pxat the boundary. For example,px(0) = 0 ⇒ a = 0.This gives us a family of lines with s lope 0. To find b, we would need another piece of information.A Dirichlet bou ndary condition would pick out one of the lines with slope 0, thus determining1the solution. But observe that specifying two Neumann conditions could lead to no solution. Forexample,px(0) = 0 px(1) = 1.These two boundary conditions are inconsistent, hence there is no solution. Another example ispx(0) = 0 px(1) = 0.In this case, the given boundary conditions are consistent, but incomplete. We still do not haveenough information to identify a unique solution. The above examples illustrate the fact that in1D, for the Lap lace equation, we can determine th e solution if we have two Dirichlet boundaryconditions or one Neumann and one Dirichlet boundary condition, but will have either no solutionor an un derdetermined solution in the case of two Neumann boundary conditions.1.2 Discretizing Laplacian of PresssureWe need to numerically solving Poisson’s equationpxx= f(x).We will also need the gradient to apply the pressure. We use second order central differencing forboth. At each cell face, we approximate th e pressure gradient with(px)i+1/2=pi+1− pi∆x+ O(∆x2).From this we use central differencing again to express the Laplacian at each grid node(pxx)i=(px)i+1/2− (px)i−1/2∆x+ O(∆x2)=pi+1− 2pi+ pi−1∆x2+ O(∆x2)= fi.The result is a coupled linear system that we need to solve in order to determine p on the entiredomain. However, we cannot write this equation as is for the grid points near the boundary since itwill involve points outside of the domain. For example, assume that our domain is the interval [0, 1]and that we have grid points 0, 1, . . . , M, M + 1 uniformly spaced on the domain. The equation forp1isp2− 2p1+ p0∆x2= f1.If we have a Dirichlet boundary condition specified on the left of the domainp0= β,then the equation for p1becomesp2− 2p1∆x2= f1−β∆x2.If we have a Neumann boundary condition specified at the half grid point12(px)12= α,2we write the equation for p1asp2−p1∆x−p1−p0∆x∆x= f1.Since(px)12=p1− p0∆x+ O(∆x2),the equation for p1becomesp2− p1∆x2= f1+α∆x.Let’s look at the matrix equation for the case where we have two Dirichlet boundary conditions.−2 11 −2 1.........1 −2 1.........1 −2 11 −2p1p2...pi...pM− 1pM=∆x2f1− p0∆x2f2...∆x2fi...∆x2fM− 1∆x2fM− pM+ 1The matrix is sym metric negative definite. This is advantageous because there are fast linearsolvers for such systems, e.g. th e conju gate gradients method.In the case with two Neumann boundary conditions, the matrix equation is−1 11 −2 1.........1 −2 1.........1 −2 11 −1p1p2...pi...pM− 1pM=∆x2f1+ ∆x(px)12∆x2f2...∆x2fi...∆x2fM− 1∆x2fM− ∆x(px)M+12.Notice that the matrix has changed. In particular, it is singular since it has a non-empty null spacewhich is spanned by the vector (1, . . . , 1)T. T his is problematic, but workable. It can be solved forp up to a constant, since for any solution, ~p, ~p + c(1, . . . , 1)Tis also a solution.In multiple dimension Poisson’s equation is∆p = f.In 2D the equation ispxx+ py y= f.We again use the second order accurate central differencing to obtain the gradient components atthe cell faces(px)i+1/2,j=pi+1,j− pi, j∆x+ O(∆x2)(py)i,j+1/2=pi,j+1− pi, j∆y+ O(∆y2)3from which we apply central d ifferencing again to obtain the Laplacian(∆p)i,j=(px)i+1/2,j− (px)i−1/2,j∆x+(py)i,j+1/2− (py)i,j−1/2∆y+ O(∆x2) + O(∆y2)=pi+1,j− 2pi,j+ pi−1,j∆x2+pi,j+1− 2pi,j+ pi,j−1∆y2+ O(∆x2) + O(∆y2)= fi,j.In 2D we need boundary conditions specified around the entire domain. If at least one boundarycondition is Dirichlet, then the resulting matrix will be a banded s ymmetric positive definite matrix.We can us e an iterative solver such as p reconditioned conjugate gradients. If all the boundaryconditions are Neumann, then the matrix will have a null space, and we must ensure we have acompatible system.If the density is s patially varying, then we must use a discretization with variable coefficients.To simplify the notation slightly, let β =1ρ. Then, we can write∇ ·1ρ∇pi,j= (∇ · β∇p)i,j=βi+1/2,j(px)i+1/2,j− βi−1/2,j(px)i−1/2,j∆x+ O(∆x2)+βi,j+1/2(py)i,j+1/2− βi,j−1/2(py)i,j−1/2∆y+ O(∆y2)= fi,j.Note that this will require densities at the cell walls.1.3 Compatibility ConditionPoisson’s equation with all Neumann boundary conditions must satisfy a compatibility conditionfor a solution to exist. The problem is given by∆p = f in Ω∇p · n = g on ∂Ωwhere n is the unit normal to the boun dary. From the equation we have the relationsZΩf dV =ZΩ∆p dV =ZΩ∇ · ∇p dV =Z∂Ω∇p · n dS =Z∂Ωg dSwhere the third equality follows from the divergence theorem. The compatibility condition isZΩf dV =Z∂Ωg dS.The right hand side f will be of the form f = ∇ · u⋆, and g = 0. Therefore, the compatibilitycondition isZΩ∇ · u⋆dV =Z∂Ωu⋆· n dS = 0where the first equality follows from the divergence theorem. This condition needs to be satisfiedwhen specifying the boundary