Open boundary conditions for Navier Stokes

1

I’m trying to solve Navier Stokes equations in an ALE formulation with the Oasis solver (GitHub - mikaem/Oasis), with FEniCS 2019.2. The problem is in 3D and I’m using a pressure correction scheme, where a splitting of pressure and velocity (Chorin type) is used for reasons of computational costs, compared to a monolithic solver.

I am observing blow ups at the exit corners of the domain, where anomalous ricirculation phenomena occur. therefore I would like to apply open boundary conditions on the outer boundary of the domain, in order to allow vortices to leave the domain, instead of the classic zero Neumann for velocity and zero Dirichlet for pressure.

I don’t know of any implementations of open boundary conditions in the finite element weak formulation, only for finite volumes or spectral elements. Furthermore, some solutions seem very involved and difficult to implement in FEniCS. A promising method would be, for example, the one in the article from Eberhard Bänsch:

A finite element pressure correction scheme for the Navier–Stokes equations with traction boundary condition, CMAME, 2014

These conditions would be suitable because they don’t couple the three equations for the velocity components. However, I don’t know how to implement these in FEniCS. In particular, to solve for the pressure first I have to solve an equation on the domain boundary \Gamma of this kind:

(\rho_h, \text{tr}_{\Gamma} q_h)_{\Gamma} +\epsilon(\nabla_{\Gamma} \rho_h, \nabla_{\Gamma} \text{tr}_{\Gamma} q_h)_{\Gamma} = \int_{\Omega} \nabla \Phi\cdot \nabla q_h + \frac{1}{\Delta t} \int_{\Omega} \nabla \cdot \mathbf{u} q_h

And, once \rho_h has been solved for, I need to calculate:

\langle l_h, v_h \rangle := \int_{\Gamma}\rho_h \nabla_{\Gamma}\cdot (\text{tr}_{\Gamma} v_h)

where a tangential divergence \nabla_{\Gamma} has to be calculated. I don’t know how to implement this divergence form in FEniCS, it seems a bit involved, because my problem in in 3D (suppose the boundary is the surface of a sphere).

Does anyone know how to implement/had some success with Open boundary BCs for Navier Stokes in FEniCS? Do you know of any examples?

2

I would try the outflow boundary condition given in Section 2.2.1 of this paper, which I have used many times myself. I think the easiest way to understand the motivation behind this boundary condition is to first take a look at the formulation of stable Neumann BCs for advection–diffusion, given in this report.

2 Likes
3

Thank you @kamensky for the reply. I am reading the report, it seems that this techniques work better for stabilized finite elements, since the inf-sup condition makes it harder to respect conservation.
Would this still be beneficial for the classic P2/P1 elements?

4

The stability of Neumann BCs on inflow portions of the boundary is a well-posedness issue of the continuous problem. It is independent of the numerical discretization used. While one might walk away from the first paper I linked with the impression that the extra boundary term proposed in Section 2.2.1 is a stabilization technique for use in numerical simulations, you should really think of it as a feature of the continuous problem, to ensure that it is well-posed.

5

Dear @kamensky, after trying to implement the stabilized outlet boundary condition, I still observe blowup on the outer boundary because of spurious ricirculations. I guess my implementation is not correct with respect to my model. I am solving in ALE framework a Navier-Stokes system of the kind

\frac{\partial\mathbf{u}}{\partial t} + \boldsymbol{\Omega}\times\mathbf{u} + \left((\mathbf{u} - \mathbf{V} - \boldsymbol{\Omega}\times \mathbf{r} )\cdot\nabla \right)\mathbf{u}= \nu \Delta \mathbf{u} -\nabla p,

so the control volume is translating and rotating with a rigid motion of velocity \mathbf{V} +\boldsymbol{\Omega}\times \mathbf{r}, following a rigid body moving in the fluid, whose surface is the inner boundary in the middle of the control volume, and where a Dirichlet condition for velocity is applied. Currently the rest of the boundary, that is the outer boundary, is a unique spherical surface. I am trying this kind of boundary to have as smooth a boundary as possible, with no wedges.

I added to the weak form the following implicit term for each component of the velocity (I solve component for component) to try enforce outflow only:

-1 \int_{\Gamma} min ( \mathbf{u}^o \cdot \mathbf{n}, 0 ) u_i d\Gamma ,

where I use the old velocity in the dot product to linearize. In terms of form assembly I add to the matrix

-1 * assemble( inner ( min( \mathbf{u}^o \cdot \mathbf{n}, 0 ) * u , v ) ds ),

where u and v are trial and testfunction of the scalar velocity space. Also, on the outer boundary I want an open boundary condition like (-p + \nabla\mathbf{u}) \cdot \mathbf{n} = 0 \rightarrow \frac{\partial \mathbf{u}}{\partial n} = p, so I add to the RHS

assemble ( p^o * v * n[i] * ds),

where again I use old pressure since I’m using a projection scheme and pressure and velocity are solved for separately.

On the outer boundary I also fix pressure to 0 with a Dirichlet BC when I solve for the Poisson equation of the pressure correction. The simulation runs fine for long times, but spurious velocities arise at the outer boundary like in the image. These spurious flows increase in magnitude over time disrupting the entire simulation, which runs fine at the beginning. This is what happens to the z velocity (the body is rising in z direction):

screen1
screen12638×1600 190 KB

I treat the viscous terms implicitly with Crank Nicolson, and convective terms explicitly with Runge Kutta.

Is it possible that the ALE formulation might prevent the stability of the outflow? I am worried in particular by the rotation of the control volume. It’s strange that blowup doesn’t occur where the wake leaves the domain though.

1 Like
6

To apply this in the ALE description, you need to use the ALE advection velocity to determine inflow/outflow portions of the domain. You can use the term given in Equations (22) and (23) of this paper, where \hat{\mathbf{u}}^h is the mesh motion velocity field.

1 Like
7

Dear @kamensky, after adding the form

assemble( -gamma*ufl.Min(inner((\mathbf{u} - \mathbf{u}_{MESH}),n),Constant(0.0))*inner(u,v)*ds)

to the LHS of the weak form of Navier-Stokes I am solving, I can confirm that the instabilities at outflow are greatly reduced.

I tried the value of the coefficient \gamma = 0.2, but in my experience it was not enough to prevent an instability at ouflow eventually.

With the value \gamma=1 however, I was able to run the simulation for long times with no outflow instability at all.
I have yet to try intermediate values that maybe could still solve the problem with a less intrusive term, if I understood the article correctly.

Thanks!

8

I’m currently trying to implement a Runge-Kutta scheme to solve the wave equation, could you share an example of your approach?
Here my latest trial: Resolution of the wave equation using Runge-Kutta
Do you have any advice to give me?