Modes of the Circular Acoustic Cavity

Solving the differential equation

The modes are found by solving the Helmholtz equation, same as for the rectangular cavity: 

	(∇^2)p((Overscript[x, →])) + (k^2)p((Overscript[x, →])) == 0

where p((Overscript[x, →])) is the pressure at position (Overscript[x, →]) and k is the wave number of the mode.

The symmetry of the cavity suggests that this equation be solved using cylindrical (polar) coordinates {r, θ}.

Needs["VectorAnalysis`"]
SetCoordinates[Cylindrical[r,θ,z]];

The Laplacian operator (∇^2) separates in this coordinate system. Assume a solution of the form:

p[r_,θ_,z_]:=pr[r]pθ[θ]

Then:  (r^2)(((∇^2)pr(r)pθ(θ))/(pr(r)pθ(θ)))+(r^2)(k^2)==0

Apart[(r^2)Simplify[Laplacian[p[r,θ,z]]/p[r,θ,z]]]+(r^2)(k^2)== 0

The first two terms form a function of r only; the last term is a function of θ only. Since the sum vanishes, each must equal a constant, ±(m^2):

PolarEquation=((pθ)''[θ]+(m^2) pθ[θ]==0);
RadialEquation=(((k^2) (r^2)-(m^2))pr[r]+r (pr'[r]+r pr''[r])==0);

Solve the polar equation first:

PolarSolution=DSolve[{PolarEquation},pθ[θ],θ]

Let's use the boundary conditions pθ(0)==pθ(2π)==1, so that the speaker is positioned at θ==0. Then we only need the cosine term, and m must be an integer:

pθ[m_,θ_]:=Evaluate[pθ[θ]/.(PolarSolution[[1]]/.{C[1]->1,C[2]->0})]

The radial equation is known as Bessel's equation. Its solutions are the Bessel functions of the first and second kind, (Subscript[J, m]) (x) and  (Subscript[Y, m]) (x):

RadialSolution=DSolve[RadialEquation,pr[r],r]

The functions (Subscript[Y, m]) (x) diverge at the origin x=0, so we use only the (Subscript[J, m]) (x):

pr[m_,r_]:=Evaluate[pr[r]/.(RadialSolution[[1]]/.{C[1]->1,C[2]->0})]

The solution so far is:

pr[m,r]pθ[m,θ]

The final boundary condition is that the cylindrical wall at r=1/2 is a free boundary, so
	 (d/dr)p(r,θ,z)(Subscript[|, r=1/2])==0. 
(We have chosen the diameter of the cavity as our unit of length) This boundary condition determines the wave numbers k, Because it requires:
	(Subscript[J, m]) '(k/2)==0.

Here are plots of (Subscript[J, m]) (x) for m=0…5:

Plot[Evaluate@Table[Tooltip[BesselJ[i,x],(Subscript[J, i])[x]],{i,0,5}],{x,0,15},PlotStyle->ColorData[41,"ColorList"]]

Only (Subscript[J, 0]) (x) doesn't vanish at the origin; for the others, (Subscript[J, m>0]) (0)==0, and the first m-1 derivatives vanish at the origin as well. Far from the origin the functions look like sinusoids with amplitudes that decrease as (x^(-1/2)).

Mathematica provides a function which finds the zeros of a Bessel function; here is a custom function which uses Mathematica's function to numerically finds zeros of the derivative of a Bessel function:

Find an initial estimate of the (n^th) zero of the derivative of (Subscript[J, m]) (x) (not counting the origin)

JPrimeZeroInit[0,n_]:=.5(BesselJZero[0,n]+BesselJZero[0,n+1])
JPrimeZeroInit[m_,1]:=m/(m+1.)BesselJZero[m,1]
JPrimeZeroInit[m_,n_]:=.5(BesselJZero[m,n-1]+BesselJZero[m,n])

Use the initial estimate of the (n^th) zero of the derivative to numerically find that zero

JPrimeZero[m_,n_]:=Block[{x},
x/.FindRoot[D[BesselJ[m,x],x]==0,{x,JPrimeZeroInit[m,n]}]
]

Here's a simple function which finds the coordinates {x,(Subscript[J, m])(x)} of the (n^th) zero of the derivative of (Subscript[J, m])(x):

point[m_,n_]:=Block[{x=JPrimeZero[m,n]},{x,BesselJ[m,x]}]

Show[
ListPlot[Evaluate@Table[point[m,1],{m,0,8}]],
ListPlot[Evaluate@Table[point[m,2],{m,0,5}]],
ListPlot[Evaluate@Table[point[m,3],{m,0,2}]],
Plot[Evaluate@Table[Tooltip[BesselJ[i,x],(Subscript[J, i])[x]],{i,0,8}],{x,0,11},PlotStyle->ColorData[41,"ColorList"]],
PlotLabel->Style["Zeros of the derivatives of some Bessel Functions","Subsubtitle"]
]

Determining the normal modes

Using the diameter of the cavity as our unit of length, then the (n^th) zero of the derivative
of (Subscript[J, m]) (x=k r) determines a normal mode wave number (Subscript[k, m n])=2x. 

Recall that the angular variation for the mode goes as cos ( m θ ), so the first mode index  m determines the number of nodal lines along diameters of the cavity for that mode.

The second nodal index n determines the number of circular nodal lines within the cavity, since it determines how many times (Subscript[J, m]) (x) passes through 0 before reaching the (n^th) zero of its derivative. In the case of (Subscript[J, 0]) (x), n is the number of circular nodal lines; for all other Bessel functions the number of circular nodal lines is n-1.

k[m_,n_]:=2 JPrimeZero[m,n]

Generate a list of the 15 modes with the lowest frequencies:

modes=Take[Sort[Flatten[Table[{k[m,n],{m,n}},{m,0,8},{n,1,3}],1]],15];

ListPlot[
Evaluate@Table[
Tooltip[modes[[i,1]]/modes[[1,1]],{2π/modes[[i,1]],modes[[i,2]]}],{i,1,Length[modes]}],
Filling->Axis,AxesOrigin-> {0,0}, Ticks->{None,Automatic},
PlotLabel->
Style["15 lowest frequency modes\nwavelength and mode index in tooltip","Subsubtitle"],AxesLabel->{None,k/(Subscript[k, lowest])}
]

Plots of nodal patterns

Clear[p]

A pressure function using Cartesian coordinates:

p[m_,n_]:=
(p[m,n]=With[{K=k[m,n]},
Function[{x,y},pr[m,Norm[x+I y]] pθ[m,Arg[x+I y]]/.k->K]])
p[{m_,n_}]:=p[m,n]

Use the pressure function to plot the nodal pattern for a given pair of mode indices:

plot[{m_,n_}]:=Show[
Graphics[{Thick,Circle[{0,0},1/2]}],
ContourPlot[p[m,n][x,y]==0,{x,-1/2,1/2},{y,-1/2,1/2},RegionFunction->Function[{x,y},Norm[x+I y]<1/2],AspectRatio->Automatic,Frame->False,PlotPoints->40]
]

Display the nodal patterns of the lowest frequency modes

Print[Style["The lowest frequency modes (freqs relative to first mode):","Subsubtitle"]];
Do[
Print[];
Print[
Grid[{
Table[modes[[i,1]]/modes[[1,1]],{i,Min[3j+1,Length[modes]],Min[3j+3,Length[modes]]}],
Table[With[{i=i},Tooltip[plot[modes[[i,2]]],modes[[i,2]]]],{i,Min[3j+1,Length[modes]],Min[3j+3,Length[modes]]}]
}]
],
{j,0,Floor[Length[modes]/3]-1}]

Accidental near-degeneracy of a pair of modes

Note that the two modes { 1,2 } and { 4, 1 } are very nearly degenerate, that is, they have almost the same frequency (2.89 times the lowest resonant frequency); the difference in frequency is only about 1 part in 400. Since the Q of the cavity is probably no more than 80 at this frequency, we would not expect to be able to resolve the two resonances. Consequently, we will see a nodal pattern which results from the sum of the pressure patterns from the two modes:

Grid[{{Round[Mean[{modes[[5,1]],modes[[6,1]]}]/modes[[1,1]],.001]},
{Show[
Graphics[{Thick,Circle[{0,0},1/2]}],
ContourPlot[p[1,2][x,y]+p[4,1][x,y]==0,{x,-1/2,1/2},{y,-1/2,1/2},RegionFunction->Function[{x,y},Norm[x+I y]<1/2],AspectRatio->Automatic,Frame->False]
]}
}]

A similar degeneracy occurs for the two modes { 1,3 } and { 7, 1 }:

Grid[{{Round[Mean[{modes[[12,1]],modes[[13,1]]}]/modes[[1,1]],.001]},
{Show[
Graphics[{Thick,Circle[{0,0},1/2]}],
ContourPlot[p[1,3][x,y]+p[7,1][x,y]==0,{x,-1/2,1/2},{y,-1/2,1/2},RegionFunction->Function[{x,y},Norm[x+I y]<1/2],AspectRatio->Automatic,Frame->False]
]}
}]