# Violin with f-holes and FEM simulation of Helmholtz resonance in 3D

This question is about FEM simulation of violin vibration modes in 3D. There are several problem around. One of them is Helmholtz resonance. Air inside the violin body has own frequencies dependent on geometry of f-holes. Practically this is the main reason why violin has so pronounced sound. We start from the body geometry definition (this is my design with given volume and area of f-holes, but it taken from the real violin):

xy = {{3.805405405405406, 3.34954954954955}, {3.8252252252252257,
6.6990990990991}, {3.9441441441441443,
7.9081081081081095}, {4.47927927927928,
8.601801801801802}, {5.014414414414414,
8.264864864864865}, {4.816216216216216,
7.8882882882882885}, {4.895495495495496,
7.630630630630631}, {5.232432432432433,
7.432432432432433}, {5.47027027027027,
7.491891891891892}, {5.648648648648649,
7.8882882882882885}, {5.668468468468468,
8.046846846846847}, {5.56936936936937,
8.403603603603605}, {5.252252252252252,
8.681081081081082}, {4.855855855855856,
8.780180180180182}, {4.518918918918919,
8.8}, {3.9639639639639643,
8.522522522522523}, {3.567567567567568,
7.967567567567568}, {3.3693693693693696,
7.372972972972973}, {3.2306306306306305,
6.67927927927928}, {3.1513513513513516,
3.3693693693693696}, {3.1513513513513516,
2.655855855855856}, {2.9729729729729732,
1.783783783783784}, {2.8738738738738743,
1.4666666666666668}, {2.100900900900901,
0.7927927927927928}, {1.7243243243243245,
1.3081081081081083}, {2.021621621621622,
1.7639639639639642}, {2.0414414414414415,
2.0414414414414415}, {1.9621621621621623,
2.23963963963964}, {1.6648648648648652,
2.4378378378378383}, {1.4666666666666668,
2.5171171171171176}, {1.10990990990991,
2.338738738738739}, {0.891891891891892,
1.9423423423423425}, {0.9315315315315316,
1.4072072072072073}, {1.5657657657657658,
0.7927927927927928}, {2.081081081081081,
0.6342342342342343}, {2.5963963963963965,
0.7927927927927928}, {3.0918918918918923,
1.2090090090090093}, {3.5081081081081082,
1.902702702702703}, {3.706306306306306, 2.6954954954954955}};

reg1 = RegionUnion[Disk[{0, 19.5/2}, 19.5/2],
Disk[{0, 36 - 15.5/2}, 15.5/2],
Rectangle[{-10, 15}, {10, 25}]]; reg2 =
RegionDifference[reg1,
RegionUnion[Disk[{-10, 20}, 9.5/2], Disk[{10, 20}, 9.5/2]]];
c0 = {0, 36 - 15.5/2}; c1 = {7.03562, 25};
f[x_] := c0[[2]] + x (c1[[2]] - c0[[2]])/(c1[[1]] - c0[[1]]); r1 =
Norm[c1 - {10, f[10]}];
reg3 = RegionDifference[reg2, Disk[{10, f[10]}, r1]];
f1[x_] := c0[[2]] - x (c1[[2]] - c0[[2]])/(c1[[1]] - c0[[1]]);
reg4 = RegionDifference[reg3, Disk[{-10, f1[-10]}, r1]]; c10 = {0,
19.5/2}; c11 = {8.215838362577491, 15};
f11[x_] := c10[[2]] + x (c11[[2]] - c10[[2]])/(c11[[1]] - c10[[1]]);
r2 = Norm[c11 - {10, f11[10]}];
reg5 = RegionDifference[reg4, Disk[{10, f11[10]}, r2]];
f12[x_] := c10[[2]] - x (c11[[2]] - c10[[2]])/(c11[[1]] - c10[[1]]);
reg6 = RegionDifference[reg5, Disk[{-10, f12[-10]}, r2]]; p6 =
RegionPlot[reg6, AspectRatio -> Automatic];
fh[xf_, yf_] :=
RegionUnion[
Polygon[Table[{xy[[i, 1]] - xf, xy[[i, 2]] + yf}, {i,
Length[xy]}]],
Polygon[Table[{-xy[[i, 1]] + xf, xy[[i, 2]] + yf}, {i,
Length[xy]}]]];


General view of the body from the front and back side

    Show[p6, Graphics[fh[7, 12], AspectRatio -> Automatic]]
dz = 3.79; reg8 =
ImplicitRegion[Element[{x, y}, reg6] && 0 <= z <= dz, {x, y, z}];
mesh3d1 = DiscretizeRegion[reg8, {{-10, 10}, {0, 36}, {0, dz}}]


Next step is the computation of air modes with

ca = 34321(*T=20C*); L = -Laplacian[u[x, y, z], {x, y, z}];

{vals, funs} =
NDEigensystem[{L,
DirichletCondition[u[x, y, z] == 0,
Element[{x, y}, fh[7, 11.49]] && z == dz]}, u,
Element[{x, y, z}, mesh3d1], 15];


Finally we plot 5 modes in 2D and the main mode in 3D

{Table[DensityPlot[funs[[i]][x, y, dz/2], {x, -10, 10}, {y, 0, 36},
PlotRange -> All, PlotLabel -> ca Sqrt[vals[[i]]]/(2 Pi),
ColorFunction -> "Rainbow", AspectRatio -> Automatic], {i, 1,
5}],
DensityPlot3D[
funs[[1]][x, y, z], {x, -10, 10}, {y, 0, 36}, {z, 0, dz},
PlotRange -> All, PlotLabel -> ca Sqrt[vals[[1]]]/(2 Pi),
ColorFunction -> "Rainbow", AspectRatio -> Automatic,
PlotLegends -> Automatic, PlotPoints -> 100, BoxRatios -> Automatic,
OpacityFunction -> None, Boxed -> False]}


