Issue with eigenmodes concentration in outlet

mathematics
1

Hi everyone,

I’m building a generalized FEM framework for linear stability (eigenvalue) analysis of incompressible Navier–Stokes flows.

My test case is the classic 2D flow around a cylinder (Taylor–Hood, P2–P1). I impose Dirichlet BCs at the inlet (input velocity) and cylinder (no-slip), and have tried both Neumann and Robin at the outlet. Top and bottom are set as free-flow (homogeneous Neumann).

The problem:
After assembling and solving the eigenvalue problem (SLEPc, shift-invert, LU preconditioned), all computed eigenmodes are highly localized at the outlet or at the bottom right corner.
The baseflow is physically reasonable and matches literature. Mesh refinements at the outlet and switching between Neumann and Robin have not helped.

Has anyone seen similar eigenmodes concentrating at an outlet/corner, or thinks that the issue sounds familiar?
Any suggestions would be very helpful.

Framework and case setup are at /ferdean/lsa-fw in GitHub if anyone wants to dig in further (ferdean-solver-fix branch). Still in development (documentation not updated)

Thanks for any ideas or pointers.

— Ferran

2

Example eigenmode:

em
em829×229 22.1 KB

With the following baseflow:

bf
bf965×259 22.8 KB

3

Are you solving a generalized eigenvalue problem with the mass on the right-hand-side?

4

Yes, exactly, generalized non-Hermitian eigenvalue problem (Lx = λMx).

Thanks for the answer!

5

Could you state the weak form that you’ve implemented (code or math) that defines L? Some integration by parts choices might imply awkward BCs at the outflow (that for a transient simulation don’t matter, but for eigenmodes would).

6

Of course. The weak form I implemented is the linearized incompressible Navier–Stokes operator around a stationary base flow (u_b):

image
image379×479 18.8 KB

Which corresponds to this code:

class VariationalForms:
    """Collector for variational forms for linearized incompressible Navier-Stokes equations."""

    @staticmethod
    def mass(u: Argument, v: Argument) -> Form:
        return inner(u, v) * dx

    @staticmethod
    def convection(u: Argument, v: Argument, u_base: dfem.Function) -> Form:
        return inner(dot(u_base, nabla_grad(u)), v) * dx

    @staticmethod
    def shear(u: Argument, v: Argument, u_base: dfem.Function) -> Form:
        return inner(dot(u, grad(u_base)), v) * dx

    @staticmethod
    def pressure_gradient(p: Argument, v: Argument) -> Form:
        return -inner(p, div(v)) * dx

    @staticmethod
    def viscous(u: Argument, v: Argument, re: float) -> Form:
        return (1.0 / re) * inner(grad(u), grad(v)) * dx

    @staticmethod
    def divergence(u: Argument, q: Argument) -> Form:
        return inner(div(u), q) * dx

    @staticmethod
    def stiffness(u: Argument, v: Argument) -> Form:
        return inner(grad(u), grad(v)) * dx

# Later in the code...

self._u, self._p = TrialFunctions(self._spaces.mixed)
self._v, self._q = TestFunctions(self._spaces.mixed)

# ...

shear = VariationalForms.shear(self._u, self._v, self._u_base)
convection = VariationalForms.convection(self._u, self._v, self._u_base)
viscous = VariationalForms.viscous(self._u, self._v, self._re)
gradient = VariationalForms.pressure_gradient(self._p, self._v)
divergence = VariationalForms.divergence(self._u, self._q)

full_form = shear + convection + viscous + gradient + divergence

for marker, alpha, _ in self._robin_bcs:
    full_form += alpha * inner(self._u, self._v) * self._ds(marker)

mat = assemble_matrix(dfem.form(full_form, dtype=Scalar), bcs=self._dirichlet_bcs)
mat.assemble()

If no pressure Dirichlet BCs are present, I’m also attaching the pressure nullspace.

Note that I voluntarily flip the sign of the pressure term, which helped me ensuring symmetry in other assemblers (like the simple Stokes flow).

Same procedure applies for mass matrix.

I appreciate a lot your time @Stein !

7

I believe the Robin BC should be unstable at outflow boundaries, is that one active?

8

I added them expecting it to work as a “sponge” that would absorb spurious perturbations at the outlet; but you are right, I was not successful.

But no, in the above figures they were not active.

For the baseflow, I used Dirichlet velocity at the inlet; and at the cylinder (no slip), and homogeneous Neumann at the rest (free-flow). Same for the perturbation flow (with the difference that the Dirichlet velocity BC at the inlet for the perturbation flow is also homogeneous).

9

One last thing that I can think of is that the shifting is not working correctly. I had a student work on a similar problem a while back, and if I remember correctly, that was quite a bit of a challenge to get to work robustly and truly give the initial ones.

Without context, but hopefully somewhat helpfull nonetheless:

# Set the options for the SLEPc solver with spectral transformation method Shift-Invert
eps = SLEPc.EPS().create(mesh.comm)
eps.setOperators(-A,M) # Solve A*x + lambda*M*x = 0
eps.setProblemType(SLEPc.EPS.ProblemType.PGNHEP) 
tol = 1e-6
eps.setTolerances(tol=tol)
st = eps.getST()
st.setType(SLEPc.ST.Type.SINVERT) 
eps.setTarget(0.1)
eps.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_REAL)
nev = 4 # nev is the # of eig. values we are looking for
eps.setDimensions(nev=nev) 
eps.solve()