2D Helmholtz problem coupling FEM and BEM#

Keywords: PDE:Helmholtz Method:FEM Method:BEM Method:Coupling FE:Lagrange Solver:Direct

We want to solve the acoustic propagation of a plane wave in a heterogeneous medium. In order to do that, we distinguish a domain $Ω$ that is heterogeneous, its boundary $Γ$ and the exterior domain $Ω_{ext}$ that is homogeneous:

Figure made with TikZ

Fig. 17 Domains for computation: $Ω$ the heterogeneous medium, $Ω_{ext}$ the homogeneous exterior domain and $Γ = \partial Ω$ .

Figure made with TikZ

Fig. 18 Domains for computation: $Ω$ the heterogeneous medium, $Ω_{ext}$ the homogeneous exterior domain and $Γ = \partial Ω$ .

We solve

\begin{array}{r} {\begin{cases} Δ u (x) + k^{2} η^{2} (x) u (x) = 0 & in R^{2} \\ u (x) = - u_{i} (x) & on Γ \end{cases} \end{array}

with $η (x) = 1$ in $Ω_{ext}$ , and $η (x)$ that can vary in $Ω$ , and finally with $u_{i} = e^{i k x}$ .

We will use: $Ω = {[- 0.5, 0.5]}^{2}$ and

\begin{array}{r} η (x) = {\begin{cases} \exp (- (x_{1}^{2} - 0.25) * (x_{2}^{2} - 0.25) / (2. * 0.05)), when max (x_{1}, x_{2}) < 0.5 . \\ 1 otherwise \end{cases} \end{array}

../_images/helmholtz2d_FEM_BEM_eta.png — Fig. 19 $η (x)$ in $Ω \cup Ω_{ext}$ #

We decompose the problem in a coupled system of two equations:

in the FEM part, the solution solves the following equation:

$Δ u + k^{2} η^{2} u = 0$

which gives the variational formulation

$\begin{array}{r} | \begin{array}{c} Find u \in H^{1} (Ω) such that : \\ \int_{Ω} \nabla u (x) \cdot \nabla \bar{v} (x) d x - k^{2} \int_{Ω} η^{2} (x) u (x) \bar{v} (x) d x - \int_{Γ} λ (x) \bar{v} (x) d x = 0, \forall v \in H^{1} (Ω) \end{array} \end{array}$

with $λ = \frac{\partial u}{\partial n}$ is the normal trace of $u$ on $Γ$ .
in the BEM part, we solve

$\begin{array}{r} {\begin{cases} Δ u + k^{2} u = 0 & in Ω_{ext} \\ u = - u_{i} & on Γ . \end{cases} \end{array}$

The scattered field verifies

$u_{s} (x) = - S_{Γ} λ (x) + K_{Γ} u (x), x \in Ω_{ext}$

with $u$ the total field solution of the equation and $λ$ the normal trace of $u$ on $Γ$ , $S_{Γ}$ and $K_{Γ}$ are the single and double layer boundary potentials

$\begin{aligned} S_{Γ} ϕ (x) & = \int_{Γ} G (x, y) ϕ (y) d y, \\ K_{Γ} ϕ (x) & = \int_{Γ} \frac{\partial G (x, y)}{\partial n_{y}} ϕ (y) d y, \end{aligned}$

and

$G (x, y) = \frac{e^{i k ‖ x - y ‖}}{4 π ‖ x - y ‖} .$

Since $u_{s} = u - u_{i}$ , and taking the limit when $x$ goes to $Γ$ , we obtain the integral equation

$(\frac{I}{2} - K_{Γ}) u (x) + S_{Γ} λ (x) = u_{i} (x), x \in Γ .$

The resulting variational formulation for the BEM part is then

$\begin{array}{r} | \begin{array}{c} Find u \in H^{1 / 2} (Γ) and λ \in H^{- 1 / 2} (Γ) such that : \\ \frac{1}{2} \int_{Γ} u (x) \bar{τ} (x) d x - \int_{Γ \times Γ} u (y) \frac{\partial G (x, y)}{\partial n_{y}} \bar{τ} (x) d y d x + \int_{Γ \times Γ} λ (y) G (x, y) \bar{τ} (x) d y d x \\ = \int_{Γ} u_{i} (x) \bar{τ} (x) d x, \forall τ \in H^{1 / 2} (Γ) . \end{array} \end{array}$

By adding the variational formulations relatives to the two linked problems, we obtain the final variational formulation.

Finally, the solution is obtained directly from $u$ for the FEM part, and we need to compute the integral representation to obtain $u_{s}$ , the scattered field, and then to add the incident field to obtain the total field for this problem.

The last step is to merge the FEM solution in $Ω$ and the BEM solution in $Ω_{ext}$ to obtain a solution on the whole domain $Ω \cup Ω_{ext}$ to simplify the visualization.

../_images/helmholtz2d_FEM_BEM_sol.png — Fig. 20 Solution of the FEM-BEM problem.#

The code of this example follows:

#include "xlife++.h"
using namespace xlifepp;
using namespace std;

// find = eta(x)
Real find(const Point & M, Parameters & pa = defaultParameters)
{
  Real res=1.;
  if (std::max(std::abs(M[0]), std::abs(M[1])) < 0.5)
    res=std::exp(-((M[0]*M[0]-0.25)*(M[1]*M[1]-0.25))/(2.*0.05));
  return res;
}
Real eta2(const Point & M, Parameters & pa = defaultParameters)
{
  Real tmp=find(M);
  return tmp*tmp;
}
Complex g1(const Point& M, Parameters& pa = defaultParameters)
{
  Real k=real(pa("k"));
  Point d(1., 0.);
  return exp(i_*(k*dot(M, d)));
}

int main(int argc, char** argv)
{
  init(argc, argv, _lang=en);   // mandatory initialization of xlifepp
  verboseLevel(10);
  Real k=10.;
  // meshing
  Real hsize=(2*pi_/k)/15.;
  SquareGeo sp(_center=Point(0., 0.), _length=1., _hsteps=hsize, _domain_name="Omega", _side_names="Gamma");
  Mesh m1=Mesh(sp, _shape=triangle, _generator=gmsh);
  Domain omega = m1.domain("Omega");
  Domain gamma = m1.domain("Gamma");
  theCout << "Mesh size = " << hsize << eol;
  theCout << "Number of Triangles = " << m1.nbOfElements() << eol;
  // defining parameter and kernel
  Parameters pars;
  pars << Parameter(k, "k");
  Kernel G=Helmholtz2dKernel(pars);
  Function finc(g1, pars);
  // defining space, unknown and test function
  Space V1(_domain=omega, _interpolation=P1, _name="V1", _notOptimizeNumbering);
  Space V0(_domain=gamma, _interpolation=P1, _name="V0", _notOptimizeNumbering);
  Unknown u1(V1, _name="u1"); TestFunction v1(u1, _name="v1");
  Unknown l0(V0, _name="l0"); TestFunction lt0(l0, _name="lt0");
  theCout << "Nb dofs BEM= " << V0.nbDofs() << " Nb dofs FEM= " << V1.nbDofs() << eol;
  // defining bilinear and linear form
  IntegrationMethods ims(_method=Duffy, _order=15, _bound=0., _quad=defaultQuadrature, _order=12, _bound=1.,
                         _quad=defaultQuadrature, _order=10, _bound=2., _quad=defaultQuadrature, _order=8);
  BilinearForm blf=intg(omega, grad(u1)|grad(v1))-k*k*intg(omega, eta2*u1*v1) - intg(gamma, l0*v1) + 0.5*intg(gamma, u1*lt0)
                  - intg(gamma, gamma, u1*ndotgrad_y(G)*lt0, _method=ims) + intg(gamma, gamma, l0*G*lt0, _method=ims);
  LinearForm lf=intg(gamma, finc*lt0);
  // computing FEM/BEM matrix and right hand side vector
  TermMatrix lhs(blf, _name="lhs");
  TermVector rhs(lf);
  // solving linear system using direct method
  TermVector sol=directSolve(lhs, rhs);

  // Representing the solution FEM and BEM
  SquareGeo Sint(_center=Point(0., 0.), _length=1, _hsteps=hsize, _domain_name="S_int");
  SquareGeo Sext(_center=Point(0., 0.), _length=3, _hsteps=1.5*hsize, _domain_name="S_ext");
  Mesh mrep(Sext+Sint, _shhape=triangle, _generator=gmsh);
  Domain S_ext=mrep.domain("S_ext"), S_int=mrep.domain("S_int");
  Domain S=merge(S_ext, S_int, "S");
  Space Vrep(_domain=S, _interpolation=P1, _name="Vrep", _notOptimizeNumbering);
  Unknown ur(Vrep, _name="ur");
  Function Find(find, pars);
  TermVector findex(ur, S, Find);
  saveToFile("findex", findex, _format=vtu); // Representing eta
  TermVector Uint=interpolate(ur, S_int, sol(u1)); // FEM solution (total field)
  saveToFile("Uint", Uint, _format=vtu);

  // Representing of the BEM part
  IntegrationMethods imr(_quad=GaussLegendre, _order=20, _bound=1., _quad=GaussLegendre, _order=10, _bound=2.,
                         _quad=GaussLegendre, _order=5);
  TermVector Uext = - integralRepresentation(ur, S_ext, intg(gamma, G*sol(l0), _method=imr))
                    + integralRepresentation(ur, S_ext, intg(gamma, ndotgrad_y(G)*sol(u1), _method=imr));

  TermVector Uinc(ur, S_ext, finc);
  saveToFile("Uinc", Uinc, _format=vtu); // Incident field
  saveToFile("Uext", Uext, _format=vtu); // scattered field in exterior domain
  TermVector Uext_t = Uext + Uinc;
  saveToFile("Uext_t", Uext_t, _format=vtu); // Total field in exterior domain
  TermVector U=merge(Uint, Uext_t); // Merged FEM and BEM solutions
  saveToFile("U", U, _format=vtu);
  theCout << "Program finished" << eol;
  return 0;
}