Navier-stokes equations on a curved manifold

General
1

Dear all,
I am trying to solve the Navier-Stokes (NS) equations for a two-dimensional fluid. The fluid does not live on a plane, but on a two-dimensional manifold (like a sphere, or ellipsoid) embedded in three-dimensional space.

The shape of the manifold is known and fixed. The ordinary derivatives appearing in the standard NS equations are replaced by covariant derivatives, which depend on the manifold shape.

I have been trying to solve this problem with the IPCS scheme, as it is done in this Fenics example for the NS equations in two dimensions. I solve the NS equations for the two components v^i of the velocity vector tangent to the manifold, as well as for the surface tension \sigma, which is a scalar on the manifold. This scheme works only when the manifold is close enough to flat. Otherwise, the method crashes and I get the error message

Starting time iteration ...
	0.78 %
	1.56 %
	2.34 %
Traceback (most recent call last):
  File "navier_stokes_membrane.py", line 280, in <module>
    solve(A1v, v_.vector(), b1v, 'bicgstab', 'hypre_amg')
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py", line 227, in solve
    return dolfin.la.solver.solve(*args)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/la/solver.py", line 72, in solve
    return cpp.la.solve(A, x, b, method, preconditioner)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to solve linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 0 iterations (PETSc reason DIVERGED_NANORINF, residual norm ||r|| = 0.000000e+00).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

after a few iterations in time. Decreasing the time step \Delta t does not help.

Are you aware of a substitute of the IPCS scheme that is adapted to my problem ? Or do you know how to pin down the source of this error ?

Best,
Michele

2

Please note that the splitting scheme in this tutorial is quite simple (and unstable). I would suggest using a more advanced scheme, like:
https://jsdokken.com/dolfinx-tutorial/chapter2/ns_code2.html
which is described in detail at
https://computationalphysiology.github.io/oasisx/splitting_schemes.html

Please also note that without a minimal reproducible example, it is really hard to give you any guidance.

I also have not seen anyone trying to use IPCS on a manifold before. Do you have any references to papers that deal with this?

3

Hello, sorry for the late reply.

Here is a minimal description of the equations. I kept this text as short as possible for clarity, please let me know if you need any additional information.

There is a two-dimensional manifold \Omega with coordinates {\bf x} = (x^1, x^2), and metric tensor g_{ij}({\bf x}). Specifically, x^1 = x, x^2 = y and z are Cartesian coordinates in a plane and \Omega is a surface in three dimensions given by the height function z= z(x,y)

Consider \Omega as a layer on which a fluid flows, with velocity tangential to the manifold \bf v ( a vector field in the tangent bundle of \Omega) and velocity normal to \Omega w (a scalar on \Omega).

The equations for the velocities read

\rho (\partial_t v^i + v^j \nabla_j v^i - 2 w v^j b^i_j -w \nabla^i w ) = \nabla ^i \sigma + \eta[- \nabla_{\rm LB} v^i + \cdots ],
\rho [\partial_t w + v^i (v^j b_{ji} +\nabla_i w )] = \kappa [2 \nabla_{\rm LB} H + \cdots ] + 2 \sigma H + 2 \eta [(\nabla^i v_j) b_{ij} + \cdots],
and the continuity equation reads
\nabla_i v^i - 2 H w =0.

Here

  • \sigma is the surface tension of the fluid, \eta the viscosity, \rho the density

  • indexes are raised an lowered according to Einstein notation with g

  • \nabla is the covariant derivative and b_{ij} the first funcamental form on \Omega, they depend on z({\bf x})

  • \nabla_{\rm LB } is the Laplace-Beltrami operator on \Omega

Finally, the manifold profile z({\bf x}) evolves according to the velocities \bf v and w with
\partial_t z = (v^i {\bf e}_i + w \hat{n})\cdot \hat{z} + \cdots ,

where {\bf e}_i = \partial / \partial x_i are the vectors tangent to the coordinate lines on \Omega and \hat n is the unit normal to \Omega.

I am trying to solve the three coupled equations above, for the fields {\bf v}, w and z.

I managed to succcesfully solve them if I start from a nearly flat manifold, with both the IPCS scheme and this more stable scheme. However, as soon as I iterate in time for too long or start with a strongly non-flat manifold, the code crashes: I get this kind of error

Starting time iteration ...
	0.78 %
	1.56 %
Traceback (most recent call last):
  File "code.py", line 267, in <module>
    solve(A1, vw_.vector(), b1, 'bicgstab', 'hypre_amg')
  File "/usr/local/lib/python3.6/dist-packages/dolfin/fem/solving.py", line 227, in solve
    return dolfin.la.solver.solve(*args)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/la/solver.py", line 72, in solve
    return cpp.la.solve(A, x, b, method, preconditioner)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to solve linear system using PETSc Krylov solver.
*** Reason:  Solution failed to converge in 0 iterations (PETSc reason DIVERGED_NANORINF, residual norm ||r|| = 0.000000e+00).
*** Where:   This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f

Do you know a scheme that would be more adapted to this kind of problem?
Thank you !