Benjamin Seibold Applied Mathematics Massachusetts Institute of Technology www-math.mit.edu/~seibold firstname.lastname@example.org March 31, 2008
On the following pages you ﬁnd a documentation for the Matlab program mit18086 navierstokes.m, as it is usedfor Course 18.086: Computational Science and Engineering II at the Massachusetts Institute of Technology. The code can be downloaded from the Computational Science and Engineering web page http://www-math.mit.edu/cse. It is also linked on the course web page http://www-math.mit.edu/18086. This code shall be used for teaching and learning about incompressible, viscous ﬂows. It is an example of asimple numerical method for solving the Navier-Stokes equations. It contains fundamental components, such as discretization on a staggered grid, an implicit viscosity step, a projection step, as well as the visualization of the solution over time. The main priorities of the code are 1. Simplicity and compactness: The whole code is one single Matlab ﬁle of about 100 lines. 2. Flexibility: The codedoes not use spectral methods, thus can be modiﬁed to more complex domains, boundary conditions, and ﬂow laws. 3. Visualization: The evolution of the ﬂow ﬁeld is visualized while the simulation runs. 4. Computational speed: Full vectorization and pre-solving the arising linear systems in an initialization step results in fast time stepping. The code provides: • variable box sizes, grid resolution,time step, and Reynolds numbers • implicit time stepping for the ﬂuid viscosity 1
• visualization of pressure ﬁeld, velocity ﬁeld and streamlines • automatic transition from central to donor-cell discretization for the nonlinear advection part • fast computation of the time-dependent solution for small to moderate mesh sizes The code does not provide: • 3d • unstructured meshes or complexgeometries • time-dependent geometries • adaptivity, such as local mesh reﬁnement, time step control, etc. • non-constant viscosity or density • higher order time stepping • turbulence models • time and memory eﬃciency for large computations The code qualiﬁes as a basis for modiﬁcations and extensions. The following extensions have already been applied to the code by MIT students and other users: •external forces • inﬂow and outﬂow boundaries • addition of a drag region • geometry change to a backwards facing step • time dependent geometries • coupling with an advection-diﬀusion equation
Incompressible Navier-Stokes Equations
We consider the incompressible Navier-Stokes equations in two space dimensions 1 (uxx + uyy ) Re 1 vt + py = −(uv)x − (v 2 )y + (vxx + vyy ) Re ux + vy = 0 ut+ px = −(u2 )x − (uv)y + 2 (1) (2) (3)
on a rectangular domain Ω = [0, lx ]×[0, ly ]. The four domain boundaries are denoted North, South, West, and East. The domain is ﬁxed in time, and we consider no-slip boundary conditions on each wall, i.e. u(x, ly ) = uN (x) u(x, 0) = uS (x) u(0, y) = 0 u(lx , y) = 0 v(x, ly ) = 0 v(x, 0) = 0 v(0, y) = vW (y) v(lx , y) = vE (y) (4) (5) (6) (7)
Aderivation of the Navier-Stokes equations can be found in . The momentum equations (1) and (2) describe the time evolution of the velocity ﬁeld (u, v) under inertial and viscous forces. The pressure p is a Lagrange multiplier to satisfy the incompressibility condition (3). Note that the momentum equations are already put into a numerics-friendly form. The nonlinear terms on the right hand sideequal (u2 )x + (uv)y = uux + vuy (uv)x + (v )y = uvx + vvy
which follows by the chain rule and equation (3). The above right hand side is often written in vector form as (u · )u. We choose to numerically discretize the form on the left hand side, because it is closer to a conservation form. The incompressibility condition is not a time evolution equation, but an algebraic condition....