CS205b CME306 Lecture 17 1 1 1 Incompressible Flow MAC Grid Supplementary Reading Osher and Fedkiw 18 1 18 2 Recall that the system of equations we must solve for incompressible flow is u 0 t u 0 p ut u u g 1 2 3 4 Harlow and Welch 2 proposed the use of a special grid for incompressible flow computations This specially defined grid decomposes the computational domain into cells with velocities defined on the cell faces and scalars defined at cell centers That is in 2D pi j i j are defined at cell centers while ui 1 j and vi j 1 are defined at the appropriate cell faces 2 2 Equation 2 is solved by first defining the cell center velocities with simple averaging ui j vi j ui 1 j ui 1 j 2 2 2 vi j 1 vi j 1 2 2 2 Then the spatial derivatives are evaluated in a straightforward manner for example using 3rd order accurate Hamilton Jacobi ENO The temporal derivative can be evaluated with a TVD RK scheme In order to update the velocity based on equation 3 we first need u and v at all the cell faces Again we obtain the values by simple averaging For example 1 vi 1 j 1 vi 1 j 1 vi j 1 vi j 1 vi 1 j 2 2 2 2 2 4 Similarly to get u values on the v faces we compute the average 1 ui 1 j 1 ui 1 j 1 ui 1 j ui 1 j ui j 1 2 2 2 2 2 4 The term MAC Grid stands for Marker And Cell It refers to what the discretization was first used for rather than describing the discretization itself 1 1 2 Discretizing Divergence of Velocity Supplementary Reading Osher and Fedkiw 18 3 23 1 Here we discuss the discretization in 8 of the term u Since we are solving for the pressure which on the MAC grid lives at the cell centers we need to discretize the term at the cell centers We have that u i j u x vy i j u i 1 j u i 1 j 2 2 x vi j 1 v i j 1 2 2 y O x2 O y 2 So we have used the intermediate velocity stored at the cell faces to get a second order accurate approximation to u at the cell centers 1 3 Semi Lagrangian Velocity Advection Evolving the momentum equation 3 is done by first advecting velocities and applying body forces with 6 to obtain u Then equation 6 is evaluated by solving the elliptic Poisson s equation followed by applying the resulting pressure to obtain a final velocity un 1 There is no time step restriction for second step so the only CFL restriction is on the velocity advection Therefore if we use the semi Lagrangian method for the velocity advection we can eliminate the remaining time step restriction For u the method is u j un xj unj t For v we must also account for gravity so we have vj v n xj unj t t g where we are computing vt u v g This is a Godunov splitting which is first order accurate 1 4 Projection Method In order to update the velocity we use the projection method due to Chorin 1 The projection method is applied by first computing an intermediate velocity field u ignoring the pressure term u un un un gn t 5 and then computing a divergence free velocity field un 1 un 1 u pn 1 n 1 0 t 6 using the pressure as a correction Note that combining equations 5 and 6 to eliminate u results in equation 3 2 Taking the divergence of equation 6 results in n 1 p u n 1 t 7 after setting un 1 0 i e after assuming that the new velocity field is divergence free equation 7 defines the pressure in terms of the value of t used in equation 5 Defining a scaled pressure of p p t leads to p un 1 u n 1 0 and p n 1 u in place of equations 6 and 7 where p does not depend on t When the density is spatially constant we can define p p t leading to un 1 u p 0 and p u 8 where p does not depend on t or This method utilizes 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 a divergence free vector field plus the gradient of a scalar field 1 5 Computing Boundary Conditions We have shown how to discretize the Poisson equation and handle Dirichlet and Neumann boundary conditions We have not yet said much about how to obtain these boundary conditions A typical scenerio that leads to a Dirichlet boundary condition is a free surface such as the surface of the water in a glass In this case the pressure will be 1 atm the pressure that the air applies to the surface of the water The other type of boundary condition that occurs in the example of a glass of water is the boundary condition between the water and the walls of the container This boundary condition may be described by requiring the velocity component at the cell face of the boundary to be equal to the velocity of the wall itself If the walls of the container are stationary then this velocity component will be zero Unfortunately this is a condition on velocity not a condition on pressure as is required by the Poisson discretization To obtain this boundary condition we begin with 6 Lets also assume that the boundary condition is on the side of the container so that we are considering the x direction so that we should consider un 1 u px 0 t 3 For this face we have un 1 uBC which simply states that we should be computing a velocity that agrees with the velocity uBC of the container uBC u px 0 t Solving this for the pressure derivative yields px uBC u t This procedure may be simplified somewhat by observing that we may simply enforce the boundary conditions on u so that it need not be computed but can just be assigned u uBC Then the corresponding pressure condition becomes px 0 1 6 Algorithm Overview We now have the following steps in updating the velocity field for incompressible flow using the projection method 1 Compute the intermediate velocity field u at cell faces u un u u g t 9 2 Solve an elliptic equation for the pressure at cell centers u 1 p t 3 Compute the divergence free velocity field un 1 10 at cell faces un 1 u p 0 t 11 We can multiply p by t to get p p t and rewrite steps 2 and 3 as 2a u 12 un 1 u p 0 13 1 p 3a If is constant then we …