3D Divergence-free Stokes equations

I would like to have a divergence-free formulation of the Navier-Stokes equations for creeping flow (a.k.a. Stokes equation, N.-S. without the inertial terms) in 3D, for purposes of tracking particle deformation in finite elements.

I would also like to be able to specify the usual kinds of boundary conditions (no-slip conditions, moving walls, slip conditions, pressure).

For calculating divergence-free flows, I have found a few options.

One is $$\vec{u} = \operatorname{grad} \psi_1 \times \operatorname{grad} \psi_2$$ The problem with this is that $\psi_1$ and $\psi_2$ are not uniquely determined, and that imposing boundary conditions for non-zero velocities is not trivial.

Another possibility I have found is to take the curl of a vector, $$\vec{u} = \operatorname{curl} \left( \begin{matrix} \psi_1\\ \psi_2\\ \psi_3 \end{matrix} \right)$$

where additional conditions (what kind?) have to be imposed to ensure uniqueness of the solution. Again, specifying boundary conditions is non-trivial.

Any advice or pointers to literature on how to proceed?

A clarification: I want to solve the Stokes equation. For my application, it is important that the divergence of the flow field vanishes identically (subject to small errors), also inside an element. If you look at the results of the standard solutions that Comsol delivers with its CFD code, you will find quite large errors for $\operatorname{div} \vec{u}$ for flow regions close to a stagnation point.

• Judging from the answers, and judging from what I have read from the literature, it seems that what I want is a more or less active research topic, so it does not appear to be doable in Comsol without a lot of effort. – tkoencov Aug 19 '15 at 15:41