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

FEM2D_NAVIER_STOKES is a FORTRAN90 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 and stored using a form of the LINPACK/LAPACK "general band" storage. The system is factored and solved using general band factorization and solution routines similar to LINPACK's DGBFA and DGBSL 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

subroutine boundary_type ( node_num, node_xy, node_boundary, node_type, node_u_condition, node_v_condition, node_p_condition )

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

subroutine dirichlet_condition ( node_num, node_xy, u_bc, v_bc, p_bc )

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

subroutine neumann_condition ( node_num, node_xy, u_bc, v_bc, p_bc )

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

subroutine rhs ( node_num, node_xy, u_rhs, v_rhs, 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.

Program Usage:

To run the program, the user must write a file, called, perhaps, myprog.f90, containing routines defining certain data, compile this file with fem2d_navier_stokes.f90, and run the executable, supplying the node and triangle files on the command line.

gfortran fem2d_navier_stokes.f90 myprog.f90

compiles fem2d_navier_stokes.f90 and your program, and creates an executable called a.out.

a.out nodes.txt triangles.txt

runs the program with the geometry defined in nodes.txt and triangles.txt.

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 FORTRAN90 library which contains the user-supplied routines necessary to run fem2d_navier_stokes on the "cavity" problem.

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

FEM2D_NAVIER_STOKES_INOUT, a FORTRAN90 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.f90, 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 finite element 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 a case insensitive comparison of two characters for equality.
• CH_TO_DIGIT returns the integer value of a base 10 digit.
• DGB_FA performs a LINPACK-style PLU factorization of an 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.
• 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.
• GET_UNIT returns a free FORTRAN unit number.
• I4_HUGE returns a "huge" I4.
• I4COL_COMPARE compares columns I and J of an I4COL.
• I4COL_SORT_A ascending sorts an I4COL.
• I4COL_SWAP swaps columns I and J of an I4COL.
• I4MAT_DATA_READ reads data from an I4MAT file.
• I4MAT_HEADER_READ reads the header from an I4MAT.
• I4MAT_TRANSPOSE_PRINT_SOME prints some of the transpose of an I4MAT.
• JACOBIAN_ADJUST_DIRICHLET adjusts the jacobian for Dirichlet conditions.
• JACOBIAN_FEM evaluates the Navier Stokes jacobian matrix.
• LVEC_PRINT prints a logical vector.
• NEUMANN_APPLY accounts for Neumann boundary conditions.
• 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.
• R8MAT_DATA_READ reads data from an R8MAT file.
• R8MAT_HEADER_READ reads the header from an R8MAT file.
• 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.
• RESIDUAL_ADJUST_NEUMANN adjusts the residual for Neumann conditions.
• RESIDUAL_FEM evaluates the Navier Stokes residual function.
• 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.
• 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 which nodes are 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 FORTRAN90 source codes.