FEM2D_NAVIER_STOKES
Steady Incompressible Navier Stokes Equations in 2D
Finite Element Solution
Banded Storage

FEM2D_NAVIER_STOKES is a C++ program which applies the finite element method to solve, over an arbitrary triangulated region, a discretized version of the nonlinear partial differential equations for steady incompressible Navier Stokes flow.

The geometry is entirely external to the program. The user specifies one file of nodal coordinates, and one file that describes the triangles in terms of six node coordinates.

The program makes a default assumption that all boundary velocities correspond to Dirichlet boundary conditions, and that one pressure is specified (for uniqueness of the pressure system). The user can adjust these boundary conditions (and even specify Dirichlet constraints on any variable at any node) by setting the appropriate data in certain user routines.

At the moment, Neumann conditions, if specified, must have a zero right hand side. The machinery to integrate a nonzero Neumann condition has not been set up yet.

The linear system is created, stored, and solved using a form of the LINPACK/LAPACK "general band" storage and versions of LINPACK's DGBFA and DGBSL factorization and solution routines.

Computational Region

The computational region is initially unknown by the program. The user specifies it by preparing a file containing the coordinates of the nodes, and a file containing the indices of nodes that make up triangles that form a triangulation of the region. For the following ridiculously small example:

       10-11-12
        |\    |\
        | \   | \
        6  7  8  9
        |   \ |   \
        |    \|    \
        1--2--3--4--5
      

         0.0 0.0
         1.0 0.0
         2.0 0.0
         3.0 0.0
         0.0 1.0
         1.0 1.0
         2.0 1.0
         3.0 1.0
         0.0 2.0
         1.0 2.0
         2.0 2.0
      

         1  3  10  2  7  6
        12 10  3  11  7  8
         3  5 12   4  9  8
      

The Navier Stokes Equations

The state variables are a velocity vector (U,V)(X,Y) and a scalar pressure P(X,Y). The state variables are constrained by the momentum and continuity equations, which apply inside the region:

        - nu * (Uxx + Uxx) + U * dU/dx + V * dU/dy + dP/dx  = u_rhs(x,y)

        - nu * (Vxx + Vyy) + U * dV/dx + V * dV/dy + dP/dy  = v_rhs(x,y)

        dU/dx + dV/dy = p_rhs(x,y)
      

Boundary Conditions

At every point on the boundary of the region, the program assumes that both components of the velocity are specified.

        U(node) = U_BC(node)
        V(node) = V_BC(node)
      

At one point in the region, the program assumes that the value of the pressure is specified.

        P(node) = P_BC(node)
      

The user routine boundary_type can be used to modify the type of the boundary conditions associated with a degree of freedom at a boundary node - or even to add constraints to variables associated with nodes in the interior.

One choice that the user can make is to reset certain boundary conditions to be of Neumann type:

        dU/dn(node) = U_BC(node)
        dV/dn(node) = V_BC(node)
      

Computational Procedure

We use linear finite elements for the pressure function, and to generate these, we only need the nodes that are the vertices of the triangles. (We will call these vertices "pressure nodes.") Because quadratic basis functions are to be used for the velocity, however, each triangle will also have three extra midside nodes available for that.

We now assume that the unknown velocity component functions U(x,y) and V(x,y) can be represented as linear combinations of the quadratic basis functions associated with each node, and that the scalar pressure P(x,y) can be represented as a linear combination of the linear basis functions associated with each pressure node.

For every node, we can determine a quadratic velocity basis function PSI(I)(x,y). For every pressure node I, we can determine a linear basis function PHI(I)(x,y). If we assume that our solutions are linear combinations of these basis functions, then we need to solve for the coefficients.

To do so, we apply the Galerkin-Petrov method. For each pressure node, and its corresponding basis function PHI(I), we generate a copy of the continuity equation, multiplied by that basis function, and integrated over the region:

Integral ( dUdx(x,y) + dVdy(x,y) ) * PHI(I)(x,y) dx dy = Integral ( P_RHS(x,y) * PHI(I)(x,y) dx dy )

Similarly, for each node and its corresponding velocity basis function PSI(I), we generate two copies of the momentum equation, one for each component. We multiply by PSI(I), integrate over the region, and use integration by parts to lower the order of differentiation. This gives us:

Integral ( nu * ( dPSIdx(I) * dUdx + dPSIdy(I) * dUdy ) + PSI(I) * ( ( U * dUdx + V * dUdy ) + dPdx - U_RHS ) dx dy ) = 0
Integral ( nu * ( dPSIdx(I) * dVdx + dPSIdy(I) * dVdy ) + PSI(I) * ( ( U * dVdx + V * dVdy ) + dPdy - V_RHS ) dx dy ) = 0

After adjusting for the boundary conditions, the set of all such equations yields a nonlinear system for the coefficients of the finite element representation of the solution. To solve such a nonlinear system, we apply Newton's method. This means we must have a starting solution. For starting solution, we solve the (linear) Stokes problem first. The value of this initial approximate solution is printed out, since the Stokes solution is of interest in its own right. Then the Newton iteration is applied to determine the solution of the nonlinear system.

User Input Routines

The program requires the user to supply the following routines:

The default boundary condition types are passed to the user, along with other information. The user modifies any data as necessary, and returns it. This is done by a user-supplied routine of the form

void boundary_type ( int node_num, double node_xy[], bool node_boundary[], int node_type[], int node_u_condition[], int node_v_condition[], int node_p_condition[] )

The right hand sides of any Dirichlet boundary conditions are determined by a user-supplied routine of the form

void dirichlet_condition ( int node_num, double node_xy[], double u_bc[], double v_bc[], double p_bc[] )

The right hand sides of any Neumann boundary conditions are determined by a user-supplied routine of the form

void neumann_condition ( int node_num, double node_xy[], double u_bc[], double v_bc[], double p_bc[] )

The right hand sides (or "source terms") of the state equations are determined by a user-supplied routine of the form

void rhs ( int node_num, double node_xy[], double u_rhs[], double v_rhs[], double p_rhs[] )

Program Output

The program writes out various node, triangle, pressure and velocity data files that can be used to create plots of the geometry and the solution.

Graphics files created include:

• nodes6.eps, an image of the nodes;
• triangles6.eps, an image of the elements;

Data files created include:

• nodes3.txt, the nodes associated with pressure;
• triangles3.txt, the linear triangles associated with pressure;
• stokes_pressure3.txt, the pressure at the pressure nodes for the initial Stokes solution (if requested);
• stokes_velocity6.txt, the velocity at the velocity nodes for the initial Stokes solution (if requested).
• pressure3.txt, the pressure at the pressure nodes for the final Navier Stokes solution;
• velocity6.txt, the velocity at the velocity nodes for the final Navier Stokes solution.

Usage:

To run the program, the user must first write a file, called, perhaps, myprog.C, containing routines defining certain data, compile this file with fem2d_navier_stokes.C, and create an executable. It would be sensible to rename the executable to fem2d_navier_stokes.

g++ fem2d_navier_stokes.C myprog.C
mv a.out fem2d_navier_stokes

Once the executable is created, it can be run as follows:

fem2d_navier_stokes nodes.txt triangles.txt

• nodes.txt contains the coordinates of nodes;
• triangles.txt lists the indices of nodes that form elements.

Licensing:

The computer code and data files described and made available on this web page are distributed under the GNU LGPL license.

Languages:

FEM2D_NAVIER_STOKES is available in a C++ version and a FORTRAN90 version and a MATLAB version.

Related Data and Programs:

FEM2D_NAVIER_STOKES_CAVITY, a C++ library which contains the user-supplied routines necessary to run fem2d_navier_stokes on the "cavity" problem.

FEM2D_NAVIER_STOKES_CHANNEL, a C++ library which contains the user-supplied routines necessary to run fem2d_navier_stokes on the "channel" problem.

FEM2D_NAVIER_STOKES_INOUT, a C++ library which contains the user-supplied routines necessary to run fem2d_navier_stokes on the "inout" problem.

Reference:

1. Max Gunzburger,
   Finite Element Methods for Viscous Incompressible Flows,
   A Guide to Theory, Practice, and Algorithms,
   Academic Press, 1989,
   ISBN: 0-12-307350-2,
   LC: TA357.G86.
2. Hans Rudolf Schwarz,
   Finite Element Methods,
   Academic Press, 1988,
   ISBN: 0126330107,
   LC: TA347.F5.S3313.
3. Gilbert Strang, George Fix,
   An Analysis of the Finite Element Method,
   Cambridge, 1973,
   ISBN: 096140888X,
   LC: TA335.S77.
4. Olgierd Zienkiewicz,
   The Finite Element Method,
   Sixth Edition,
   Butterworth-Heinemann, 2005,
   ISBN: 0750663200,
   LC: TA640.2.Z54

Source Code:

• fem2d_navier_stokes.cpp, the source code;
• fem2d_navier_stokes.sh, commands to compile the partial program;

List of Routines:

• MAIN is the main routine of FEM2D_NAVIER_STOKES.
• ASSEMBLE_STOKES assembles the system for the Stokes equations.
• BANDWIDTH determines the bandwidth of the coefficient matrix.
• BASIS_MN_T3: all bases at N points for a T3 element.
• BASIS_MN_T6: all bases at N points for a T6 element.
• CH_CAP capitalizes a single character.
• CH_EQI is true if two characters are equal, disregarding case.
• CH_TO_DIGIT returns the integer value of a base 10 digit.
• DGB_FA performs a LINPACK-style PLU factorization of a DGB matrix.
• DGB_PRINT_SOME prints some of a DGB matrix.
• DGB_SL solves a system factored by DGB_FA.
• DIRICHLET_APPLY accounts for Dirichlet boundary conditions.
• DTABLE_DATA_READ reads the data from a real TABLE file.
• DTABLE_HEADER_READ reads the header from a real TABLE file.
• FILE_COLUMN_COUNT counts the number of columns in the first line of a file.
• FILE_NAME_SPECIFICATION determines the names of the input files.
• FILE_ROW_COUNT counts the number of row records in a file.
• I4_HUGE returns a "huge" I4.
• I4_MAX returns the maximum of two integers.
• I4_MIN returns the smaller of two integers.
• I4COL_COMPARE compares columns I and J of an I4COL.
• I4COL_SORT_A ascending sorts the columns of an I4COl.
• I4COL_SWAP swaps two columns of an I4COL.
• I4MAT_TRANSPOSE_PRINT_SOME prints some of an I4MAT, transposed.
• ITABLE_DATA_READ reads data from an integer TABLE file.
• ITABLE_HEADER_READ reads the header from an integer TABLE file.
• JACOBIAN_ADJUST_DIRICHLET adjusts the jacobian for Dirichlet conditions.
• JACOBIAN_FEM evaluates the Navier Stokes jacobian matrix.
• LVEC_PRINT prints a logical vector.
• NODES3_WRITE writes the pressure nodes to a file.
• POINTS_PLOT plots a pointset.
• PRESSURE3_WRITE writes the pressures to a file.
• QUAD_RULE sets the quadrature rule for assembly.
• R8_HUGE returns a "huge" R8.
• R8_NINT returns the nearest integer to an R8.
• R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed.
• R8VEC_PRINT_SOME prints "some" of an R8VEC.
• REFERENCE_TO_PHYSICAL_T6 maps T6 reference points to physical points.
• RESIDUAL_ADJUST_DIRICHLET adjusts the residual for Dirichlet conditions.
• S_LEN_TRIM returns the length of a string to the last nonblank.
• S_TO_I4 reads an I4 from a string.
• S_TO_I4VEC reads an integer vector from a string.
• S_TO_R8 reads an R8 from a string.
• S_TO_R8VEC reads an R8VEC from a string.
• S_WORD_COUNT counts the number of "words" in a string.
• SORT_HEAP_EXTERNAL externally sorts a list of items into ascending order.
• TIMESTAMP prints the current YMDHMS date as a time stamp.
• TIMESTRING returns the current YMDHMS date as a string.
• TRIANGLE_AREA_2D computes the area of a triangle in 2D.
• TRIANGLES3_WRITE writes the pressure triangles to a file.
• TRIANGULATION_ORDER6_BOUNDARY_NODE indicates nodes on the boundary.
• TRIANGULATION_ORDER6_PLOT plots a 6-node triangulation of a set of nodes.
• VELOCITY6_WRITE writes the velocities to a file.

You can go up one level to the C++ source codes.