23
$\begingroup$

Take an empty glass, hit the side, the glass will make a sound that can be recorded using

s0=AudioCapture["C:\\Users\\...\\Desktop\\\\glass0.wav", MaxDuration -> 2]

Find the sound spectrum

Spectrogram[s0]

The photo shows a glass and a spectrum of sound Figure 1

Now we measure the dimensions of the glass, take the density, Young's modulus, glass Poisson's ratio from the reference book, compose the equations and find the eigenvalues

<< NDSolve`FEM`;
L = .14; L1 = .01; r1 = .085/2; r2 = .055/
  2; del = .006;(*cg=3962 m/s, 3980, 5100, 5640*);
reg = RegionUnion[
   ImplicitRegion[(r2 + (r1 - r2) (z - L1)/(L - L1))^2 <= 
      x^2 + y^2 <= (r2 + (r1 - r2) (z - L1)/(L - L1) + del)^2 && 
     L1 <= z <= L, {x, y, z}], 
   ImplicitRegion[
    0 <= x^2 + y^2 <= (r2 + del)^2 && 0 <= z <= L1, {x, y, z}]];
param = {Y -> 56*10^9, ν -> 25/100}; rho = 2500;
ClearAll[stressOperator];
stressOperator[
   Y_, ν_] := {Inactive[
      Div][{{0, 0, -((Y*ν)/((1 - 2*ν)*(1 + ν)))}, {0, 0, 
        0}, {-Y/(2*(1 + ν)), 0, 0}}.Inactive[Grad][
       w[t, x, y, z], {x, y, z}], {x, y, z}] + 
    Inactive[
      Div][{{0, -((Y*ν)/((1 - 2*ν)*(1 + ν))), 
        0}, {-Y/(2*(1 + ν)), 0, 0}, {0, 0, 0}}.Inactive[Grad][
       v[t, x, y, z], {x, y, z}], {x, y, z}] + 
    Inactive[
      Div][{{-((Y*(1 - ν))/((1 - 2*ν)*(1 + ν))), 0, 
        0}, {0, -Y/(2*(1 + ν)), 0}, {0, 
        0, -Y/(2*(1 + ν))}}.Inactive[Grad][
       u[t, x, y, z], {x, y, z}], {x, y, z}], 
   Inactive[
      Div][{{0, 0, 0}, {0, 
        0, -((Y*ν)/((1 - 
               2*ν)*(1 + ν)))}, {0, -Y/(2*(1 + ν)), 
        0}}.Inactive[Grad][w[t, x, y, z], {x, y, z}], {x, y, z}] + 
    Inactive[
      Div][{{0, -Y/(2*(1 + ν)), 
        0}, {-((Y*ν)/((1 - 2*ν)*(1 + ν))), 0, 0}, {0, 0, 
        0}}.Inactive[Grad][u[t, x, y, z], {x, y, z}], {x, y, z}] + 
    Inactive[
      Div][{{-Y/(2*(1 + ν)), 0, 
        0}, {0, -((Y*(1 - ν))/((1 - 2*ν)*(1 + ν))), 
        0}, {0, 0, -Y/(2*(1 + ν))}}.Inactive[Grad][
       v[t, x, y, z], {x, y, z}], {x, y, z}], 
   Inactive[
      Div][{{0, 0, 0}, {0, 
        0, -Y/(2*(1 + ν))}, {0, -((Y*ν)/((1 - 
               2*ν)*(1 + ν))), 0}}.Inactive[Grad][
       v[t, x, y, z], {x, y, z}], {x, y, z}] + 
    Inactive[
      Div][{{0, 0, -Y/(2*(1 + ν))}, {0, 0, 
        0}, {-((Y*ν)/((1 - 2*ν)*(1 + ν))), 0, 
        0}}.Inactive[Grad][u[t, x, y, z], {x, y, z}], {x, y, z}] + 
    Inactive[
      Div][{{-Y/(2*(1 + ν)), 0, 0}, {0, -Y/(2*(1 + ν)), 
        0}, {0, 0, -((Y*(1 - ν))/((1 - 
               2*ν)*(1 + ν)))}}.Inactive[Grad][
       w[t, x, y, z], {x, y, z}], {x, y, z}]};

