# How can I make sure the flow is divergence-free when I use moving mesh?

I am using projection method and P2/P1 finite element method to solve the incompressible Navier-Stokes equations while the mesh is constantly adapted as the body moves (edge swapping, splitting and collapsing is performed locally based on a certain condition). The velocity on the new vertex is interpolated locally.

The test case considered is the flow past a cylinder at Reynolds number of 40. The drag coefficient is used to validate the results; and it compares well and does converge on the theoretical result of 1.50.

However, it can be observed from the figure that the results is not smooth and there is sudden change in drag coefficient every time when re-mesh is performed (see "drag.e4"). If we do not perform re-mesh, a smooth results are generated (see 'drag.117610' and 'drag.117740').

I think this happens because the interpolated velocity does not satisfy the divergence-free condition. I am wondering if there is any way to correct this interpolated velocity and make it divergence-free.

Thank you!

• Have you computed the divergence before and after interpolation ? Convince yourself that this is actually the problem by calculating the divergence. Most numerical methods do not actually give zero divergence, so it is not obvious that this causes your problem. – cpraveen Dec 30 '17 at 2:32
• @PraveenChandrashekar I did not compute the divergence before, but I use the projection method to solve the incompressible N-S equations and there is a step in projection method to make sure that the flow velocity satisfy the divergence free condition. – shidi.yan1992 Dec 30 '17 at 15:37

It's important to realize that the original velocity field $\mathbf v_0$ is also not divergence free in a pointwise sense. Rather, it is only divergence free when tested with the pressure test functions, i.e., $(\nabla\cdot\mathbf v_0,q)=0$ for all $q\in P_1(T_0)$ where $T_0$ is the original mesh.

What you want to achieve is to get a velocity field $\mathbf v_1$ that lives on your new mesh, somehow resembles the $\mathbf v_0$ defined on the old mesh, and that satisfies $(\nabla\cdot\mathbf v_1,q)=0$ for all $q\in P_1(T_1)$ where $T_1$ is the new mesh.

You can't do this through a local interpolation, but you can do this via a projection. However, this will require you to solve a global Laplace problem on the entire mesh -- in other words, it is no longer a local operation.

As a postscript: To find this projection, you want $\mathbf v_1$ to be closest to $\mathbf v_0$ but subject to the constraint that it is (discretely) divergence free. In other words, you want to solve the following problem: $$\min{\mathbf v_1 \in P_2(T_1)} \frac 12 \| \mathbf v_1 - \mathbf v_0 \|^2 \\ \text{subject to} \ (\nabla \cdot \mathbf v_1, q_i) = 0 \qquad \forall q_i \in P_1(T_1).$$ This is a constrained optimization problem that is easiest solved using a Lagrange multiplier approach. To this end, let $L: \in P_2(T_1) \times \in P_1(T_1) \rightarrow {\mathbb R}$ be the Lagrangian: $$L(\mathbf v_1,\lambda) = \frac 12 \| \mathbf v_1 - \mathbf v_0 \|^2 + (\nabla \cdot \mathbf v_1, \lambda).$$ Then you need to find that $\mathbf v_1, \lambda$ so that $$\frac{\partial L(\mathbf v_1,\lambda)}{\partial\mathbf v_1}(\mathbf w) = (\mathbf w, \mathbf v_1 - \mathbf v_0) + (\nabla \cdot \mathbf w, \lambda) = 0, \\ \frac{\partial L(\mathbf v_1,\lambda)}{\partial\lambda}(q) = (\nabla \cdot \mathbf v_1, q) =0 .$$ for all $w\in P_2(T_1), q \in P_1(T_1)$.

This is equivalent to finding the finite element solution of the system of equations $$\mathbf v_1 - \nabla\lambda = \mathbf v_0, \\ \nabla \cdot \mathbf v_1 = 0.$$ BY substituting the second into the (divergence of the) first equation, you will see that this corresponds to solving $-\nabla \cdot \nabla \lambda = \nabla\cdot v_0$, i.e., a Laplace equation. You could then compute $\mathbf v_1 = \mathbf v_0 +\nabla\lambda$. However, all of these equations have to be understood in the (finite dimensional) weak form shown above for the optimality condition, and you will want to solve the mixed problem instead.

