From our element equation:
we apply a Penalty Method whereby the pressure, P is assumed to depend on the velocity gradients in the following form
where 
is a sufficiently large positive number. The element
equation becomes:
This yields a matrix equation
where 
is given by
and 
is
given by
For a stable solution, it turns out to be necessary that 
be singular []. This can be accomplished by using
reduced integration in the numerical integration []
[]. This simply consists of using a lower-order Gaussian
integration rule for evaluation of the integrals in 
while using a sufficiently high-order rule for the integration of
. This requirement actually reduces the computational load in
generating the matrix equations.
In addition we have also gained in the elimination of the additional unknown representing the pressure in each element. As an added benefit, the matrix is more closely banded (without the kX1 and 1Xk pressure terms, as well as the elimination of the zeroes along the diagonal of the pressure equations which effectively eliminated the possibility of an iterative solution of the matrix equations.
The penalty constant, 
, is typically chosen to be
where c is a constant on the order of 
[]. This
choice seems to yield reasonable results with regards compressibility
and pressures, without producing numerical ill conditioning.