Therefore the first mode is close to "A4" (440 Hz) tone. But we expecting "C4" (261.626 Hz), or "C#4" (277.183 Hz). The main reason of this discrepancies could be the wood plate vibration from the back side. Thus we define mesh, parameters of the wood plate and modes as follows

dreg = DiscretizeRegion[reg6, {{-10, 10}, {0, 36}},
MaxCellMeasure -> .05]
Y = 10.8*10^9; nu = 31/100; rho = 500; h = .003; d =
10^4 Sqrt[Y h^2/(12 rho (1 - nu^2))];Ld2 = {Laplacian[-d u[x, y], {x, y}] +
v[x, y], -d Laplacian[v[x, y], {x, y}]};

{vals, funs} =
NDEigensystem[{Ld2, DirichletCondition[u[x, y] == 0, True]}, {u, v},
Element[{x, y}, dreg], 5];

Table[DensityPlot[Re[funs[[i, 1]][x, y]], {x, y} \[Element] dreg,
PlotRange -> All, PlotLabel -> vals[[i]]/(2 Pi),
ColorFunction -> "Rainbow", AspectRatio -> Automatic], {i, 2,
Length[vals]}]


Ok! We have mode of 259.394 Hz and it is close to C4. The question is how we can connect this mode with Helmholtz resonance?

• How have you modelled the acoustic boundary conditions for the f-holes? I would imagine they could behave as Helmholtz resonators which could modify the acoustic natural frequencies substantially. The problem of the structural-acoustic interaction is non trivial. I would be interested in seeing a simpler problem of a cubical box with one side a diaphragm. This will require solving a coupled problem with complicated boundary conditions. One boundary condition is a common normal velocity, the other involves pressure. Pressure is difficult.
– Hugh
Dec 30, 2020 at 17:43
• @Hugh To compute modes we use DirichletCondition[u[x, y, z] == 0, Element[{x, y}, fh[7, 11.49]] && z == dz] on the surface of f-holes. Dec 30, 2020 at 17:46
• A slug of significant mass is made to oscillate in the f- holes; I think you are ignoring this. This is the phenomenon of end correction Can we calculate end corrections? by how much does this lower the frequency?
– Hugh
Dec 30, 2020 at 18:06
• @Hugh You are thinking about 0-D models, while we use 3D FEM. Dec 30, 2020 at 18:21
• In 3D there is still a slug of air that is driven in and out of the hole. It is simpler in pipes but occurs in all holes whatever their shape. This slug has inertia which you are not including. What is not clear is the importance of this slug compared with the flexibility of the walls.
– Hugh
Dec 30, 2020 at 18:57

We can solve the problem of connection plate mode with Helmholtz resonance by applying DirichletCondition[] to the Helmholtz equation on the back plate with wooden plate eigenfunction. Code to calculate plate modes (we use dreg from the code shown above)

    Y = 10.8*10^9; nu = 31/100; rho = 500; h = .003; d =
10^4 Sqrt[
Y h^2/(12 rho (1 - nu^2))]; Ld2 = {Laplacian[-d u[x, y], {x, y}] +
v[x, y], -d Laplacian[v[x, y], {x, y}]};
{vals, funs} =
NDEigensystem[{Ld2, DirichletCondition[u[x, y] == 0, True]}, {u, v},
Element[{x, y}, dreg], 10,
Method -> {"Interpolation" -> {"ExtrapolationHandler" -> \
{Automatic, "WarningMessage" -> False}}}];
Table[DensityPlot[Re[funs[[i, 1]][x, y]], {x, y} \[Element] dreg,
PlotRange -> All, PlotLabel -> vals[[i]]/(2 Pi),
ColorFunction -> "Rainbow", AspectRatio -> Automatic], {i, 2,
Length[vals]}]


Code to compute modes in the body with f-holes (we use mesh1=mesh3d1 and ca from the code above)

Do[solp[i] =
NDSolveValue[{-u[x, y, z] (vals[[i]]/ca)^2 -
Laplacian[u[x, y, z], {x, y, z}] == 0,
DirichletCondition[u[x, y, z] == Re[funs[[i, 1]][x, y]], z == 0],
DirichletCondition[u[x, y, z] == 0,
Element[{x, y}, fh[7, 11.49]] && z == dz]}, u,
Element[{x, y, z}, mesh1]];, {i, 2, Length[vals]}]

Table[Show[
ContourPlot[solp[i][x, y, dz - .1], {x, -10, 10}, {y, 0, 36},
PlotRange -> All, ColorFunction -> "Rainbow",
AspectRatio -> Automatic, PlotLegends -> Automatic,
PlotPoints -> 50, Contours -> 20, PlotLabel -> vals[[i]]/(2 Pi)],
Graphics[{Green,
Polygon[Table[{xy[[i, 1]] - 7, xy[[i, 2]] + 11.49}, {i,
Length[xy]}]],
Polygon[Table[{-xy[[i, 1]] + 7, xy[[i, 2]] + 11.49}, {i,
Length[xy]}]]}]], {i, 2, Length[vals]}]
`

• What a great work Alexander! Apr 5, 2022 at 7:35