There are many examples of eigenmodes computations for surfaces with Mathematica, such as:
- https://www.wolfram.com/mathematica/new-in-10/pdes-and-finite-elements/solve-a-wave-equation-in-2d.html,
- https://www.wolfram.com/language/11/differential-eigensystems/find-eigenvalues-that-lie-in-an-interval.html?product=mathematica
- https://www.wolfram.com/language/11/differential-eigensystems/compute-eigenfunctions-for-a-clamped-membrane.html?product=mathematica
- Eigenfunctions of the Laplacian on an arbitrary mesh
I wanted to compute the slightly more complicated case of eigenmodes of a volume, or, more precisely, a shell (i.e., a volume with one dimension much small than the two others).
The mesh is created with:
ir = ImplicitRegion[0.8 <= x^2 + y^2 <= 1 && 0 <= z <= 3, {x, y, z}];
dr = DiscretizeRegion[ir, {{-1, 1}, {-1, 1}, {0, 2}}];
The 3D wave equation for a linear elastic material is
$$(\lambda+\mu)\boldsymbol\nabla(\boldsymbol\nabla \cdot \mathbf{u})+\mu \Delta \mathbf{u} = \rho \ddot{\mathbf{u}}$$
where lambda
and mu
are the Lamé parameters. I implemented this PDE (without the RHS) as:
uVec = Through[{u, v, w}[x, y, z]];
vars = {x, y, z};
pde = Table[(lambda + mu)*
D[Sum[D[uVec[[j]], vars[[j]]], {j, 3}], vars[[i]]] +
mu*Sum[D[uVec[[i]], {vars[[j]], 2}], {j, 3}], {i, 3}]
so the goal is to compute the eigenfunctions of this system of PDEs. Let us say the bottom of the cylinder is clamped, that is $u(x,y,0)=0$ or:
dirichlet = DirichletCondition[# == 0, z == 0] & /@ uVec;
The following retunrs the five first eigenfunctions:
{vals, funs} = NDEigensystem[Join[pde,dirichlet], uVec, {x, y, z} \[Element] ir,5]
Now, I would like to vizualize to modal shapes (funs
). Each element of funs
is a list of three InterpolatingFunctions
. I tried ContourPlot3D
or VectorPlot3D
but no matter the PlotRange
, I get InterpolatingFunction :: dmval Input value lies outside the range of data
, which I don't undertand: the funs
seem well defined on the domain of interest ($[-1,1]^2\times [0,3]$).
So my question is, did I do something wrong, and how to plot the eigenmodes of the shell?
Edit Using user21's answer, I am now able to plot some mode shapes. However, that's really strange. user21's system of PDEs pde2
is the same as mine (except the sign):
pde2 = stressOperator[Y, nu] /. material
Activate@pde === -Expand@pde2
(* True *)
but the results are completely different...