Pressure solution is not unique in splitting scheme for Navier Stokes equations

[mekong]

Hello,
In this Fenics example for the Navier Stokes equations, in the second step of the splitting scheme the boundary conditions (BCs) for the pressure correction \phi are

• a Dirichlet BC in the boundary \partial \Omega_N. This results from the fact that the pressure is constrained by a Dirichlet BC on \partial \Omega_N.
• a Neumann BC in the boundary \partial \Omega_D (where Dirichlet BCs are imposed on the velocity).

Now consider a slightly different problem: Navier Stokes equations in a disk, in which I impose, on the disk boundary \partial \Omega

• a traction BC \sigma_{ij}n_j = ... where \sigma is the stress tensor and \hat n the unit normal
• a BC on the normal velocity v_i n_i = ...

Proceeding along the lines as in the Fenics example, the second step is
\nabla^2 \phi = \frac{\rho}{\delta t} \nabla \cdot u^* \qquad\text{in } \Omega
with BC
\hat{n} \cdot \vec{\nabla} \phi = 0 on \partial \Omega.

The solution is not unique, and defined modulo an additive constant. How can I get the correct pressure field and get around this non-uniqueness?

Thank you

[Stein]

@mekong

a traction BC \sigma_{ij}n_j = … where \sigma is the stress tensor and \hat n the unit normal

So this is your Neumann boundary. Then there would be a Dirichlet condition there, right?

You just write out \sigma_{ij} n_j = 2\mu[\nabla^s u]_{ij} \cdot n_j - p n_i.

(This is what you should always do, really. Often one can assume the first part is small, which is why the “free outflow” condition is often misunderstood as a pressure condition. )

[mekong]

@Stein

Then there would be a Dirichlet condition there, right?

Yes, like I said in my original post, ‘a BC on the normal velocity …’’

@Stein

You just out \sigma_{ij} n_j = 2\mu[\nabla^s u]_{ij} \cdot n_j - p

What do you mean by ‘you just out’? Also, \nabla^s is not defined.

The variational problem that I posted determines \phi up to an additive constant, and I would like to know how to fix that constant.

Thank you!

[Stein]

@mekong

What do you mean by ‘you just out’? Also, \nabla^s is not defined.

Typo, fixed now. \nabla^s is the symmetric gradient (the strain rate operator). \sigma = 2\mu \nabla^s u - p I is how the stress tensor is defined for an incompressible newtonian fluid.

@mekong

Yes, like I said in my original post, ‘a BC on the normal velocity …’’

This is not related to the normal component boundary condition though… As you indicated in your initial post, the Dirichlet condition for the pressure problem occurs at the Neumann (or really, the natural) boundaries for the Navier-Stokes problem. For the Navier-Stokes equations, the natural boundary condition is:

\boldsymbol{\sigma} \cdot \boldsymbol{n} = \boldsymbol{t}_N

i.e.,

2\mu \nabla^s \boldsymbol{u} \cdot \boldsymbol{n} - p \boldsymbol{n} = \boldsymbol{t}_N

Or, by dotting with \boldsymbol{n}

p = (2\mu \nabla^s \boldsymbol{u} \cdot \boldsymbol{n}) \cdot \boldsymbol{n} - \boldsymbol{t} \cdot \boldsymbol{n}

Which gives you your Dirichlet expression for p. That Dirichlet expression fixes the constant you’re referring to.

[mekong]

@Stein

That Dirichlet expression fixes the constant you’re referring to.

Yes. So what if I

• solve for \phi by adding a temporary Dirichlet condition on a vertex on the mesh, setting for example \phi = 0 on that vertex. This temporary boundary condition selects one out of the many degenerate solutions, and prevents the solve from going to ‘large’ solutions.
• The resulting pressure is p_{n-1/2} = p_{n-3/2} + \phi + C, and I determine C from the condition 𝑝 =(2𝜇∇𝑠𝒖 ⋅𝒏) ⋅𝒏 −𝒕_N ⋅𝒏?

[Stein]

I don’t think that works. Applying a Neumann condition for \phi (the surrogate pressure variable) on \partial\Omega_N is just inconsistent. As you’ve said yourself in your first post, that’s were there should be a Dirichlet condition on the pressure problem…

Also, now you’d have to solve another PDE for C.

Why wouldn’t you just determine the condition on \phi from that last statement, per my earlier suggestion? Am I missing something?

[mekong]

@Stein

Applying a Neumann condition for \phi (the surrogate pressure variable) on \partial\Omega_N is just inconsistent

Thank you for your reply. I don’t know what you mean by ‘surrogate pressure variable’, here \phi = p^{n-1/2} - p^{n-3/2}, see here.

\phi must satisfy the Neumann BC above on \partial \Omega: In fact, from the splitting scheme (see here, the equation after (35)) we have
\frac{\rho}{\Delta t}({\bf v}^n - \overline{{\bf v}}) = - \vec{\nabla}\phi
where {\bf v}^n \equiv {\bf v}(t = t_n) and \overline{{\bf v}} is the approximate velocity. By taking the dot product of the equation above with \hat{n} and using the BC that constrians the normal component of the velocity (see my original post)

\hat{n} \cdot {\bf v}^n = \hat{n} \cdot \overline{{\bf v}} = v_0 \text{ on } \partial \Omega,

we obtain that \phi must satisfy

\hat{n} \cdot \vec{\nabla} \phi = 0 \text{ on } \partial \Omega.

@Stein

Also, now you’d have to solve another PDE for C.

Here C is a constant, so I need to solve an algebraic equation for it, not a PDE.

[Stein]

Aah, you have \boldsymbol{v}\cdot\boldsymbol{n} on the entire domain boundary? I got confused by your leading point

@mekong

… , in which I impose, on the disk boundary \partial \Omega

• a traction BC \sigma_{ij}n_j = … where \sigma is the stress tensor and \hat n the unit normal

But I suppose you meant a condition t_i \sigma_{ij}n_j = ... (with t the tangential vector).

If you have an impermeability condition all around, then your pressure is indeed only defined up to a constant. The same thing applies to \phi. With pressure surrogate I meant that it holds some equilvalently to the pressure.

You can:

• Add the condition that \phi has zero mean (requires an additional lagrange multiplier constraint in the real function space from scifem, Real function spaces — scifem )
• Prescribe a single dof with a Dirichlet condition per your suggestion. I believe this might have some negative consequences in terms of solver performance though
• Use iterative solvers (these typically struggle less with nullspaces), and postprocess \phi to have zero mean (see these helper functions src/fx4f/fluid_mechanics/postprocess.py · main · FEniCSx4Flow / Fx4F · GitLab )

[mekong]

Great, we agree then, thank you for your help!