manufactured solutions for incompressible Navier-Stokes -- how to find divergence-free velocity fields?

Question

In the method of manufactured solutions (MMS) one postulates an exact solution, substitutes it in the equations and calculates the corresponding source term. The solution is then used for code verification.

For incompressible Navier-Stokes equations, MMS easily leads to a (non-zero) source term in the continuity equation. But not all codes allow source terms in the continuity equations, so for these codes only manufactured solutions with a divergence-free velocity fields will do. I found this example for a domain $\Omega=[0,1]^2$ \begin{align} u_1 &= -\cos(\pi x) \sin(\pi y) \\ u_2 &= \sin(\pi x) \cos(\pi y) \end{align} In general 3D cases, how does one manufacture a divergence-free velocity field?

Answer 1

Use a vector streamfunction or take the cross product of two gradients. I.e.: $$ \boldsymbol{u}=\nabla\times\boldsymbol{A} $$ where $\boldsymbol{A}$ is a vector field of your choosing, or $$ \boldsymbol{u}=\nabla f\times\nabla g $$ where $f$ and $g$ are two scalar fields of your choosing.

It's hard both have the velocity be divergence-free and prescribe the boundary conditions, but as long as your code allows you to set arbitrary functions for your boundary conditions, you should be OK.

ETA: Of course, your momentum equation will have to accept a forcing function, but I've always felt better about forcing the momentum equation than I have adding a right-hand side to the continuity equation.

Answer 2

This is not a general answer, but for the Navier-Stokes equations, there are manufactured solutions that describe real flow. For example, the Kovasznay flow field is a popular choice:

http://link.springer.com/article/10.1007/BF00948290

The original reference is: Kovasznay L.I.G., "Laminar flow behind a two-dimensional grid". Proc. Cambridge Philos. Soc., page 44, 1948.

Answer 3

That's what I usually do.

Define streamline function:

$$ \Psi = \left[ \begin{array}{c} \psi_x\\ \psi_y\\ \psi_z\\ \end{array} \right] $$

the velocity equals:

$$ \mathbf{u} = \nabla\times\Psi = \left[ \begin{array}{c} u_x =\partial_y \psi_z - \partial_z \psi_y\\ u_y =\partial_z \psi_x - \partial_x \psi_z\\ u_z =\partial_x \psi_y - \partial_y \psi_x\\ \end{array} \right]. $$

Now you can pick any reasonable zero-averaged pressure and construct a forcing term.

I post a SymPy example code for $\Omega = [0,1]^3$ and homogeneous boundary conditions, enjoy:

 from sympy import *

 x,y,z = symbols('x y z')

 X = Matrix([[x],[y],[z]])

 psi = zeros(3,1)
 psi[0,0] = sin(2*pi*x)*y**2*(1-y)**2*z**2*(1-z)**2
 psi[2,0] = x**2*(1-x)**2*y**2*(1-y)**2*sin(2*pi*z)

 curl_psi = zeros(3,1)
 curl_psi[0] = diff(psi[2],X[1]) - diff(psi[1],X[2])
 curl_psi[1] = diff(psi[0],X[2]) - diff(psi[2],X[0])
 curl_psi[2] = diff(psi[1],X[0]) - diff(psi[0],X[1])