3D Elastic waves in a glass

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.


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