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
 2 real kc2 = 1.;
 3 func g = y*(1.-y);
 4 
 5 // Mesh
 6 border a0(t=0., 1.){x=5.; y=1.+2.*t;}
 7 border a1(t=0., 1.){x=5.-2.*t; y=3.;}
 8 border a2(t=0., 1.){x=3.-2.*t; y=3.-2.*t;}
 9 border a3(t=0., 1.){x=1.-t; y=1.;}
10 border a4(t=0., 1.){x=0.; y=1.-t;}
11 border a5(t=0., 1.){x=t; y=0.;}
12 border a6(t=0., 1.){x=1.+4.*t; y=t;}
13 
14 mesh Th = buildmesh(a0(20) + a1(20) + a2(20)
15     + a3(20) + a4(20) + a5(20) + a6(20));
16 
17 // Fespace
18 fespace Vh(Th, P1);
19 Vh u, v;
20 
21 // Solve
22 solve 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
34 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 // Parameters
 2 real sigma = 20; //value of the shift
 3 
 4 // Problem
 5 // OP = A - sigma B ; // The shifted matrix
 6 varf op(u1, u2)
 7    = int2d(Th)(
 8         dx(u1)*dx(u2)
 9       + dy(u1)*dy(u2)
10       - sigma* u1*u2
11    )
12    ;
13 
14 varf b([u1], [u2])
15    = int2d(Th)(
16         u1*u2
17    )
18    ; // No Boundary condition see note \ref{note BC EV}
19 
20 matrix OP = op(Vh, Vh, solver=Crout, factorize=1);
21 matrix B = b(Vh, Vh, solver=CG, eps=1e-20);
22 
23 // Eigen values
24 int nev=2; // Number of requested eigenvalues near sigma
25 
26 real[int] ev(nev);  // To store the nev eigenvalue
27 Vh[int] eV(nev);    // To store the nev eigenvector
28 
29 int k=EigenValue(OP, B, sym=true, sigma=sigma, value=ev, vector=eV,
30    tol=1e-10, maxit=0, ncv=0);
31 
32 cout << ev(0) << " 2 eigen values " << ev(1) << endl;
33 v = eV[0];
34 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=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