Negative viscosity stabalized by fourth order terms

Question

I am trying to solve a "Navier-Stokes"-type problem where the viscosity is negative. Of course this renders the equation unstable and thus I add a fourth order term, so the entire equation becomes:

dv/dt + (v*nap)v = - nap p + nu * n^2 v - n2 * nap^4 v

where n < 0 and n2 > 0 and nap = del operator.

Now the fourth order term should stabalize the solution with negative viscosity - I can show this by stability analysis - but I can't get fenics to solve it. The solutions just blow up.

I followed this guide: http://fenicsproject.org/documentation/dolfin/dev/python/demo/pde/navier-stokes/python/documentation.html to get fenics to solve the basic Navier Stokes. I then added the term

n2*inner( div(grad(u0)), div(grad(v)))*dx

and I also tried with

n2*inner( div(grad(div(grad(u0)))),v)*dx

Neither works - and surprisingly they also give different results. As you can probably infer from my questions, I am new to this CFM, but would really like some help.

Thanks in advance.

Answer 1

Interesting question and I only can give some hints.

Because of the pressure-velocity interaction, the time integration of the Navier-Stokes Equations has its own stability issues. So I suggest you first try your stabilization for a convection-diffusion equation for a scalar quantity $\rho$.

Then you should give some more thoughts on what is the right weak formulation of the fourth order term and what are the right finite-elements to use. E.g., if your ansatz functions are not smooth enough, you need to take a discontinuous Galerkin framework.

Regarding the weak formulation: The commonly used identity $(\nabla \cdot v, q ) = -(v,\nabla q)$, where $v$...vector and $q$...scalar, is only valid under particular conditions. In your case, the '2nd' application is not justified, because $\nabla \cdot (\nabla \rho)$ will not be zero at the boundary. That's also the reason, why you get different results depending on your formulation.

Furthermore, $\nabla v$ is actually a tensor and the relation of $\nabla \cdot (\nabla v) = \Delta v$ includes some assumptions (like $\nabla \cdot v=0$) that may not hold in a repeated application.

All in all, I do think it is a question of numerical analysis rather than of fenics implementation. If you have the equations ready, you may try to get help at http://scicomp.stackexchange.com/

Answer 2

The vanilla finite element method will not provide the regularity that you need for this problem.

For some ideas on how to solve problems with fourth-order spatial derivatives, for scalar equations take a look at the demo/documented/cahn-hilliard for a operator splitting approach and at demo/documented/biharmonic for a $C^{0}$ discontinuous Galerkin (DG) method. You can find a $C^{0}$ discontinuous Galerkin for fourth-order, nonlinear elasticity equations in the supporting material to the paper http://dx.doi.org/10.1016/j.jmps.2011.04, which is at https://www.repository.cam.ac.uk/handle/1810/236975. The model in the paper can have negative stiffness, which is very similar to what you're looking at.