FEM2D_STOKES_SPARSE is a MATLAB program which applies the finite element method to solve a form of the steady incompressible Stokes's equations over an arbitrary triangulated region in 2D.
The linear system is created and stored using MATLAB's sparse matrix storage. The system is factored and solved by MATLAB, using sparse matrix solution techniques.
This program is derived from a similar program, FREE_FEM_STOKES, which uses banded storage, factorization and solution methods.
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 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
the node file could be:
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
and the triangle file would be
1 3 10 2 7 6
12 10 3 11 7 8
3 5 12 4 9 8
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 + Uyy ) + dP/dx = U_RHS(x,y)
- nu * ( Vxx + Vyy ) + dP/dy = V_RHS(x,y)
dU/dx + dV/dy = P_RHS(x,y)
where, typically, the right hand side functions are zero.
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)
This is known as a "Dirichlet boundary condition". The right hand side
of the boundary condition is left unspecified until the user
routine is called. If a wall is intended, then the appropriate
condition has U_BC and V_BC zero. An inlet might have a
particular flow profile function used to assign nonzero values.
At one point in the region, the program assumes that the value of the pressure is specified.
P(node) = P_BC(node)
Such a condition must be supplied; otherwise the pressure cannot
be uniquely determined, since it is essentially a potential function,
and so is unique only "up to a constant". Note that the program allows
the user to specify pressure conditions anyway, and these can be
of Dirichlet or Neumann type. In general, however, this is not
a physically or mathematically or computationally good thing to do!
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)
The right hand side of the boundary condition is left unspecified
until the user routine is called. As mentioned earlier, the program
cannot currently handle Neumann conditions with nonzero right hand side.
(A nonzero value is simply ignored, but won't actually cause the
program to fail.)
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 ( Ux(x,y) + Vy(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 * ( Ux(x,y) * PSIx(I)(x,y) + Uy(x,y) * PSIy(I)(x,y) ) + Px(x,y) * PSI(I)(x,y) dx dy = Integral ( U_RHS(x,y) * PSI(I)(x,y) dx dy )
Integral nu * ( Vx(x,y) * PSIx(I)(x,y) + Vy(x,y) * PSIy(I)(x,y) ) + Py(x,y) * PSI(I)(x,y) dx dy = Integral ( V_RHS(x,y) * PSI(I)(x,y) dx dy )
After adjusting for the boundary conditions, the set of all such equations yields a linear system for the coefficients of the finite element representation of the solution.
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 MATLAB M file of the form:
function [ node_u_condition, node_v_condition, node_p_condition ] = boundary_type ( node_num, node_xy, node_boundary, node_type, node_u_condition, node_v_condition, node_p_condition )
The value of the kinematic viscosity is determined by a user-supplied MATLAB M file of the form
function nu = constants ( )
The right hand side of any Dirichlet boundary conditions are determined by a user-supplied MATLAB M file of the form
function [ u_bc, v_bc, p_bc ] = dirichlet_condition ( n, xy )
The right hand side of any Neumann boundary conditions are determined by a user-supplied MATLAB M file of the form
function [ u_bc, v_bc, p_bc ] = neumann_condition ( n, xy )
The right hand side source term functions are determined by a user-supplied MATLAB M file of the form
function [ u_rhs, v_rhs, p_rhs ] = rhs ( m, xy )
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:
Data files created include:
fem2d_stokes_sparse ( 'prefix' )where 'prefix' is the common input filename prefix:
The computer code and data files described and made available on this web page are distributed under the GNU LGPL license.
FEM2D_STOKES_SPARSE_BIG_CAVITY, a MATLAB library which contains the user-supplied routines necessary to run fem2d_stokes_sparse on the "big cavity" problem.
FEM2D_STOKES_SPARSE_CAVITY, a MATLAB library which contains the user-supplied routines necessary to run fem2d_stokes_sparse on the "cavity" problem.
FEM2D_STOKES_SPARSE_CHANNEL, a MATLAB library which contains the user-supplied routines necessary to run fem2d_stokes_sparse on the "channel" problem.
FEM2D_STOKES_SPARSE_INOUT, a MATLAB library which contains the user-supplied routines necessary to run fem2d_stokes_sparse on the "inout" problem.
FEM2D_STOKES_SPARSE_SPIRAL, a MATLAB library which contains the user-supplied routines necessary to run fem2d_stokes_sparse on the "spiral" problem.
You can go up one level to the MATLAB source codes.