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"]
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.
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}*)