{vals, funs} = 
 NDEigensystem[
  stressOperator[56*10^9, 1/4] + 
    rho {D[u[t, x, y, z], {t, 2}], D[v[t, x, y, z], {t, 2}], 
      D[w[t, x, y, z], {t, 2}]} == {0, 0, 0}, {u, v, w}, 
  t, {x, y, z} ∈ reg, 15];

Frequencies in Hertz

Abs[vals ]/(2 Pi)

Out[9]= {0.000389602, 0.000865814, 0.000865814, 0.000921462, \
0.000921462, 0.00136215, 0.00136215, 0.00152256, 0.00152256, \
0.0015598, 0.0015598, 2140.67, 2140.67, 2144.36, 2144.36}

And so we see that frequencies 2140-2144 explain the result of our experiment (in the spectrogram, the peak is about 2000 H). Build 3D functions u,v,w for frequency 2144.36

DensityPlot3D[Re[funs[[15, 1]][x, y, z]], {x, y, z} ∈ reg, 
 ColorFunction -> "Rainbow", OpacityFunction -> None, Boxed -> False, 
 PlotLabel -> Row[{"f = ", Abs[vals [[15]]]/2/Pi}], 
 BoxRatios -> Automatic, PlotPoints -> 50]

DensityPlot3D[Re[funs[[15, 2]][x, y, z]], {x, y, z} ∈ reg, 
 ColorFunction -> "Rainbow", OpacityFunction -> None, Boxed -> False, 
 PlotLabel -> Row[{"f = ", Abs[vals [[15]]]/2/Pi}], 
 BoxRatios -> Automatic, PlotPoints -> 50]
DensityPlot3D[Re[funs[[15, 3]][x, y, z]], {x, y, z} ∈ reg, 
 ColorFunction -> "Rainbow", Boxed -> False, 
 PlotLabel -> Row[{"f = ", Abs[vals [[15]]]/2/Pi}], 
 BoxRatios -> Automatic, PlotPoints -> 50]

Fig2

OK! Problems arise if we put del=0.003 (real glass wall thickness). First, the desired frequencies 2140-2144H disappear. Secondly, the 3D functions u,v,w look as if there are holes in the glass Figure3

Is it possible to get the desired result for del=.003?

Update 1. We use the algorithm proposed by user21 with a small modification and with the boundary condition DirichletCondition[{u[t, x, y, z] == 0, v[t, x, y, z] == 0, w[t, x, y, z] == 0}, z == 0]. Then the first 5 modes are consistent with the experiment (15 modes can be calculated with an error):

<< NDSolve`FEM`;
L = 0.14; L1 = 0.01; r1 = 0.085/2; r2 = 0.055/2; del = 0.003;


reg = RegionUnion[
   ImplicitRegion[(r2 + (r1 - r2) (z - L1)/(L - L1))^2 <= 
      x^2 + y^2 <= (r2 + (r1 - r2) (z - L1)/(L - L1) + del)^2 && 
     L1 <= z <= L, {x, y, z}], 
   ImplicitRegion[
    0 <= x^2 + y^2 <= (r2 + del)^2 && 0 <= z <= L1, {x, y, z}]];
(mesh = ToElementMesh[reg, 
    "BoundaryMeshGenerator" -> {"BoundaryDiscretizeRegion", 
      Method -> {"MarchingCubes", PlotPoints -> 31}}, 
    "MeshOrder" -> 1])["Wireframe"]

Modes

{vals, funs} = 
 NDEigensystem[{stressOperator[56*10^9, 1/4] + 
     rho {D[u[t, x, y, z], {t, 2}], D[v[t, x, y, z], {t, 2}], 
       D[w[t, x, y, z], {t, 2}]} == {0, 0, 0}, 
   DirichletCondition[{u[t, x, y, z] == 0, v[t, x, y, z] == 0, 
     w[t, x, y, z] == 0}, z == 0]}, {u, v, w}, 
  t, {x, y, z} \[Element] mesh, 5];

Modes in Hz

Abs[vals]/(2 Pi)

Out[]= {2047.63, 2048.03, 2048.03, 2336.35, 2336.35}

