This section follows the derivation of Bill Schmidt, as presented in a three day workshop in 1978 at the University of Maine.
Conservation of mass and the incompressibility condition require that
or
Conservation of linear momentum requires that
or
or
Here 
are the components of stress, 
are the
components of the body force, 
are the components of the
acceleration (zero for steady-state), 
is the density, m is the
mass, and 
are the coordinates.
Conservation of angular momentum requires that
which means that the problem is symmetric.
At this point there are more unknowns than equations (12-15 unknowns, only 7 equation). We must add the following material law or constitutive equation to generate a solvable set of equations.
or
etc., where P is the hydrostatic pressure, 
is the effective
viscosity, and 
is the kronika-delta (=1 for i=j, =0
for i not equal j).
The effective viscosity 
, depends on all the stresses. Following
Glen's Flow Law:
where 
is the effective strain rate given by the following:
or
We have one more kinematic equation that relates strain rates to gradients of velocity components:
or
where 
are the velocity components.
There is one more conservation law, conservation of energy, which would determine temperatures (B is a function of temperature). If we consider B independent of temperature, we don't need the energy equation.
At this point we apply the following assumptions or restrictions on the set of equations:
because
.
and
(the hydrostatic
equation).
To simplify formulation set up equations all in terms of velocity
components. All stresses and strain rates are then derive quantities
from the above equations. There are three nonlinear, coupled, partial
differential equations in terms of 
, 
, and P.
Applying the FEM (Galerkin Variational Technique) we construct the functional
where
Forming the variational of this functional, 
, yields the
equations encompassing the above conservation laws and constitutive
relations, as well as the boundary conditions. On the portion of the
boundary where velocity components are specified, 
, we have
On the portion of the boundary where an applied surface stress or
traction is specified, 
, we have
The domain of the problem is broken into subdomains, called elements. The global functional becomes the sum of element functionals
where for a typical element
For an element with 4 nodes (a linear quadrilateral), the x-component of the velocity at any point within the element can be expressed as a sum of values at the nodes times shape functions evaluated at the point.
where 
is the value of 
an node i in element e.
Similarly for the y-component of the velocity:
Here the N's are shape functions or interpolating polynomials, and
are defined such that each 
is equal to one at the ith node, and
equal to zero at all the other nodes in the element. With these
substitutions for 
and 
it is apparent that the functional,
J, depends only on the values of the velocity components at the nodes
and the pressures in each element.
The necessary condition for a solution is the set of u's and P that
minimize the functional, i.e. 
and
. For convenience let us define for each
node in the element a velocity vector, 
(not to be confused with
the acceleration, which is zero).
Let us similarly define a velocity vector at any point within the
element, 
.
Without subscripts this is simply
To minimize the functional we now need
for 
and
The first of these yields
From the chain rule
Recall our expression for strain rates
which we can express as a matrix product with the differential matrix
operator 
From the definition of w we have the first term of the chain rule
With the definition of 
we have for the second term
and for the third term
and hence
or finally
The other variational with respect to P gives us
where 
. For a four-node quadrilateral element the
unknowns will consist of the eight velocity components and the
pressure, (
, 
, 
, 
, 
, 
,
, 
, H), where 
to reduce round-off in the
matrix.
