### Eigenvalues of cantilever beam in linear elasticity problem

317
views
0
3 months ago by
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:

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):

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
1
Are you sure it's the static problem's eigenvalues you're interested in? What about the mass matrix?
written 3 months ago by klunkean
Thanks for pointing this out. I thought for a continuous beam the mass matrix could be diagonal and would not affect the eigenvalues. Now I tried to include a mass matrix but I am stuck with how to do this in FEniCS. I thought I can define the mass matrix as

L = rho * u * v * dx​

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

Invalid ranks 1 and 1 in product.​
Could you point me into the right direction of how to assamble the mass matrix correctly?
written 3 months ago by Daniel
1
The notion of * 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:
a_M = rho*inner(u,v)*dx​
Looking at this wikipedia article, I think the Eigenvalue-Problem looks like
$Ku=\omega^2Mu$Ku=ω2Mu
since 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 .
written 3 months ago by klunkean
Thanks for your help. Just modified the code such that:
a = inner(sigma(u), epsilon(v))*dx
b = rho * inner(u, v) * dx

A = PETScMatrix()
B = PETScMatrix()
assemble(a, tensor=A)
assemble(b, tensor=B)
bc.zero(A)
bc.zero(B)

eigensolver = SLEPcEigenSolver(A, B)
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​
However, it looks like FEniCS does not compute the eigenvalues now. Since it returns the message
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-37-5a111b7fb220> in <module>()
2 w = np.zeros(n)
3 for i in range(n):
----> 4     r, c, rx, cx = eigensolver.get_eigenpair(i)
5     w[i] = r
6     u_tmp = Function(V)

/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
***
***
*** 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
*** -------------------------------------------------------------------------​
Am I on the wrong track how to use SLEPcEigenSolver?
written 3 months ago by Daniel
1
print eigensolver.get_number_converged()​

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

written 3 months ago by klunkean
Thanks for your help. I understand much better now what is going on. I will work on a solution and post it here.
written 3 months ago by Daniel

4
3 months ago by
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):

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

written 11 days ago by carlos lugo
2
3 months ago by
Dear all,
I have just written a week ago a detailed commented demo on this very problem: