0.10 Numerical simulation

The classical exact solutions of the Navier-Stokes equations are valuable for the insight they give and as an exact result that any approximate method should be able to duplicate to an acceptable precision. However, if we limit ourselves to exact solutions or even approximate perturbation or series solutions, a vast number of important engineering problems will be beyond reach. To meet this need, a number of computational fluid mechanics (CFM) codes have become commercially available. A student could use the CFM code as a virtual experiment to learn about fluid mechanics. While this may be satisfactory for someone who will become a design engineer, this is not sufficient for the researcher who may wish to solve a problem that has not been solved before. The user usually has to treat the CFM code as a "black box" and accept the result of the simulation as correct. The researcher must always be questioning if any computed result is valid and if not, what must be done to make it valid. Therefore I believe Ph.D. students should know how to develop their own numerical simulators rather than be just a user.

We will start with a working FORTAN and MATLAB code for solving very simple, generic problems in 2-D. The student should be able to examine each part of the code and understand everything with the exception of the algorithm for solving the linear system of equations. Then the student will change boundary conditions, include transient and nonlinear capabilities, include curvilinear coordinates, and compute pressure and stresses.

The stream function - vorticity method

Two-dimensional flow of an incompressible, Newtonian fluid can be formulated with the stream function and vorticity as dependent variables.

We will begin developing a generic code by assuming steady, creeping (zero Reynolds number) flow with specified velocity on the boundaries in a rectangular domain. The PDE are now,

This is a pair of linear, elliptic PDEs. The boundary conditions with specified velocity are

The boundary conditions at a plane of symmetry are

The PDE and definitions are made dimensionless with respect to reference quantities such that the variables are of the order of unity.

We have four unspecified dimensionless groups and two unspecified reference quantities. We can specify two of the groups to equal unity and the other two are equal to the aspect ratio or its square.

The PDE and definitions with the * dropped are now as follows.

Finite difference approximation

The PDE and BC will be solved using a finite difference method. A grid point formulation rather than a grid block formulation will be used since the dependent variables are specified at the boundaries rather than their flux or normal derivative. The computation domain will be discretized such that the boundary conditions and dependent variables are evaluated at $x_1,x_2, ,x_{JM1},x_{JMAX}$ and $y_1,y_2, ,y_{JM1},y_{JMAX}$ . The finite difference grid appears as follows.