Redefining velocity field to make it divergence free on the grid level

Question

I am currently working with convection equation $\frac{\partial C}{\partial t}+u\frac{\partial C}{\partial x}+v\frac{\partial C}{\partial y}=0$. I assume that the velocity field is incompressible i.e. $\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$. However in most scenarios, I don't have the explicit expression for $u$ and $v$ at all points. But instead, I only have the data about the velocities at certain locations only.

More precisely, suppose I have a velocity field $(u,v)$ defined only at the grid points $(i\Delta x,j\Delta y)$ with $i=0,1,\cdots,N$ and $j=0,1,\cdots,N$ (let's just say that the domain is a square). Is there a some sort of algorithm to "redefine" the velocity field so that it is close to divergence free at the grid level (in other words $\frac{u_{i+1,j}-u_{i,j}}{\Delta x}+\frac{v_{i,j+1}-v_{i,j}}{\Delta y}$ is very close to zero) and the redefined vector field is "close" to the original one? Would be helpful if there is a paper or reference that I can look at on this matter.

Answer 1

In finite element simulations, we have similar issues for simulations where the spatial mesh changes, e.g. for adaptive meshes that change between time steps. There you have that the solution from one mesh is no longer discretely divergence-free. In the finite element context, we also have basically the same two solutions as in the answer by NNN.

Firstly, you could solve an auxillary Stokes problem. The math behind this is desribed in the paper by Besier and Wollner. Two possible ways are the $L^2$-projection from equation (4.6) and the $H^1$-projection from equation (4.12). In the paper it is all written with weak formulations, but basically you solve for the $L^2$-projection the PDE: $$ v + \nabla p = v_{data} \\ \nabla \cdot v = 0 $$ or for the $H^1$-projection the Stokes problem: $$ -\Delta v + \nabla p = -\Delta v_{data} \\ \nabla \cdot v = 0 $$ Here, I denoted by $v_{data}$ your velocity field that might not yet be divergence free. This should be related to the last optimization problem from NNN.

Secondly, you could also use special finite elements that are divergence-free and for this you could use the Helmholtz decomposition. A lot of work on this issue has been done by Alexander Linke et al. You can find this e.g. in Linke's paper from 2014 or one of his newer works. The only limitation here is that this is mostly being done for triangular elements and for quadrilateral elements there is much less literature.

Answer 2

I'm surprised that this question hasn't received an authoritative answer. I'd like to know an authoritative answer to this question myself.

This is my stab at it, which I think works in 3D. In 2D the curl would be along the $k$ direction and your velocity in $i,j$, which probably wouldn't lead to solutions that you want.

From wikipedia we know that

$$ \nabla\cdot(\nabla\times \mathbf{A}) = 0 $$

So, how about minimizing wrt $\mathbf{A}$:

$$ \pi(\mathbf{A}) = \frac{1}{2}\|\mathbf{u}-(\nabla\times \mathbf{A} )\|^2 $$

Then $\nabla\times \mathbf{A}$ is your divergence free field. Note: I've omitted the spatial dependence of $\mathbf{A}$ and $\mathbf{u}$ (your original field that you want to make divergence free).

Leray Projection also comes to mind. But apart from the Wikipedia link, I don't have the ability to tell you anything more about it.

You might also want to try the Helmholtz-Hodge decomposition, though the integrals seem complicated. See also the section on "Fields with prescribed divergence and curl".

Finally, you might simply want to minimize

$$ \pi(\mathbf{v}) = \frac{1}{2}\|\mathbf{u}-\mathbf{v}\|^2 + \lambda\|\nabla\cdot \mathbf{v}\|^2 $$

Where $\lambda$ is a Lagrange multiplier or a large penalty constant. $\mathbf{v}$ would be your divergence free field.