3D Elastic waves in a glass

Question

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

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]

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

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

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

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]

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. — user21, Commented May 8, 2020 at 7:48
@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. — Alex Trounev, Commented May 8, 2020 at 14:24
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, — user21, Commented May 11, 2020 at 5:43
@user21 Are you talking about last code with "OpenCascade Link"` using? — Alex Trounev, Commented May 11, 2020 at 11:19

user21 · Accepted Answer · 2020-03-26 05:58:15Z

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

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

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]

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]

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

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.

Thank you for the mesh (+1) and links. They also use a thick cylinder. I will try to create a thin glass manually. — Alex Trounev, Commented Feb 7, 2020 at 10:01
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? — Hugh, Commented Feb 7, 2020 at 16:21
@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. — user21, Commented Mar 24, 2020 at 9:20
@user21 thank you. This is a very interesting update in version 12.1. I'll check how it works. — Alex Trounev, Commented Mar 24, 2020 at 10:48

Pinti · Accepted Answer · 2020-02-07 14:39:15Z

11

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

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

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}*)

edited Feb 7, 2020 at 14:39

answered Feb 7, 2020 at 13:26

Pinti

6,5231 gold badge18 silver badges48 bronze badges

$\begingroup$ (+1) nice. How many elements does the final mesh have? And I assume it is first order? $\endgroup$
– user21
Commented Feb 7, 2020 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
Commented Feb 7, 2020 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
Commented Feb 7, 2020 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
Commented Feb 7, 2020 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
Commented Feb 10, 2020 at 18:09

| Show 1 more comment

Stack Exchange Network

3D Elastic waves in a glass

2 Answers 2

Your Answer

Not the answer you're looking for? Browse other questions tagged
equation-solving
finite-element-method
or ask your own question.

Linked

Hot Network Questions

3D Elastic waves in a glass

2 Answers 2

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged equation-solvingfinite-element-method or ask your own question.

Linked

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
equation-solving
finite-element-method
or ask your own question.