I'm exploring the use of Mathematica (or one of it's online cousins) for finite element analysis of a piezoelectric resonator. If successful, I would be able to input the 3D resonator geometry, with boundary conditions and anisotropic material properties, and get out the resonant frequency. (I've successfully done similar analyses before using dedicated FEM software, but I've never attempted it in Wolfram language.)
While that broader goal is a long way off, there's one hurdle I know I'd need to overcome first. In such a simulation, it's necessary not to only pull the lowest n eigenvalues (as NDEigensystem easily does), but to pull eigenvalues "near" a known target eigenvalue (or frequency). This is critical, since the meaningful resonances are often 5th overtones of a specific mode of vibration, and thus it is computationally infeasible to calculate all eigenmodes/values below this quite high target mode. Does anyone know a way (using NDEigensystem or otherwise) to pull higher eigenvalues than the lowest n values?
To add specifics let's discuss a simple case. I've modified the beam example in the Wolfram FEM Tutorial to work with NDEigensystem to show the eigenvalues of the beam problem (although I'm missing density, so there is more to be refined here):
Needs["NDSolve`FEM`"]
op = {Inactive[
Div][({{0, -((Y \[Nu])/(1 - \[Nu]^2))}, {-((Y (1 - \[Nu]))/(
2 (1 - \[Nu]^2))), 0}}.Inactive[Grad][
v[x, y], {x, y}]), {x, y}] +
Inactive[
Div][({{-(Y/(1 - \[Nu]^2)),
0}, {0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}}.Inactive[
Grad][u[x, y], {x, y}]), {x, y}],
Inactive[
Div][({{0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}, {-((
Y \[Nu])/(1 - \[Nu]^2)), 0}}.Inactive[Grad][
u[x, y], {x, y}]), {x, y}] +
Inactive[
Div][({{-((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2))),
0}, {0, -(Y/(1 - \[Nu]^2))}}.Inactive[Grad][
v[x, y], {x, y}]), {x, y}]} /. {Y -> 10^3, \[Nu] -> 33/100};
Subscript[\[CapitalGamma], D] =
DirichletCondition[{u[x, y] == 0., v[x, y] == 0.}, x == 0];
\[CapitalOmega] =
ImplicitRegion[ True, {{x, 0, 5}, {y, 0,
1}}];
RegionPlot[\[CapitalOmega], AspectRatio -> Automatic]
{lam2,uvfun2} =
NDEigensystem [{op,
Subscript[\[CapitalGamma], D]} , {u, v}, {x, y} \[Element] \[CapitalOmega],5]
whiche=5;
Sqrt[lam2[[whiche]]]/(2*Pi)
mesh = uvfun2[[whiche,1]]["ElementMesh"];
Show[{
mesh["Wireframe"[ "MeshElement" -> "BoundaryElements"]],
ElementMeshDeformation[mesh, uvfun2[[whiche]]][
"Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]
This gives a nice picture of the vibrating cantilevered beam (see attached). But how could one arrive at higher-level eigenmodes without calculating each lower one?
Attachments:
