18
$\begingroup$

There are many examples of eigenmodes computations for surfaces with Mathematica, such as:

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}}];

enter image description here

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...

Improve this question
$\endgroup$

1 Answer 1

Reset to default
17
$\begingroup$

I did not get it to work with your Navier-Cauchy operator, not sure if it is correct. Here is an approach with a standard stress operator:

Needs["NDSolve`FEM`"] 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}]} material = {nu -> 1/3, Y -> 200*10^9}; rls = {lambda -> nu*Y/((1 + nu)*(1 - 2 nu)), mu -> Y/(2 (1 + nu))} /. material; uVec = Through[{u, v, w}[x, y, z]]; vars = {x, y, z}; pde = stressOperator[Y, nu] /. material; (*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}]/.rls;*)

Use a somewhat finer mesh than the default. I use a first order mesh because a second order mesh would take more time to compute.

ir = ImplicitRegion[0.8 <= x^2 + y^2 <= 1 && 0 <= z <= 3, {x, y, z}]; mesh = ToElementMesh[ir, {{-1, 1}, {-1, 1}, {0, 2}}, "MeshOrder" -> 1, MaxCellMeasure -> 0.0005] nrEigen = 12; {vals, vecs} = NDEigensystem[ Join[pde, DirichletCondition[# == 0, z == 0] & /@ uVec], {u, v, w}, {x, y, z} \[Element] mesh, nrEigen]; GraphicsGrid[Partition[Show[ MeshRegion[ ElementMeshDeformation[mesh, vecs[[#]], "ScalingFactor" -> 1/5]], ElementMeshEdgeframe3D[mesh]] & /@ Range[nrEigen], 3]]

enter image description here

If you get it to work with a Navier-Cauchy operator, let me know.

There are a few other questions related to this, for example here and here.

enter image description here

Update:

Let's reconsider the Cauchy-Navier operator. Here is a different form of it:

pde2 = -((lambda + mu) Grad[ Div[{u[x, y, z], v[x, y, z], w[x, y, z]}, {x, y, z}], {x, y, z}] + mu Laplacian[{u[x, y, z], v[x, y, z], w[x, y, z]}, {x, y, z}]) /. rls

In a first step one might compare this to the stress operator above:

Thread[Activate[pde] == pde2] // Simplify {True, True, True}

This, however, is not the whole story. Let's look at what actually gets parsed for the inactive stress operator:

{state} = NDSolve`ProcessEquations[ Join[Thread[pde1 == {0, 0, 0}], DirichletCondition[# == 0, z == 0] & /@ uVec], {u, v, w}, {x, y, z} \[Element] mesh]; pdec = state["FiniteElementData"]["PDECoefficientData"]; pdec["DiffusionCoefficients"] // MatrixForm

enter image description here

Now, compare that to the activated or Navier-Cauchy operator:

{state} = NDSolve`ProcessEquations[ Join[Thread[pde2 == {0, 0, 0}], DirichletCondition[# == 0, z == 0] & /@ uVec], {u, v, w}, {x, y, z} \[Element] mesh]; pdec = state["FiniteElementData"]["PDECoefficientData"]; pdec["DiffusionCoefficients"] // MatrixForm

enter image description here

Now, closely look at the off diagonal elements and compare them. The inactive stress operator is able to retain the the skew off-diagonal elements of the operator while this form of the Cauchy-Navier operator is not. I currently do not see a way to make the Cauchy-Navier operator behave the same way as the inactive stress operator. So my advice is to use the inactive stress operator for 3D stress problems.

Update:

This behavior is expected and, in fact, was one of the main motivations for introducing Inactive. You can find this documented and explained with a 2D example in the section Formal Partial Differential Equations.

Improve this answer
$\endgroup$
14
  • 1
    $\begingroup$ Hope you don't mind the edit, couldn't resist :) $\endgroup$ Commented Jan 31, 2018 at 9:15
  • $\begingroup$ @Kuba Yea! Great! Love it - reminds me of that dancing baby. We need sound on this page! Thx. $\endgroup$ Commented Jan 31, 2018 at 9:36
  • $\begingroup$ This is what I get using my pde and GraphicsGrid[ Partition[ Show[MeshRegion[ ElementMeshDeformation[mesh, vecs[[#]], "ScalingFactor" -> 1/5]]] & /@ Range[nrEigen], 3]]: imgur.com/a/YIeDO. Looks different but looks "a bit" like modes too... I'll double check later. (and thanks!) $\endgroup$ Commented Jan 31, 2018 at 10:41
  • $\begingroup$ @anderstood, what platform and M version are you using, what happens with a finer mesh? $\endgroup$ Commented Jan 31, 2018 at 10:48
  • $\begingroup$ I'm using MMA 11.0 with Ubuntu 16.04. The result does not converge easily, even with MaxCellMeasure -> 0.0001. I'm not sure what is wrong with my pde. I'll try on a plate to verify. $\endgroup$ Commented Jan 31, 2018 at 15:46

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.