There are radial and azimuthal modes Figure4

Update 2. We use the algorithm proposed by Pinti with a modification and with the boundary condition DirichletCondition[{u[t, x, y, z] == 0, v[t, x, y, z] == 0, w[t, x, y, z] == 0}, y == 0]. Then the first 9 modes are consistent with the experiment (modes can be calculated without an error):

Get["MeshTools`"]

L = 0.14; L1 = 0.01; r1 = 0.085/2; r2 = 0.055/2; del = 0.003;

n1 = 5;
n2 = 31;
n3 = 5;
n4 = 12;
mesh2D = StructuredMesh[{{{r2, 0}, {r1, L}}, {{r2 - del, 
     0}, {r1 - del, L}}}, {n2, n1}]

mesh2D["Wireframe"[Axes -> True, AxesOrigin -> {0, 0}]]

Modes

{vals, funs} = 
 NDEigensystem[{stressOperator[56*10^9, 1/4] + 
     rho {D[u[t, x, y, z], {t, 2}], D[v[t, x, y, z], {t, 2}], 
       D[w[t, x, y, z], {t, 2}]} == {0, 0, 0}, 
   DirichletCondition[{u[t, x, y, z] == 0, v[t, x, y, z] == 0, 
     w[t, x, y, z] == 0}, y == 0]}, {u, v, w}, 
  t, {x, y, z} \[Element] mesh, 9];

vals in Hz

     Abs[vals]/(2 Pi)

Out[]= {23.1411, 1806.36, 1806.36, 1806.36, 1806.36, 1970.47, \
1970.47, 1970.58, 1970.58}

There are radial and azimuthal modes too Figure 5

Update 3. We use the algorithm proposed by user21 for version 12.1 with a small modification

<< NDSolve`FEM`;
L = 0.14; L1 = 0.01; del = 0.003; r1 = 0.085/2; r2 = 0.055/2;

polygon = 
  Polygon[{{0, 0, 0}, {r2 + del, 0, 0}, {r2 + del, 0, L1}, {r1 + del, 
     0, L}, {r1, 0, L}, {r2, 0, L1}, {0, 0, L1}}];

Needs["OpenCascadeLink`"]
shape = OpenCascadeShape[polygon];
axis = {{0, 0, 0}, {0, 0, 3/2 L}}; sweep = 
 OpenCascadeShapeRotationalSweep[shape, axis, 2 Pi];
bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[sweep, 
   "ShapeSurfaceMeshOptions" -> {"LinearDeflection" -> 0.0003}];


mesh = ToElementMesh[bmesh, AccuracyGoal -> 5, PrecisionGoal -> 5, 
  "MeshOrder" -> 1];


param = {Y -> 56*10^9, \[Nu] -> 25/100}; rho = 2500; cg = 
 Sqrt[56.*10^9/rho]; 


