Skip to content

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:

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{eqnarray} k^{2}v + c^{2}\Delta v &= 0 &\hbox{ in } \Omega\\ \frac{\p v}{\p n}|_\Gamma &= g. \end{eqnarray}

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
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
// Parameters
real kc2 = 1.;
func g = y*(1.-y);

// Mesh
border a0(t=0., 1.){x=5.; y=1.+2.*t;}
border a1(t=0., 1.){x=5.-2.*t; y=3.;}
border a2(t=0., 1.){x=3.-2.*t; y=3.-2.*t;}
border a3(t=0., 1.){x=1.-t; y=1.;}
border a4(t=0., 1.){x=0.; y=1.-t;}
border a5(t=0., 1.){x=t; y=0.;}
border a6(t=0., 1.){x=1.+4.*t; y=t;}

mesh Th = buildmesh(a0(20) + a1(20) + a2(20)
    + a3(20) + a4(20) + a5(20) + a6(20));

// Fespace
fespace Vh(Th, P1);
Vh u, v;

// Solve
solve sound(u, v)
    = int2d(Th)(
          u*v * kc2
        - dx(u)*dx(v)
        - dy(u)*dy(v)
    )
    - int1d(Th, a4)(
          g * v
    )
    ;

// Plot
plot(u, wait=1, ps="Sound.eps");

Results are on figure 1. 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
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
// Parameters
real sigma = 20; //value of the shift

// Problem
// OP = A - sigma B ; // The shifted matrix
varf op(u1, u2)
    = int2d(Th)(
          dx(u1)*dx(u2)
        + dy(u1)*dy(u2)
        - sigma* u1*u2
    )
    ;

varf b([u1], [u2])
    = int2d(Th)(
        u1*u2
    )
    ; // No Boundary condition see note \ref{note BC EV}

matrix OP = op(Vh, Vh, solver=Crout, factorize=1);
matrix B = b(Vh, Vh, solver=CG, eps=1e-20);

// Eigen values
int nev=2; // Number of requested eigenvalues near sigma

real[int] ev(nev);  // To store the nev eigenvalue
Vh[int] eV(nev);    // To store the nev eigenvector

int k=EigenValue(OP, B, sym=true, sigma=sigma, value=ev, vector=eV,
    tol=1e-10, maxit=0, ncv=0);

cout << ev(0) << " 2 eigen values " << ev(1) << endl;
v = eV[0];
plot(v, wait=true, ps="eigen.eps");
Fig. 1: Amplitude of an acoustic signal coming from the left vertical wall. First eigen state (\lambda=(k/c)^2=19.4256) close to 20 of eigenvalue problem : -\Delta \varphi = \lambda \varphi and \frac{\partial \varphi}{\partial n} = 0 on \Gamma}
Acoustics Acoustics