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 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
 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 Fig. 11. But when \(kc2\) is an eigenvalue of the problem, then the solution is not unique:

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");
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=19.4256\)) close to \(20\) of eigenvalue problem: \(-\Delta \varphi = \lambda\varphi\) and \(\frac{\partial \varphi}{\partial n} = 0\) on \(\Gamma\)}

Acoustics

Table of content