Next: Appendix F: Ellipse
Previous: Appendix D: Exact
Finite element methods  are well suited to find the extrema of 'energy functionals'. In the physical cases considered in this paper, where variational principles exist from which the relevant equations can be derived, the energy functional to be minimized is the actual energy of the system, where perhaps Lagrange multiplier terms have been added to enforce constraints. The appropriate functional for elastic problems is the elastic energy ,
where εij is the strain tensor and Cijkl is the elastic stiffness tensor. The appropriate functional for the conductivity problem  is the dissipated electrical power,
where Ei is the electric field vector and σij is the conductivity tensor.
The elastic-fluid mapping mentioned previously [17,102,103,147] implies that a similar variational principle exists for the linearized fluid problem (Stokes equation), where only terms linear in the velocities are retained. This fact has been known for some time, being first discussed by Helmholz  and proven and elaborated on by Korteweg . Millikan , however, proved that if one is restricted to action integrals involving only fluid velocities and first order spatial derivatives, then there is no variational principle whose Euler equations give the full Navier-Stokes equations.
Two recent discussions of the variational principle for viscous fluids were given by de Veubeke  and Keller et. al. , where the pressure is introduced as a Lagrange multiplier :
where F is the rate of energy dissipation, u is the fluid velocity, and p is the pressure. Carrying out the variation indicated in eqn. (105) gives the linearized steady-state Navier-Stokes equation in the creeping flow limit, commonly called the Stokes equation. To enforce the incompressibility condition · u = 0 requires an extra effort on the part of the would-be solver of eqn. (105).
In the finite element solution used in this paper we have used a formulation of the problem described by Zienkiewicz  in which the pressure is ignored and the incompressibility condition is only approximately maintained via a "penalty method" . A term is added to the energy dissipation F of the form (1/2) β ( · u)2 where taking the 'penalty parameter' large (β → ∞) corresponds to the incompressible limit ( · u = 0). This method works extremely well, although run times become longer as the value of β is increased. Finite values of β imply some degree of compressibility of the simulated fluid.
However, it turns out that the intrinsic viscosity, or more specifically, the extra dissipated energy terms in eqn. (89) from Appendix C do not depend, (under steady-state conditions) on the compressibility of the fluid, which means that the intrinsic viscosity does not depend on the compressibility either . This allowed run times to be considerably shortened, as much snaller values of β could be used with no change in the final results. The argument goes as follows. The stress tensor of a compressible isotropic fluid is a function of both the shear viscosity η and the bulk viscosity ξ . The total stress tensor σij is the sum of an ideal fluid pressure component and a viscosity component,
where p is the macroscopic pressure and the viscosity contribution to the stress tensor, , reflects the explicit dependence on η and ξ ,
For an incompressible fluid, the viscous stress tensor evidently reduces to θij in eqn. (106). The compressibility-related terms in Eqn. (106) can be absorbed into a redefinition of the pressure Pcomp which then varies with the local spatial position,
We then observe that the integral of the gradient of Pcomp over any closed circuit C in the fluid equals zero,
since Pcomp is a single-valued function of position. Consequently, there is no extra dissipation  associated with the fluid compressibility ( · u 0) . Notably eqn. (111) is independent of the surface boundary condition.
Finally, in actually carrying out the finite element computation for the viscosity, the "no-slip" boundary condition at a fluid-solid boundary is enforced by holding the velocity at all computational nodes either at the boundary or inside the particle fixed at zero during the conjugate gradient  cycles that find the minimum energy. To treat a plate geometry one simply holds the velocities of all the nodes in the appropriate planar region fixed at zero. In the coresponding electrical problem one simply holds the voltage constant at the same nodes as in the fluid problem.