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:

• 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. 11 Amplitude of an acoustic signal coming from the left vertical wall.

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