FreeFEM Documentation on GitHub

stars - forks

Acoustics

Summary : Here we go to grip with ill posed problems and eigenvalue problems

Pressure variations in air at rest are governed by the wave equation:

\[{\partial^2 u \over \partial t^2} - c^2 \Delta u =0\]

When the solution wave is monochromatic (and that depends on the boundary and initial conditions), \(u\) is of the form \(u(x,t)=Re(v(x) e^{ik t})\) where \(v\) is a solution of Helmholtz’s equation:

\[\begin{split}\begin{array}{rcl} k^{2}v + c^{2}\Delta v &= 0 &\hbox{ in } \Omega\\ \frac{\partial v}{\partial\boldsymbol{n}}|_\Gamma &= g & \end{array}\end{split}\]

where \(g\) is the source.

Note the “+” sign in front of the Laplace operator and that \(k > 0\) is real. This sign may make the problem ill posed for some values of \(\frac c k\), a phenomenon called “resonance”.

At resonance there are non-zero solutions even when \(g=0\). So the following program may or may not work:

 1// Parameters
 2real kc2 = 1.;
 3func g = y*(1.-y);
 4
 5// Mesh
 6border a0(t=0., 1.){x=5.; y=1.+2.*t;}
 7border a1(t=0., 1.){x=5.-2.*t; y=3.;}
 8border a2(t=0., 1.){x=3.-2.*t; y=3.-2.*t;}
 9border a3(t=0., 1.){x=1.-t; y=1.;}
10border a4(t=0., 1.){x=0.; y=1.-t;}
11border a5(t=0., 1.){x=t; y=0.;}
12border a6(t=0., 1.){x=1.+4.*t; y=t;}
13
14mesh Th = buildmesh(a0(20) + a1(20) + a2(20)
15    + a3(20) + a4(20) + a5(20) + a6(20));
16
17// Fespace
18fespace Vh(Th, P1);
19Vh u, v;
20
21// Solve
22solve sound(u, v)
23   = int2d(Th)(
24        u*v * kc2
25      - dx(u)*dx(v)
26      - dy(u)*dy(v)
27   )
28   - int1d(Th, a4)(
29        g * v
30   )
31   ;
32
33// Plot
34plot(u, wait=1, ps="Sound.eps");

Results are on Fig. 11. But when \(kc2\) is an eigenvalue of the problem, then the solution is not unique:

  • if \(u_e \neq 0\) is an eigen state, then for any given solution \(u+u_e\) is another solution.

To find all the \(u_e\) one can do the following :

 1// Parameters
 2real sigma = 20; //value of the shift
 3
 4// Problem
 5// OP = A - sigma B ; // The shifted matrix
 6varf op(u1, u2)
 7   = int2d(Th)(
 8        dx(u1)*dx(u2)
 9      + dy(u1)*dy(u2)
10      - sigma* u1*u2
11   )
12   ;
13
14varf b([u1], [u2])
15   = int2d(Th)(
16        u1*u2
17   )
18   ; // No Boundary condition see note \ref{note BC EV}
19
20matrix OP = op(Vh, Vh, solver=Crout, factorize=1);
21matrix B = b(Vh, Vh, solver=CG, eps=1e-20);
22
23// Eigen values
24int nev=2; // Number of requested eigenvalues near sigma
25
26real[int] ev(nev);  // To store the nev eigenvalue
27Vh[int] eV(nev);    // To store the nev eigenvector
28
29int k=EigenValue(OP, B, sym=true, sigma=sigma, value=ev, vector=eV,
30   tol=1e-10, maxit=0, ncv=0);
31
32cout << ev(0) << " 2 eigen values " << ev(1) << endl;
33v = eV[0];
34plot(v, wait=true, ps="eigen.eps");
Amplitude of an acoustic signal coming from the left vertical wall

Fig. 11 Amplitude of an acoustic signal coming from the left vertical wall.

First eigen state

Fig. 12 First eigen state (\(\lambda=(k/c)^2=14.695\)) close to \(15\) of eigenvalue problem: \(-\Delta \varphi = \lambda\varphi\) and \(\frac{\partial \varphi}{\partial\boldsymbol{n}} = 0\) on \(\Gamma\)}

Acoustics

Table of content