I’m trying to simulate porous absorbers inside a waveguide(Parallel baffle silencer) and there are a few things that I can not figure out by myself.
Does the introduction of the porous absorber mean that I need a more sophisticated matched layer than
Z_s = rho_0 * c0
inner(1j*omega/Z_s*u, v) * ds(2)
- I want to calculate the transmission loss through my domain.
2a. How do I assemble the solution at surface 2?
2b. How do I project or interpolate my solution as to calculate the transmission loss?
import time
start_time = time.time()
import numpy as np
from mpi4py import MPI #Message Passing Interface
from petsc4py import PETSc #PETSc (the Portable, Extensible Toolkit for Scientific Computation)
#is a suite of data structures and routines for the scalable (parallel)
#solution of scientific applications modeled by partial differential equations.
from dolfinx import geometry
from dolfinx.io import gmshio, XDMFFile
from dolfinx.fem import FunctionSpace, Function, Constant
from dolfinx.fem.petsc import LinearProblem
import ufl
from ufl import dx, grad, inner, ds, Measure
import os
###change to complex mode source /usr/local/bin/dolfinx-complex-mode
##
# Extract base name of the input mesh file
mesh_file = "duct_cube_25_partition.msh"
mesh_base_name = os.path.splitext(mesh_file)[0]
# import the gmsh generated mesh
msh, cell, facet_tags = gmshio.read_from_msh(mesh_file, MPI.COMM_WORLD, 0, gdim=3) #From documentation: Returns:A triplet (mesh, cell_tags, facet_tags) with meshtags for associated physical groups for cells and facets.
# Construct output folder name
output_folder = "Results_" + mesh_base_name
# Create the output folder if it doesn't exist
if not os.path.exists(output_folder):
os.makedirs(output_folder)
# Construct output filenames using the base name of the input mesh file
output_prefix = mesh_base_name + "_"
#Creating frequency axis
#f_axis = np.arange(100, 8000, 100)
f_axis = [63, 125, 250, 500, 1000,2000,4000]
#f_axis = [8000]
#Define trial and test function
V = FunctionSpace(msh,("CG", 2)) #Continious Galerkin, Second order trial and test functions
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
#Defining air characteristics
c0= 340
rho_0 =1.225
omega = Constant(msh, PETSc.ScalarType(0)) #Constant over the domain, but changes with frequency_loop
k_0 = Constant(msh, PETSc.ScalarType(0))
k_c = Constant(msh, PETSc.ScalarType(0))
c_c = Constant(msh, PETSc.ScalarType(0))
Z_c = Constant(msh, PETSc.ScalarType(0))
rho_c = Constant(msh, PETSc.ScalarType(0))
#Delany-Bazley constants
C1 = 10.8
C2 = 0.70
C3 = 10.3
C4 = 0.59
C5 = 9.08
C6 = 0.75
C7 = 11.9
C8 = 0.73
#Defining absorber
rho_rw = 20.
R_f = 95.35 * (rho_rw **1.37)
# Normal velocity boundary condition
v_n = 0.01 #m/s
# Surface Impedance
Z_s = rho_0 * c0 #\rho*c means the surface is a perfect absorber
#Defining bilinear and linear forms
#
dx = Measure("dx", domain=msh, subdomain_data=cell) #All volumes need to be accounted for.
ds = Measure("ds", domain=msh, subdomain_data=facet_tags)
dI = Measure("dS", domain=msh, subdomain_data=facet_tags)
#Because dx is not defined by cell tags the whole domain is a single volume
#Air and porous volume
a = ( # Bilinear form definition
#Air and porous volume
(inner(grad(u)/rho_c, grad(v)) * dx - k_c**2/rho_c * inner(u, v) * dx(5))
-k_0**2 * inner(u/rho_0, v) * dx(4)
- inner(1j*omega/Z_s*u, v) * ds(2) # Far end of duct made a perfect absorber
- inner(1j*omega/Z_s*u, v) * ds(1) # Velocity boundary face made to be perfect absorber
)
'''a = ( # Bilinear form definition
#All volumes considered to be air
(1/rho_0)*(inner(grad(u), grad(v)) * dx - k_0**2 * inner(u, v) * dx
- inner(1j*omega/c0*u, v) * ds(2) # Far end of duct made a perfect absorber
- inner(1j*omega/c0*u, v) * ds(1)) # Velocity boundary face made to be perfect absorber
)'''
L = (1/rho_0)*inner(1j*-omega*rho_0*v_n, v) * ds(1) #Linear form, velocity boundary(sound source)
uh = Function(V) # uh is a function that belongs to the FunctionSpaceV ?
uh.name = "pressure"
solver = LinearProblem(a, L, u=uh, petsc_options={"ksp_type" : "preonly", "pc_type" : "lu", "pc_factor_mat_solver_type": "mumps"})
def frequency_loop(nf):
omega.value = 2*np.pi*f_axis[nf]
k_0.value = (2*np.pi*f_axis[nf]/c0)
Re_k = (1 + C1 * ((1000*rho_0*f_axis[nf])/R_f) ** (-C2)) * (2*np.pi*f_axis[nf]/c0)
Im_k = (1j * (C3 * ((1000*rho_0*f_axis[nf])/R_f) ** (-C4))) * (2*np.pi*f_axis[nf]/c0)
Re_Z = (1 + C1 * ((1000*rho_0*f_axis[nf])/R_f) ** (-C2)) * (rho_0*c0)
Im_Z = (1j * (C7 * ((1000*rho_0*f_axis[nf])/R_f) ** (-C8))) * (rho_0*c0)
Z_c.value = complex(Re_Z - Im_Z)
k_c.value = complex(Re_k + Im_k)
rho_c.value = complex((Z_c*k_c)/omega)
print (f_axis[nf])
c_c.value = complex(omega/k_c)
solver.solve()
for freq in f_axis:
with XDMFFile(msh.comm, os.path.join(output_folder, output_prefix + str(f_axis[nf])) + "_with_porous" ".xdmf", "w") as xdmf:
xdmf.write_mesh(msh)
xdmf.write_function(uh)
from multiprocessing import Pool
nf = range(0,len(f_axis))
if __name__== '__main__':
print("Computing...")
pool = Pool(3) #Something regarding splitting up processes
pool.map(frequency_loop,nf)
pool.close()
pool.join()
assemble(uh, ds(2))
end_time = time.time()
elapsed_time = end_time - start_time
print("Elapsed time:", elapsed_time, "seconds")
Mesh can be found here.
Running this code on a dolfinx container version 0.6.0.
docker pull dolfinx/dolfinx:v0.6.0
Regards
Ólafur Marteinsson