2
$\begingroup$

I am trying to visualise the eigenvectors of a thick disk. I manage to calculate them and look at one but fail to look at more -probably due to using too much memory. I start by making a mesh and then use NDEigensystem.

Needs["NDSolve`FEM`"]
ClearAll[stressOperator];
stressOperator[
  Y_, \[Nu]_] := {Inactive[
     Div][{{0, 0, -((Y*\[Nu])/((1 - 2*\[Nu])*(1 + \[Nu])))}, {0, 0, 
       0}, {-Y/(2*(1 + \[Nu])), 0, 0}}.Inactive[Grad][
      w[x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{0, -((Y*\[Nu])/((1 - 2*\[Nu])*(1 + \[Nu]))), 
       0}, {-Y/(2*(1 + \[Nu])), 0, 0}, {0, 0, 0}}.Inactive[Grad][
      v[x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{-((Y*(1 - \[Nu]))/((1 - 2*\[Nu])*(1 + \[Nu]))), 0, 
       0}, {0, -Y/(2*(1 + \[Nu])), 0}, {0, 
       0, -Y/(2*(1 + \[Nu]))}}.Inactive[Grad][
      u[x, y, z], {x, y, z}], {x, y, z}], 
  Inactive[Div][{{0, 0, 0}, {0, 
       0, -((Y*\[Nu])/((1 - 
              2*\[Nu])*(1 + \[Nu])))}, {0, -Y/(2*(1 + \[Nu])), 
       0}}.Inactive[Grad][w[x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{0, -Y/(2*(1 + \[Nu])), 
       0}, {-((Y*\[Nu])/((1 - 2*\[Nu])*(1 + \[Nu]))), 0, 0}, {0, 0, 
       0}}.Inactive[Grad][u[x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{-Y/(2*(1 + \[Nu])), 0, 
       0}, {0, -((Y*(1 - \[Nu]))/((1 - 2*\[Nu])*(1 + \[Nu]))), 0}, {0,
        0, -Y/(2*(1 + \[Nu]))}}.Inactive[Grad][
      v[x, y, z], {x, y, z}], {x, y, z}], 
  Inactive[Div][{{0, 0, 0}, {0, 
       0, -Y/(2*(1 + \[Nu]))}, {0, -((Y*\[Nu])/((1 - 
              2*\[Nu])*(1 + \[Nu]))), 0}}.Inactive[Grad][
      v[x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{0, 0, -Y/(2*(1 + \[Nu]))}, {0, 0, 
       0}, {-((Y*\[Nu])/((1 - 2*\[Nu])*(1 + \[Nu]))), 0, 0}}.Inactive[
       Grad][u[x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{-Y/(2*(1 + \[Nu])), 0, 0}, {0, -Y/(2*(1 + \[Nu])), 0}, {0,
        0, -((Y*(1 - \[Nu]))/((1 - 2*\[Nu])*(1 + \[Nu])))}}.Inactive[
       Grad][w[x, y, z], {x, y, z}], {x, y, z}]}
r = 0.3; (* radius of cylinder *)
L = 0.2; (* length of cylinder *)
Y = 10^6 ;(* Ratio of modulus of elasticity to density *)
\[Nu] = 0.2; (* Poission ratio *)

mesh = ToElementMesh[Cylinder[{{0, 0, 0}, {0, 0, L}}, r], 
   MaxCellMeasure -> .5];
mesh["Wireframe"]

Mathematica graphics

Here I calculate the eigenvalues and eigenvectors. The first 6 are zero (approximatly) because they are rigid body modes.

{vals, vecs} = 
  NDEigensystem[
   stressOperator[Y, \[Nu]], {u, v, w}, {x, y, z} \[Element] mesh, 14];
freqs = 1/(2 \[Pi]) Sqrt[Abs[Chop[vals]]];
TableForm[freqs]

Mathematica graphics

Now I have the problem of animating the eigenvectors. The following module does this but is wrong because I am animating every node. Only the boundary nodes can be seen and thus only these should be animated. If I make too many animations using this DynamicModule then Mathematica freezes.

    ClearAll[animateEigenVectors];
animateEigenVectors[{freqs_, vecs_}, n_] := DynamicModule[{dmesh},
  Animate[
   dmesh = 
    ElementMeshDeformation[mesh, vecs[[n]], 
     "ScalingFactor" -> 0.01 Cos[t]];
   Show[{dmesh[
      "Wireframe"[
       "ElementMeshDirective" -> 
        Directive[EdgeForm[Gray], FaceForm[Lighter[Blue]]]]]},
    PlotLabel -> 
     "Natural Frequency = " <> ToString[freqs[[n]]] <> 
      " Hz  Mode = " <> ToString[n], Boxed -> False,
    PlotRange -> 1.5 {{-r, r}, {-r, r}, {0, L}}
    ],
   {t, 0, 2 \[Pi]}, SaveDefinitions -> True]
  ]


animateEigenVectors[{freqs, vecs}, 7]

Mathematica graphics

(I can't do the fancy animation bit in SE yet.) Is it possible to extract and animate just the boundary elements? Thanks

$\endgroup$
  • $\begingroup$ Could you provide the code you used to make the mesh? $\endgroup$ – user21 Dec 9 '15 at 17:00
  • $\begingroup$ @user21 Rats I forgot two lines of code. Now edited in. Thanks for looking at this. $\endgroup$ – Hugh Dec 9 '15 at 21:31
1
$\begingroup$

Try this:

ClearAll[animateEigenVectors];
animateEigenVectors[{freqs_, vecs_}, n_] := 
 DynamicModule[{dmesh}, Animate[
   (*dmesh=ElementMeshDeformation[mesh,vecs[[n]],
   "ScalingFactor"\[Rule]0.01 Cos[t]];*)

   scalefactor = 0.01 Cos[t];
   move = 
    scalefactor*
     Transpose[
      Developer`ToPackedArray[#["ValuesOnGrid"] & /@ vecs[[n]]]];
   dmesh = 
    ToBoundaryMesh["Coordinates" -> mesh["Coordinates"] + move, 
     "BoundaryElements" -> mesh["BoundaryElements"], 
     "CheckIntersections" -> False];

   Show[{dmesh[
      "Wireframe"[
       "ElementMeshDirective" -> 
        Directive[EdgeForm[Gray], FaceForm[Lighter[Blue]]]]]}, 
    PlotLabel -> 
     "Natural Frequency = " <> ToString[freqs[[n]]] <> 
      " Hz  Mode = " <> ToString[n], Boxed -> False, 
    PlotRange -> 1.5 {{-r, r}, {-r, r}, {0, L}}], {t, 0, 2 \[Pi]}, 
   SaveDefinitions -> True]]

This generates a boundary mesh only and switches off the expensive intersection checking. Since I could not reproduce the original problem, let me know if this helps to solve the issue you see.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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