ClearAll[stressOperator];
stressOperator[
   Y_, \[Nu]_] := {Inactive[
      Div][{{0, 0, -((Y*\[Nu])/((1 - 2*\[Nu])*(1 + \[Nu])))}, {0, 0, 
        0}, {-Y/(2*(1 + \[Nu])), 0, 0}}.Inactive[Grad][
       w[t, x, y, z], {x, y, z}], {x, y, z}] + 
    Inactive[
      Div][{{0, -((Y*\[Nu])/((1 - 2*\[Nu])*(1 + \[Nu]))), 
        0}, {-Y/(2*(1 + \[Nu])), 0, 0}, {0, 0, 0}}.Inactive[Grad][
       v[t, x, y, z], {x, y, z}], {x, y, z}] + 
    Inactive[
      Div][{{-((Y*(1 - \[Nu]))/((1 - 2*\[Nu])*(1 + \[Nu]))), 0, 
        0}, {0, -Y/(2*(1 + \[Nu])), 0}, {0, 
        0, -Y/(2*(1 + \[Nu]))}}.Inactive[Grad][
       u[t, x, y, z], {x, y, z}], {x, y, z}], 
   Inactive[
      Div][{{0, 0, 0}, {0, 
        0, -((Y*\[Nu])/((1 - 
               2*\[Nu])*(1 + \[Nu])))}, {0, -Y/(2*(1 + \[Nu])), 
        0}}.Inactive[Grad][w[t, x, y, z], {x, y, z}], {x, y, z}] + 
    Inactive[
      Div][{{0, -Y/(2*(1 + \[Nu])), 
        0}, {-((Y*\[Nu])/((1 - 2*\[Nu])*(1 + \[Nu]))), 0, 0}, {0, 0, 
        0}}.Inactive[Grad][u[t, x, y, z], {x, y, z}], {x, y, z}] + 
    Inactive[
      Div][{{-Y/(2*(1 + \[Nu])), 0, 
        0}, {0, -((Y*(1 - \[Nu]))/((1 - 2*\[Nu])*(1 + \[Nu]))), 
        0}, {0, 0, -Y/(2*(1 + \[Nu]))}}.Inactive[Grad][
       v[t, x, y, z], {x, y, z}], {x, y, z}], 
   Inactive[
      Div][{{0, 0, 0}, {0, 
        0, -Y/(2*(1 + \[Nu]))}, {0, -((Y*\[Nu])/((1 - 
               2*\[Nu])*(1 + \[Nu]))), 0}}.Inactive[Grad][
       v[t, x, y, z], {x, y, z}], {x, y, z}] + 
    Inactive[
      Div][{{0, 0, -Y/(2*(1 + \[Nu]))}, {0, 0, 
        0}, {-((Y*\[Nu])/((1 - 2*\[Nu])*(1 + \[Nu]))), 0, 
        0}}.Inactive[Grad][u[t, x, y, z], {x, y, z}], {x, y, z}] + 
    Inactive[
      Div][{{-Y/(2*(1 + \[Nu])), 0, 0}, {0, -Y/(2*(1 + \[Nu])), 
        0}, {0, 0, -((Y*(1 - \[Nu]))/((1 - 
               2*\[Nu])*(1 + \[Nu])))}}.Inactive[Grad][
       w[t, x, y, z], {x, y, z}], {x, y, z}]};

{vals, funs} = 
 NDEigensystem[{stressOperator[56*10^9, 1/4] + 
     rho {D[u[t, x, y, z], {t, 2}], D[v[t, x, y, z], {t, 2}], 
       D[w[t, x, y, z], {t, 2}]} == {0, 0, 0}, 
   DirichletCondition[{u[t, x, y, z] == 0, v[t, x, y, z] == 0, 
     w[t, x, y, z] == 0}, z == 0]}, {u, v, w}, 
  t, {x, y, z} \[Element] mesh, 12];

vals in Hz

     Abs[vals]/(2 Pi)

{1973.97, 1973.97, 1974.86, 1974.86, 2169.47, 2169.47, 2250.23, 2250.23, 4183.69, 4183.69, 5532.12, 5532.12} Visualisation of 3 modes

DensityPlot3D[Re[funs[[1, 1]][x, y, z]], {x, y, z} \[Element] mesh, 
 ColorFunction -> "Rainbow", OpacityFunction -> None, Boxed -> False, 
 PlotLabel -> Row[{"f = ", Abs[vals [[1]]]/2/Pi}], 
 BoxRatios -> Automatic, PlotPoints -> 50]
DensityPlot3D[Re[funs[[5, 1]][x, y, z]], {x, y, z} \[Element] mesh, 
 ColorFunction -> "Rainbow", OpacityFunction -> None, Boxed -> False, 
 PlotLabel -> Row[{"f = ", Abs[vals [[5]]]/2/Pi}], 
 BoxRatios -> Automatic, PlotPoints -> 50]
DensityPlot3D[Re[funs[[7, 1]][x, y, z]], {x, y, z} \[Element] mesh, 
 ColorFunction -> "Rainbow", OpacityFunction -> None, Boxed -> False, 
 PlotLabel -> Row[{"f = ", Abs[vals [[7]]]/2/Pi}], 
 BoxRatios -> Automatic, PlotPoints -> 50]

Figure 7

