22
$\begingroup$

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}}]

Figure 1

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]}

Figure 2

Figure 3

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]}]

Figure 4

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?

Improve this question
$\endgroup$
8
  • $\begingroup$ 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. $\endgroup$
    – Hugh
    Dec 30, 2020 at 17:43
  • $\begingroup$ @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. $\endgroup$
    – Alex Trounev
    Dec 30, 2020 at 17:46
  • $\begingroup$ 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? $\endgroup$
    – Hugh
    Dec 30, 2020 at 18:06
  • $\begingroup$ @Hugh You are thinking about 0-D models, while we use 3D FEM. $\endgroup$
    – Alex Trounev
    Dec 30, 2020 at 18:21
  • $\begingroup$ 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. $\endgroup$
    – Hugh
    Dec 30, 2020 at 18:57

1 Answer 1

Reset to default
7
$\begingroup$

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]}]

Figure 1 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]}]

Figure 2

Improve this answer
$\endgroup$
1
  • $\begingroup$ What a great work Alexander! $\endgroup$
    – Vangsnes
    Apr 5, 2022 at 7:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.