### Eigenvalues of cantilever beam in linear elasticity problem

317

views

0

To get started with FEniCS I am trying to compute eigenvalues/eigenfrequencies of a cantilever beam by modifying the linear elasticity example from the FEniCS tutorial. I came up with the following code:

The computed eigenvalues normalized to the first eigenvalue are

This is far away from the analytic values. For example the second value in this table should be around 6.2. Can someone point me into the right direction how to compute the eigenvalues of a linear elasticity problem?

```
from fenics import *
import numpy as np
L = 120e-6
b = 30e-6
h = 3.5e-6
E = 170e9
nu = 0.28
rho = 2329
mu = E / (2 * (1 + nu))
lambda_ = E * nu / ((1 + nu) * (1 - 2 * nu))
mesh = BoxMesh(Point(0, 0, 0), Point(L, b, h), 20, 10, 3)
File('export/mesh.pvd') << mesh
V = VectorFunctionSpace(mesh, 'P', 1)
tol = 1E-14
def clamped_boundary(x, on_boundary):
return on_boundary and x[0] < tol
bc = DirichletBC(V, Constant((0, 0, 0)), clamped_boundary)
def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
#return sym(nabla_grad(u))
def sigma(u):
return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)
u = TrialFunction(V)
v = TestFunction(V)
d = u.geometric_dimension()
a = inner(sigma(u), epsilon(v))*dx
A = PETScMatrix()
assemble(a, tensor=A)
bc.zero(A)
eigensolver = SLEPcEigenSolver(A)
eigensolver.solve()
n = 10
w = np.zeros(n)
for i in range(n):
r, c, rx, cx = eigensolver.get_eigenpair(i)
w[i] = r
u_tmp = Function(V)
u_tmp.vector()[:] = rx
File('export/mode_%01i.pvd' % i) << u_tmp
for el in w:
print(el / w[0])
```

The computed eigenvalues normalized to the first eigenvalue are

```
1.0
0.998557497337
0.996217990034
0.99307813208
0.989280003902
0.986539123292
0.985217624586
0.985037096289
0.982968521832
0.98101003069
```

This is far away from the analytic values. For example the second value in this table should be around 6.2. Can someone point me into the right direction how to compute the eigenvalues of a linear elasticity problem?

Community: FEniCS Project

### 2 Answers

4

Building on klunkean's comments above, I tinkered with the parameters a bit and the second frequency (i.e., the square root of the second eigenvalue) is now around 6.2, as expected:

```
from fenics import *
import numpy as np
L = 120e-6
b = 30e-6
h = 3.5e-6
E = 170e9
nu = 0.28
rho = 2329
mu = E / (2 * (1 + nu))
lambda_ = E * nu / ((1 + nu) * (1 - 2 * nu))
# Used fewer, higher-order elements.
mesh = BoxMesh(Point(0, 0, 0), Point(L, b, h), 10, 5, 2)
File('export/mesh.pvd') << mesh
V = VectorFunctionSpace(mesh, 'P', 3)
tol = 1E-14
def clamped_boundary(x, on_boundary):
return on_boundary and x[0] < tol
bc = DirichletBC(V, Constant((0, 0, 0)), clamped_boundary)
def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
def sigma(u):
return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)
u = TrialFunction(V)
v = TestFunction(V)
d = u.geometric_dimension()
a = inner(sigma(u), epsilon(v))*dx
A = PETScMatrix()
assemble(a, tensor=A)
# Shift eigenvalues corresponding to BCs to the high end of the
# spectrum by applying BCs with a very large number on the diagonal.
bc.zero(A)
bc.zero_columns(A,Function(V).vector(),1e20)
B = PETScMatrix()
assemble(rho*inner(u,v)*dx, tensor=B)
bc.apply(B)
eigensolver = SLEPcEigenSolver(A,B)
# Most reliable for small systems and useful for debugging, but does not
# scale up well (LAPACK is for dense matrices)
#eigensolver.parameters["solver"] = "lapack"
# Various choices of iterative solvers listed here:
#
# https://fenicsproject.org/olddocs/dolfin/1.3.0/python/programmers-reference/cpp/la/SLEPcEigenSolver.html
#
eigensolver.parameters["solver"] = "arnoldi"
# You usually want the lowest frequencies in vibrational problems:
eigensolver.parameters["spectrum"]="smallest magnitude"
# Since we are interested in the low frequencies, it improves convergence
# of iterative methods to look for eigenvalues of $A^{-1}$ instead of $A$.
eigensolver.parameters["spectral_transform"] = "shift-and-invert"
eigensolver.parameters["spectral_shift"] = 0.0
# Can specify the number of modes to look for:
n = 10
eigensolver.solve(n)
w = np.zeros(n)
import math
for i in range(n):
r, c, rx, cx = eigensolver.get_eigenpair(i)
# Eigenproblem gives you frequency squared:
w[i] = math.sqrt(r)
u_tmp = Function(V)
u_tmp.vector()[:] = rx
File('export/mode_%01i.pvd' % i) << u_tmp
for el in w:
print(el / w[0])
```

