6
$\begingroup$

We are trying to use NDEigensystem to solve for the vibrational modes of a cantilever with triangular cross-section. However, many of the solutions provided by NDEigensystem have imaginary eigenfrequencies and are therefore unphysical. We have determined that the imaginary eigenfrequencies result from NDEigensystem passing a non-Hermitian system matrix to Eigensystem to solve. We do not know why NDEigensystem is doing this or how to correct it. Any help will be most appreciated.

The relevant PDE is

\begin{equation} \mu\nabla^2 \vec u + (\lambda + \mu) \nabla(\nabla\cdot \vec u) =-\rho \omega^2 \vec u, \end{equation}

where

\begin{equation} \vec u(x,y,z) \end{equation}

is a 3D vector function representing the vibrational displacement field, $\rho$ is the material density, and $\lambda$ and $\mu$ are elastic constants.

The geometry is depicted below.

TrianglePillar[W_, L_] = 
Prism[{{0, 
 W*Sqrt[3]/3, -L/2}, {-W*1/2, -W*Sqrt[3]/6, -L/2}, {W*1/2, -W*
  Sqrt[3]/6, -L/2}, {0, W*Sqrt[3]/3, L/2},
{-W*1/2, -W*Sqrt[3]/6, L/2}, {W*1/2, -W*Sqrt[3]/6, L/2}}];

enter image description here

The long axis of the cantilever extends from z=-1/2 to z=1/2. The cantilever has an equilateral triangular cross-section with a sidelength of 1/2.

Dirichlet boundary conditions apply at the clamped end of the cantilever (i.e. $\vec u(x,y,-1/2)=0$) and Neumann boundary conditions apply to all other boundaries.

Our Mathematica code to solve for the first 50 mode of the cantilever is:

LIxx = FullSimplify[( {
{λ + 2 μ, 0, 0},
{0, μ, 0},
{0, 0, μ}
} )];

LIxy = FullSimplify[( {
{0, λ, 0},
{μ, 0, 0},
{0, 0, 0}
} )];

LIxz = FullSimplify[( {
{0, 0, λ},
{0, 0, 0},
{μ, 0, 0}
} )];

LIyy = FullSimplify[( {
{μ, 0, 0},
{0, λ + 2 μ, 0},
{0, 0, μ}
} )];

LIyx = FullSimplify[( {
{0, λ, 0}, 
{μ, 0, 0},
{0, 0, 0}
} )];

LIyz = FullSimplify[( {
{0, 0, 0},
{0, 0, λ},
{0, μ, 0}
} )];

LIzz = FullSimplify[( {
{μ, 0, 0},    
{0, μ, 0},
{0, 0, λ + 2 μ}
} )];

LIzx = FullSimplify[( {
{0, 0, μ},
{0, 0, 0},
{λ, 0, 0}
} )];

LIzy = FullSimplify[( {
{0, 0, 0},
{0, 0, μ},
{0, λ, 0}
} )];

PhononEquation3D = {
Inactive[Div][(LIxx.Inactive[Grad][u[x, y, z], {x, y, z}]), {x, y, z}]+
Inactive[Div][(LIxy.Inactive[Grad][v[x, y, z], {x, y, z}]), {x, y, z}]+
Inactive[Div][(LIxz.Inactive[Grad][w[x, y, z], {x, y, z}]), {x, y, z}],

Inactive[Div][(LIyy.Inactive[Grad][v[x, y, z], {x, y, z}]), {x, y, z}]+
Inactive[Div][(LIyx.Inactive[Grad][u[x, y, z], {x, y, z}]), {x, y, z}]+
Inactive[Div][(LIyz.Inactive[Grad][w[x, y, z], {x, y, z}]), {x, y, z}],

Inactive[Div][(LIzz.Inactive[Grad][w[x, y, z], {x, y, z}]), {x, y, z}]+
Inactive[Div][(LIzx.Inactive[Grad][u[x, y, z], {x, y, z}]), {x, y, z}]+
Inactive[Div][(LIzy.Inactive[Grad][v[x, y, z], {x, y, z}]), {x, y, z}]} /. {λ -> 85, μ -> 536};

{ω, f} = NDEigensystem[{PhononEquation3D,
DirichletCondition[{u[x, y, z] == 0, v[x, y, z] == 0, w[x, y, z] ==0},z == -1/2]}, {u, v, w},Element[{x, y, z}, TrianglePillar[0.5, 1]],
50,
Method -> {"SpatialDiscretization" -> {"FiniteElement", {"MeshOptions" -> {"MaxCellMeasure" -> (1/100)}}}}];

The eigenvalues we obtain are:

ω
{-132.962 + 0. I, -134.418 + 0. I, -2835.36 + 0. I, -3417.48 +
0. I, -3667.57 + 0. I, -15998.9 + 0. I, -16938.7 + 0. I, -25987.1 +
0. I, -29751.5 + 0. I, -42989.6 - 4594.81 I, -42989.6 +
4594.81 I, -44516.6 + 0. I, -60178.9 + 0. I, -62091. +
0. I, -67288.5 - 1358.7 I, -67288.5 + 1358.7 I, -68955.9 +
0. I, -75728.8 + 0. I, -78441.7 + 0. I, -81429. + 0. I, -96662.2 +
0. I, -101594. + 0. I, -106433. - 2621.72 I, -106433. +
2621.72 I, -126487. + 0. I, -136803. + 0. I, -146696. +
0. I, -151486. + 0. I, -158087. + 0. I, -164934. + 0. I, -169051. +
0. I, -183941. + 0. I, -185686. + 0. I, -201331. + 0. I, -219622. + 
0. I, -236595. + 0. I, -239297. + 0. I, -251981. + 0. I, -265628. +
0. I, -274303. - 4060.39 I, -274303. + 4060.39 I, -287874. +
0. I, -306104. + 0. I, -320046. + 0. I, -345800. + 0. I, -355459. +
0. I, -373953. + 0. I, -380989. + 0. I, -393936. +
16685.2 I, -393936. - 16685.2 I}

which are a mixture of purely real (expected) and imaginary values (erroneous). In attempting to diagnose the cause of the imaginary values, we discovered that the system matrix passed to Eigensystem to solve is non-Hermitian. This is also the case if we remove the Dirichlet boundary condition.

Any ideas?

$\endgroup$

closed as off-topic by Jens, MarcoB, xzczd, user9660, Jason B. Jan 28 '16 at 8:27

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question arises due to a simple mistake such as a trivial syntax error, incorrect capitalization, spelling mistake, or other typographical error and is unlikely to help any future visitors, or else it is easily found in the documentation." – Jens, MarcoB, xzczd, Community, Jason B.
If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – Michael E2 Jan 27 '16 at 2:14
  • 1
    $\begingroup$ The imaginary eigenvalues turn real if I set Dirichlet conditions on the entire boundary. So I guess it's got something to do with the Neumann conditions. Looks like a bug, but needs further investigation... I was able to suppress the imaginary parts in the first approach (with your Dirichlet condition) by increasing "MaxCellMeasure" -> (1/800). But they don't go away completely... $\endgroup$ – Jens Jan 27 '16 at 4:22
  • $\begingroup$ Do you know WHY you obtain a non-hermitian system matrix? It's kind of weird, the mechanical problem is purely conservative and the problem is selbst-adjoint. Have you tried to analyse a normal beam with a square surface? Maybe for your triangle this is an artefact of the discretization, but I am not sure. $\endgroup$ – Mauricio Fernández Jan 27 '16 at 8:27
  • 1
    $\begingroup$ @MauricioLobos, while in this case things worked out just fine, I am not sure if a continuous, self-adjoined operator implies that it's (FEM) discretized version also has to be self-adjoined under all circumstances. $\endgroup$ – user21 Jan 28 '16 at 0:31
  • $\begingroup$ Thank you everyone for your help, it was most appreciated. User21 correctly identified our typo in the matrix LIyx and the code now works well. $\endgroup$ – user37263 Jan 28 '16 at 2:44
5
$\begingroup$

This due to a simple typo. If you look at

{{LIxx, LIxy, LIxz}, {LIyx, LIyy, LIyz}, {LIzx, LIzy, 
   LIzz}} // MatrixForm

enter image description here

You'll see that the LIyx entry is not quite right. If you change that to

LIxx = FullSimplify[({{λ + 2 μ, 0, 0}, {0, μ, 0}, {0, 0, μ}})];
LIxy = FullSimplify[({{0, λ, 0}, {μ, 0, 0}, {0, 0, 0}})];
LIxz = FullSimplify[({{0, 0, λ}, {0, 0, 0}, {μ, 0, 0}})];
LIyx = FullSimplify[({{0, μ, 0}, {λ, 0, 0}, {0, 0, 0}})];
LIyy = FullSimplify[({{μ, 0, 0}, {0, λ + 2 μ, 0}, {0, 0, μ}})];
LIyz = FullSimplify[({{0, 0, 0}, {0, 0, λ}, {0, μ, 0}})];
LIzx = FullSimplify[({{0, 0, μ}, {0, 0, 0}, {λ, 0, 0}})];
LIzy = FullSimplify[({{0, 0, 0}, {0, 0, μ}, {0, λ, 0}})];
LIzz = FullSimplify[({{μ, 0, 0}, {0, μ, 0}, {0, 0, λ + 2 μ}})];

{{LIxx, LIxy, LIxz}, {LIyx, LIyy, LIyz}, {LIzx, LIzy, 
   LIzz}} // MatrixForm

enter image description here

things work as you intended.

As a side note, you may want to use more than 6 elements to compute the eigenvalues.

$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.