Modal analysis of an open cavity : How can I implement non reflective boundary conditions?

General
1

Greetings everyone,

I am currently trying to perform an acoustic modal analysis on an open cavity. I am using the slepc4py module with FenicsX and I am trying to adapt this tutoriel :
Introduction to Acoustic Modal Analysis - Undabit Acoustic Simulations ( solving the helmholtz equation)

I am starting simple with a rectangular open cavity generated with gmsh.
However I cannot figure out how to implement non reflective boundary conditions for the opened boundary. I thought about using simple PMLs (using a complex wave number near the boundary to damp the wave ) but I couldn’t find a way to implement it with the eigensolver. Does anyone know how to implement it ? Or any other implementation of a non reflective boundary condition ?

I also tried to generate a big external mesh but it creates a lot of undesirable low frequency modes and I don’t know how to get rid of them. If you know how to get only the modes that resonnate into the cavity that could help too.

My last resort was to try Robin boundary conditions as presented in the FenicsX tutoriel :
https://jsdokken.com/dolfinx-tutorial/chapter3/robin_neumann_dirichlet.html
but it doesn’t seems to get the desired result.
Here is a copy of my code with Robin condition applied :

```python
# Cool imports

import dolfinx as dfx
from dolfinx.fem import Function, FunctionSpace, form
from dolfinx.fem.petsc import assemble_matrix
from slepc4py import SLEPc
from ufl import dx, grad, inner, TrialFunction, TestFunction, FacetNormal, Measure
from dolfinx.io import XDMFFile, gmshio
from mpi4py import MPI
import numpy as np

from dolfinx.mesh import  locate_entities, meshtags






#   Resolution parameters

N_ext = 600  #number of acoustic modes extracted1
deg = 1  #degree of the polynomial resolution

mesh, cell_tags, facet_tags = gmshio.read_from_msh(mesh3, MPI.COMM_WORLD, 0, gdim = 2)
n = FacetNormal(mesh)

#   Air parameters

rho = 1.225  #air density
c = 341  #soundspeed
Z_inv = 1/(rho*c)



#   Test & trial functions


V = FunctionSpace(mesh, ("CG", deg))
u = TrialFunction(V)
v = TestFunction(V)



def B_air(x) :
    
	return np.isclose(x[0]**2 + x[1]**2, 0.3**2)

n = FacetNormal(mesh)
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
marker = 1
locator = B_air
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)




k_form = form(inner(grad(u), grad(v))*dx + Z_inv*inner(u, v)*ds(marker))
m_form = form(inner(u,v)*dx )
bcs = []
K = assemble_matrix(k_form, bcs)  #where we input the dirichlet bc if we had one in the bracket
M = assemble_matrix(m_form, bcs)

K.assemble()
M.assemble()




#   Setup of the eigensolver

solver = SLEPc.EPS().create()
solver.setDimensions(N_ext)  #number of acoustic modes extracted
solver.setProblemType(SLEPc.EPS.ProblemType.GHEP)  #setting up the type of problem to solve, here a General Hermitian Eigenvalue Problem

st = SLEPc.ST().create()  #initializing spectral transformation module
st.setType(SLEPc.ST.Type.SINVERT) #setting up shift and invert algorithm
st.setShift(0.1)  #setting up the shift
st.setFromOptions()  #finalizing setting up from options

solver.setST(st)  #assigning these settups to our solver

solver.setOperators(K,M)  #adjusting to stiffness and mass matrixs

Thanks you all for reading, have a nice day !