Hello everybody,
I am currently trying to solve acoustic scattering problems and seek some advice from the community.
The goal is to solve the Helmholtz equation in a domain bounded by a PML layer, excited by a plane wave, or a spherical wave outside the simulation domain. The simulation domain can contain multiple penetrable obstacles with layered media of distinct properties (sound speed, density).
To streamline the simulation, I’m attempting a Total-Field/Scattered-Field (TF/SF) approach. This method should, theoretically, simplify handling the media inside the simulation domain by applying the incident field in a specific zone between the PML and the Simulation Domain. Here’s an illustrative test case that should result in a plane wave traveling through the whole domain (see attached image).
I understand the necessity of ensuring total pressure continuity at the TF/SF interface. It appears that employing Lagrange multipliers could effectively enforce this condition. Despite this, the practicalities of setting up the correct spaces and applying the continuity condition with the presence of an incident field remain elusive to me.
I have been using the FETI.ipynb notebook in the iFEM repo, as well as the PML example for the Helmholtz case in the documentation, to try to solve the problem. Unfortunately, integrating these insights into a functional model has been challenging.
I’m reaching out for any guidance, insights, or examples you might share. Specifically, I’m looking for advice on:
- Correctly setting up the simulation spaces for the TF/SF approach.
- Applying continuity conditions with Lagrange multipliers accounting for the incident field.
I would very much appreciate any insights that could guide me in the right direction. I have been struggling on this topic for far longer than I would like to admit, moving away from scratch-built FEM solutions and starting to work with NGSolve’s elegant framework.
Thank you for any support,
Luiz Alvim
Do you mean something like this?
from netgen.occ import *
from ngsolve import *
obj = Rectangle(0.2, 0.2).Face().Rotate(Axis((0,0,0), Z), 30)
c1 = Circle((0,0),1).Face()
c2 = Circle((0,0),1.2).Face()
c2.faces.name = "pml"
c2 -= c1
c1.faces.name = "inner"
obj.edges.name = "object"
c1 -= obj
shape = Glue([c1,c2])
geo = OCCGeometry(shape, dim=2)
mesh = Mesh(geo.GenerateMesh(maxh=0.05))
Draw(mesh)
omega = 2 * pi * 3
incident = exp(1j * omega * y)
Draw(incident, mesh, "incident")
mesh.SetPML(pml.Radial(rad=1.1, alpha=1j, origin=(0,0)), "pml")
fes = H1(mesh, order=3, complex=True, dirichlet="object")
u,v = fes.TnT()
a = BilinearForm(fes)
a += (grad(u) * grad(v) - omega**2 * u * v) * dx
# a += 1j * lam * u * v * ds("pml")
gfu = GridFunction(fes)
gfu.Set(-incident, definedon=mesh.Boundaries("object"), dual=True)
r = gfu.vec.CreateVector()
with TaskManager():
a.Assemble()
r.data = -a.mat * gfu.vec
gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * r
Draw(gfu + incident, mesh, "total")
Draw(gfu, mesh, "scattered")
not exactly sure how you mean it with this incident domain. Can you write down the formulation for it or give a reference?
Hi Christopher,
I really appreciate the response. As far as I understand, defining the incident field at the scatterer’s boundary requires that the analytical source term be changed depending on the problem at hand. To avoid that, the incident field could be applied to a region that always has the same properties and let it radiate into the simulation domain. This way, I could solve layered penetrable objects without having to mess with the source term. Does it make sense?
I’ve been using as the primary reference COMSOL’s Background Pressure Field
Specifically addressing the phrase, “In a model where the background pressure field is not defined on all acoustic domains (or it is different), continuity is automatically applied in the total field p_t on interior boundaries between domains.”
For this enforced continuity, a space is created at the discontinuity interface, and the following expression is evaluated:
Constraint:
up(p_t) - down(p_t)
Constraint Force:
test(up(p_t) - down(p_t))
where up and down are operators that evaluate the variable at each side of the interface, which I imagine is equivalent to using .Other() or something like (u - uhat) where uhat is defined on the interface. In this case, the contribution from the incident field should be accounted for in p_t at the region with the incident field.
My attempt at solving the problem is attached, but clearly, the “b” matrix is wrong, and I would assume other things are too.
tf_sf_acoustic_scattering.ipynb (11.5 KB)
Another source of information for a problem similar to mine is 8.3.4.3. Analytic Wave Sources, in this case, the incident field is applied to the PML/SimulationDomain interface, a different approach that could also work on my case.
Total pressure can be split into scattered pressure and incident pressure on the interface between the general acoustic region and the PML region. This can be expressed as:
Thank you for your support.