• Thank you for your reply. When you say projection, do you mean that I do the integration on the new mesh i.e. $\int _{\Omega_1} \nabla^2 u_0 = \int _{\Omega_1} \nabla^2 u_1$ – shidi.yan1992 Dec 30 '17 at 16:07
• No. I added some text to the answer after the horizontal line. – Wolfgang Bangerth Dec 30 '17 at 21:12
• When you do edge swapping/splitting, you should be able to do the projection as explained by @WolfgangBangerth in a local manner. I.e. apply it only the triangles which have been modified. – cpraveen Dec 31 '17 at 5:56
• Hm, are you sure this is possible? For triangles, a local projection onto $P_2$ would only have the new edge midpoint plus the two triangle midpoints as velocity unknowns to play with (a total of 3x2 = 6) but you have have 4 pressure shape functions. Do you still get $O(h^3)$ approximation for the new velocity field in this case? – Wolfgang Bangerth Jan 1 '18 at 0:13

What sort of field is it? You've said "divergence free", but do you mean harmonic (zero divergence and curl), solenoidal (non-curl), or a mix of both? The distinction is important because it will affect the complexity of the algorithm. The goal is to do a Helmholtz decomposition numerically, interpolate the potentials, and then recalculate the velocities from the potentials. From the application, I'm going to guess that you're dealing with a velocity field that is a mix of both solenoidal and harmonic. So, you'll want to find a way to compute these integrals numerically, \begin{align} \Phi(\mathbf{r}) & = \hphantom{\frac{1}{4\pi} \int_V \frac{\nabla'\times \mathbf{v}(\mathbf{r}')}{\left|\mathbf{r}-\mathbf{r}'\right|} \operatorname{d}r'^3 }-\frac{1}{4\pi} \oint_{\partial V} \hat{n}' \cdot \frac{\mathbf{v}(\mathbf{r}')}{\left|\mathbf{r} - \mathbf{r}'\right|} \operatorname{d}^2r' \\ \mathbf{A}(\mathbf{r}) & = \frac{1}{4\pi} \int_V \frac{\nabla'\times \mathbf{v}(\mathbf{r}')}{\left|\mathbf{r}-\mathbf{r}'\right|} \operatorname{d}r'^3 - \frac{1}{4\pi} \oint_{\partial V} \hat{n}'\times \frac{\mathbf{v}(\mathbf{r}')}{\left|\mathbf{r}-\mathbf{r}'\right|}\operatorname{d}r'^2 \end{align}

Assuming that's the case, you're going to need to decide how local you want the computation to be (how large the volume $V$ is). In principle, you could make $V$ the entire space, or so small that it contains no points but has them on the surface. From there, you'd need to numerically calculate the integrals at enough points that you can numerically calculate $$\mathbf{v}(\mathbf{r}) = -\nabla \Phi(\mathbf{r}) + \nabla \times \mathbf{A}(\mathbf{r})$$ at the newly desired mesh points.

Of course, you could use these potentials as the basis for the computation, eliminating the need to step back and forth. Then you'd only need to keep the function $\Phi$ harmonic ($\nabla^2\Phi=0$) to maintain the divergence free condition. Of course, you'd also be introducing an unphysical degree of freedom in $\mathbf{A}$, and I have no idea what that would do.

If your mesh is fine enough (small compared to the size of any vortices in the calculation), you could make the local approximation that the velocity is the gradient of a harmonic scalar function. When I say local, I mean that the scalar function would only be temporarily defined for the interpolation step on a suitably small region without any expectation of global consistency (e.g. the magnetic scalar potential). Calculating new velocities based on that harmonic function will guarantee that the results remain divergence free. Deriving explicit formulae to do this is left to others, though looking information on biharmonic splines should prove useful.

• Thank you for your answer. I have not considered introducing the potentials before, I will look into it. – shidi.yan1992 Dec 30 '17 at 16:20

How about the Killing Vector Field? It is divergence-free (volume preserving), generates isometric motion but does not regularize angular motion.

Specifically, a vector field $X$ is a Killing field if the Lie derivative with respect to $X$ of the metric $g$ vanishes:

$${\mathcal {L}}_{{X}}g=0$$ It is invariant of the coordinate frame and each Killing vector corresponds to a quantity which is conserved along geodesics.