$\endgroup$
  • $\begingroup$ I was wondering if you would like to contribute this model to the FEMAddOns. There is a section for Application examples. Might be a nice contribution. Just a thought. $\endgroup$ – user21 May 8 at 7:48
  • $\begingroup$ @user21 Thank you for your proposal. I have tried several times to contribute something on github with .nb extension, but never successful. I don't understand how it is working. If you can, please, do it. $\endgroup$ – Alex Trounev May 8 at 14:24
  • $\begingroup$ We could do it together if you want. If you could send me an initial draft of your work notebook we can take it from there, $\endgroup$ – user21 May 11 at 5:43
  • $\begingroup$ @user21 Are you talking about last code with "OpenCascade Link"` using? $\endgroup$ – Alex Trounev May 11 at 11:19
  • $\begingroup$ Whatever you think is appropriate to get a good result. $\endgroup$ – user21 May 12 at 4:58
10
$\begingroup$

You get a better mesh with a different boundary mesh generator:

(mesh = ToElementMesh[reg, 
    "BoundaryMeshGenerator" -> \
{"BoundaryDiscretizeRegion",
      Method -> {"MarchingCubes", PlotPoints -> 33}}, 
    "MeshOrder" -> 1,
    "MaxCellMeasure"\[Rule]0.000000005])["Wireframe"]

enter image description here

For that mesh I get

Abs[vals]/(2 Pi)
(*{0.000502385, 0.000502385, 0.00072869, 0.00072869, \
0.000733392, 0.000733392, 0.0010404, 0.0010404, 0.00150767, \
0.00150767, 0.00151325, 0.00151325, 0.308656, 2238.88, 2238.88}*)

And the 14th mode looks like:

MeshRegion[
 ElementMeshDeformation[mesh, Re[Through[funs[[14]]["ValuesOnGrid"]]],
   "ScalingFactor" -> 10^9]]

enter image description here

Two other comments: the fact that NDEigensystem gives messages suggests to me that this mesh is still not good enough; as you see I also used MeshOrder->1 as I did not want to wait for a second order mesh to finish. But you might want to try that and a finer mesh. Probably by using more plot points. Perhaps generate the boundary mesh manually?

A second thing that come to mind is that I think you should have some rigid body modes because the glass stands on the table. Maybe experiment with

DirichletCondition[{u[t, x, y, z] == 0, v[t, x, y, z] == 0, 
  w[t, x, y, z] == 0}, x == 0]

Also, there is a nice Bell Acoustics customer example in the FEMAddOns. You can install that with

ResourceFunction["FEMAddOnsInstall"][]

and find it on the Applications guide page

FEMAddOns/guide/FEMApplications

or have a look at the cloud version of that notebook.

Hope this helps.

Update: 12.1

Another way to generate the mesh is to make use of the OpenCascadeLink. For this we generate a flat cross section of the glass in 3D.

polygon = 
  Polygon[{{0, 0, 0}, {r2 + del, 0, 0}, {r2 + del, 0, L1}, {r1 + del, 
     0, L}, {r1, 0, L}, {r2, 0, L1}, {0, 0, L1}}];
Graphics3D[{FaceForm[], EdgeForm[Black], polygon}, Boxed -> False]

enter image description here

We load the link

Needs["OpenCascadeLink`"]

and convert the polygon into an OCCT shape:

shape = OpenCascadeShape[polygon];

We set up an axis of revolution and sweep the polygon.

axis = {{0, 0, 0}, {0, 0, 3/2 L}};
sweep = OpenCascadeShapeRotationalSweep[shape, axis, 2 \[Pi]];

Here is a visual of the result:

bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[sweep, 
   "ShapeSurfaceMeshOptions" -> {"LinearDeflection" -> 0.00125}];
Show[Graphics3D[{{Red, polygon}, {Blue, Thick, Arrow[axis]}}], 
 bmesh["Wireframe"], Boxed -> False]

enter image description here

You see the original polygon in red and the blue arrow is the rotational axis. From here we can generate the mesh in the same way:

mesh = ToElementMesh[bmesh, "MeshOrder" -> 1(*,
  "MaxCellMeasure"\[Rule]0.000000005*)]

mesh["Wireframe"[
  "MeshElementStyle" -> 
   Directive[Opacity[0.2], Specularity[White, 17], FaceForm[White], 
    EdgeForm[]]]]

