A compact and fast matlab code solving the incompressible navier-stokes equations on rectangular domains

A compact and fast Matlab code solving the incompressible Navier-Stokes equations on rectangular domains mit18086 navierstokes.m
Benjamin Seibold Applied Mathematics Massachusetts Institute of Technology www-math.mit.edu/~seibold seibold@math.mit.edu March 31, 2008

1

Introduction

On the following pages you find 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 flows. 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 file of about 100 lines. 2. Flexibility: The codedoes not use spectral methods, thus can be modified to more complex domains, boundary conditions, and flow laws. 3. Visualization: The evolution of the flow field 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 fluid viscosity 1

• visualization of pressure field, velocity field 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 refinement, time step control, etc. • non-constant viscosity or density • higher order time stepping • turbulence models • time and memory efficiency for large computations The code qualifies as a basis for modifications and extensions. The following extensions have already been applied to the code by MIT students and other users: •external forces • inflow and outflow boundaries • addition of a drag region • geometry change to a backwards facing step • time dependent geometries • coupling with an advection-diffusion equation

2

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 fixed 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 [2]. The momentum equations (1) and (2) describe the time evolution of the velocity field (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
2

(8) (9)

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....