Thanks for you for your reply. I just took your code and tried to run it but I got the error message

```
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
<ipython-input-1-3e4f67e5e70d> in <module>()
73 import math
74 for i in range(n):
---> 75 r, c, rx, cx = eigensolver.get_eigenpair(i)
76
77 # Eigenproblem gives you frequency squared:
/usr/lib/python3/dist-packages/dolfin/cpp/la.py in get_eigenpair(self, i, r_vec, c_vec)
4410 r_vec = r_vec or PETScVector()
4411 c_vec = c_vec or PETScVector()
-> 4412 lr, lc = self._get_eigenpair(r_vec, c_vec, i)
4413 return lr, lc, r_vec, c_vec
4414
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to extract eigenpair from SLEPc eigenvalue solver.
*** Reason: Requested eigenpair (0) has not been computed.
*** Where: This error was encountered inside SLEPcEigenSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2017.2.0
*** Git changeset: 4c59bbdb45b95db2f07f4e3fd8985c098615527f
*** -------------------------------------------------------------------------
```

I am running FEniCS in Jupyter in a docker container on Windows 10. Could that influence the behaviour of the solver?
written
3 months ago by
Daniel

1

I initially ran it on Ubuntu, with version 2017.2 (installed via apt-get). I tried the stable:latest docker container, and get the same error as you. I went back and tried the stable:2016.2.0 container

quay.io/fenicsproject/stable:2016.2.0

and ran the same script (with python, not python3) and it worked... It may depend on what version of SLEPc is installed. My computer and the 2016.2.0 container have SLEPc 3.7, whereas the 2017.2 container has SLEPc 3.8.

quay.io/fenicsproject/stable:2016.2.0

and ran the same script (with python, not python3) and it worked... It may depend on what version of SLEPc is installed. My computer and the 2016.2.0 container have SLEPc 3.7, whereas the 2017.2 container has SLEPc 3.8.

written
3 months ago by
David Kamensky

I can confirm this behaviour. For stable:2017.1.0 the script runs through and gives the correct results, for stable:2017.2.0 it does not.

written
3 months ago by
Daniel

It does depend on the version of SLEPc if I recall correctly. A change was made a while ago where the spectral shift could not produce a singularity. Try a shift close to where you expect the first eigenvalue instead of 0.0.

Also make sure the shift is working correctly (run with

Also, if you're shifting the spectrum, experiment with not sorting the eigenvalues by magnitude.

Also make sure the shift is working correctly (run with

`eps_view`

).Also, if you're shifting the spectrum, experiment with not sorting the eigenvalues by magnitude.

written
3 months ago by
Nate

Has anybody found a solution for this issue? I can confirm that the the eigenvalue problem for a mode analysis of an euler-buckling beam does not converge (no eigenpairs) for stable 2017.2.0.

Simpler eigensystems do converge.

Simpler eigensystems do converge.

written
11 days ago by
carlos lugo

2

Dear all,

I have just written a week ago a detailed commented demo on this very problem:

http://comet-fenics.readthedocs.io/en/latest/demo/modal_analysis_dynamics/cantilever_modal.py.html

Best regards
I have followed your solution as it is in the post yet the issue of not convergence remains. However it appears it has to do with the version used to solve the eigenvalue problem.

Here https://www.allanswered.com/post/erpjl/#erpxl passing arguments to the solver explicitly using the slepc4py API seems to work!

I have just written a week ago a detailed commented demo on this very problem:

http://comet-fenics.readthedocs.io/en/latest/demo/modal_analysis_dynamics/cantilever_modal.py.html

Best regards

Here https://www.allanswered.com/post/erpjl/#erpxl passing arguments to the solver explicitly using the slepc4py API seems to work!

written
11 days ago by
carlos lugo

Please login to add an answer/comment or follow this question.

However, when I run this line of code I get the error message

Could you point me into the right direction of how to assamble the mass matrix correctly?

`*`

is ambiguuous for vectorial quantities and thus not properly defined. You'll want to use`dot`

or`inner`

which for vectors both refer to the standard scalar product:Looking at this wikipedia article, I think the Eigenvalue-Problem looks like

$Ku=\omega^2Mu$

Ku=ω^{2}Musince for the displacement field an exponential ansatz proportional to $\text{exp}\left(i\omega t\right)$exp(

iωt) is employed. The eigenvalues should be $\lambda=-\omega^2$λ=−ω^{2}.However, it looks like FEniCS does not compute the eigenvalues now. Since it returns the message

Am I on the wrong track how to use SLEPcEigenSolver?

after invoking the solve method returns zero. So the solver doesn't converge, which means there are no eigenpairs to extract. Why it doesn't converge, well I don't know...