Low Frequency Analysis for Small Room Acoustics

dolfinx
1

Hi all,

I am trying to implement a frequency domain study in FEniCSx 2021.1.0 (within Docker) for small room acoustics. For now I am giving it a soft start with the computation of the pressure field at a single frequency (first eigenmode, 34.3 Hz) within a cuboid room 5x4x3 m, for which the analytical solution is known. All the room surfaces are considered to be rigid, sound source is a monopole.

The formulation of the Helmholtz problem would be the following:

\int_{\Omega} \nabla u \nabla v \,dx \,-k^2\int_{\Omega} u v \,dx \, = \, \int_{\Omega} f v \,dx

Where f = j\omega\rho_0Q\, \delta(x-x_0) is the monopole source term, k is the wavenumber, \omega is the angular frequency, \rho_0 is air density, Q is the volume velocity of the source.

So far, I have managed to get a similar spatial distribution of the pressure field when compared to the analytical solution, however I notice that the pressure values change dramatically with the element number.

Here is the code:

import dolfinx
import numpy as np
import math
import meshio
import pyvista 
import warnings
warnings.filterwarnings('ignore')
import dolfinx.plot
from mpi4py import MPI
import ufl
import gmsh

# Room Geometry
Lx = 5.0
Ly = 4.0
Lz = 3.0
room_dim = [Lx, Ly, Lz]

# Monopole Source Position
Sx = 1
Sy = 1
Sz = 1
src = [Sx, Sy, Sz]

# Frequency parameters
freq = 34.3
omega = freq*2*np.pi
c = 343
rho_0 = 1.21
k = 2*np.pi*freq/c 

# Mesh parameters
n_el = 15 #number of elements per wavelength
mesh_size = c/freq/n_el

## Create Mesh
gmsh.initialize()
gmsh.model.occ.add_box(0, 0, 0, Lx, Ly, Lz, 1)
gmsh.model.occ.add_point(Sx, Sy, Sz, mesh_size, 28)
gmsh.model.occ.synchronize()
gmsh.model.mesh.embed(0, [28], 3, 1)
gmsh.model.addPhysicalGroup(3, [1], 1)
gmsh.option.setNumber("Mesh.MeshSizeMax", mesh_size)
gmsh.model.mesh.generate(3)

## Convert mesh to dolfinx
if MPI.COMM_WORLD.rank == 0:
    geometry_data = extract_gmsh_geometry(gmsh.model)
    topology_data = extract_gmsh_topology_and_markers(gmsh.model)
    gmsh_cell_type = [key for key in list(topology_data.keys())][0]
    properties = gmsh.model.mesh.getElementProperties(gmsh_cell_type) 
    name, dim, order, num_nodes, local_coords, _ = properties
    cells = topology_data[gmsh_cell_type]["topology"]
    cell_id, num_nodes = MPI.COMM_WORLD.bcast([gmsh_cell_type, num_nodes], root=0)

from dolfinx.cpp.io import perm_gmsh
from dolfinx.cpp.mesh import to_type
from dolfinx.mesh import create_mesh

ufl_domain = ufl_mesh_from_gmsh(cell_id, 3)
gmsh_cell_perm = perm_gmsh(to_type(str(ufl_domain.ufl_cell())), num_nodes)
cells = cells[:, gmsh_cell_perm]

mesh = create_mesh(MPI.COMM_WORLD, cells, geometry_data[:, :3], ufl_domain)

## Define monopole source
def delta_dirac(src, rec):
    return np.any(np.isclose(src,rec))

## Define FEM problem
V = dolfinx.FunctionSpace(mesh, ("CG", 1))
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
f = dolfinx.Function(V)
f.interpolate(lambda x: [1j*rho_0*omega*Q*delta_dirac(src, item) for item in np.transpose(x)])

a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx - k**2 * ufl.inner(u, v) * ufl.dx
L = ufl.inner(f, v) * ufl.dx

# Compute solution
uh = dolfinx.fem.Function(V)
uh.name = "u"
problem = dolfinx.fem.LinearProblem(a, L, u=uh)
problem.solve()

See the images below, the first one obtained with 15 elements / wavelegth, the second with 10 elements / wavelength. It looks like that the number of elements acts as scale factor on pressure values, but can’t figure out in what extents.

10
10990×399 180 KB

Perhaps is there any issue with my code? Or is there a better alternative to formulate the problem?

Apologies in advance should I have missed useful pieces of information.

Many thanks,
Andrea

1 Like
2

Assuming Q is constant, the magnitude of the Dirac delta in the source term is mesh-dependent. Right now, the function delta_dirac is the characteristic function of a singleton set containing the point src, not the Dirac delta at src. What you want is for the integral of the Dirac delta to be 1, not its point evaluation. You could try normalizing your current implementation, with something like the following:

def src_char_func(src, rec):
    return np.any(np.isclose(src,rec))
dirac_delta_tmp = dolfinx.Function(V)
dirac_delta_tmp.interpolate(lambda x: [src_char_func(src, item)
                                                for item in np.transpose(x)])
int_dirac_delta_tmp = dolfinx.fem.assemble_scalar(dirac_delta_tmp*ufl.dx)
dirac_delta = dirac_delta_tmp / dolfinx.Constant(mesh,int_dirac_delta_tmp)
f = 1j*rho_0*omega*Q*dirac_delta
4 Likes