enter image description here

This is a much better approximation of the geometry. Nevertheless finding the eigenvalues remains challenging as there is a strong dependency of the eigenvalues on the mesh.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you for the mesh (+1) and links. They also use a thick cylinder. I will try to create a thin glass manually. $\endgroup$ – Alex Trounev Feb 7 at 10:01
  • $\begingroup$ Finite element codes often have shell elements. It looks like they would be useful here. Is there any chance of getting such elements available in Mathematica? $\endgroup$ – Hugh Feb 7 at 16:21
  • $\begingroup$ @Hugh, sorry I saw the comment but then forgot to reply. Shell elements are currently not planed. There is still too much other stuff I can improve before going down this alley. $\endgroup$ – user21 Mar 24 at 9:20
  • $\begingroup$ @user21 thank you. This is a very interesting update in version 12.1. I'll check how it works. $\endgroup$ – Alex Trounev Mar 24 at 10:48
  • $\begingroup$ @user21 See update to my post. $\endgroup$ – Alex Trounev Mar 25 at 18:28
10
$\begingroup$

MeshTools package can help insituations where we need fine control of mesh density and shape.

First we define a 2D mesh for glass outline and revolve it around vertical axis. Then we merge it with cylinder mesh for glass bottom. We get 1st order mesh, but it can be converted to 2 order with MeshOrderAlteration from "NDSolve`FEM`" built-in package.

Get["MeshTools`"]

L = 0.14; L1 = 0.01; r1 = 0.085/2; r2 = 0.055/2; del = 0.003;

n1 = 2;
n2 = 40;
n3 = 5;
n4 = 12;
mesh2D = MergeMesh[{
   StructuredMesh[{{{r2, L1}, {r1, L}}, {{r2 - del, L1}, {r1 - del, L}}}, {n2, n1}],
   StructuredMesh[{{{r2, 0}, {r2, L1}}, {{r2 - del, 0}, {r2 - del,  L1}}}, {n3, n1}]
 }]

mesh2D["Wireframe"[Axes -> True, AxesOrigin -> {0, 0}]]

mesh2D

mesh = MergeMesh[{
   CylinderMesh[{{0, 0, 0}, {0, L1, 0}}, r2 - del, {n4, n1}],
   RevolveMesh[mesh2D, {0, 2 Pi}, 4*n4]
 }]
(* ElementMesh[{{-0.0425,0.0425},{0.,0.14},{-0.0425,0.0425}}, {HexahedronElement["<"4896">"]}]*)


mesh["Wireframe"["MeshElementStyle" -> FaceForm@LightBlue]]

mesh3D

For the calculated frequencies we get the following list.

Abs[vals]/(2 Pi)
(*{0.000290029, 0.000355687, 0.000355687, 0.000584401, 0.000584401, 0.000724522, 0.000724522, 0.000903912, 0.000903912, 0.000903912, 0.000903912, 1907.22, 1907.22, 1907.6, 1907.6}*)
| improve this answer | |
$\endgroup$
  • $\begingroup$ (+1) nice. How many elements does the final mesh have? And I assume it is first order? $\endgroup$ – user21 Feb 7 at 13:31
  • $\begingroup$ @Pint (+1) Looks good, but how did you find vals? My kernel is breaking on this mesh with using NDEigenvalues[] or NDEigensystem[] $\endgroup$ – Alex Trounev Feb 7 at 13:55
  • $\begingroup$ @AlexTrounev Before I have copied wrong set of mesh division settings. These current settings work for me on Windows, MMa 12. For many other combinations of element division kernel also dies on me. I don't know why. $\endgroup$ – Pinti Feb 7 at 14:41
  • $\begingroup$ @Pinti Unfortunately, this mesh also does not work in version 12 for Windows 10 64 bit. However, I found many other combinations of n1,n2,n3,n4 that work, but none of them give the desired frequencies. $\endgroup$ – Alex Trounev Feb 7 at 15:50
  • $\begingroup$ @Pinti I managed to find a combination of mesh parameters on which NDEigensystem[] finds the modes without errors - see update2. $\endgroup$ – Alex Trounev Feb 10 at 18:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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