Beam modal shape extraction based on 'Modal analysis of an elastic structure' example

[Matt141]

Hi everyone,

I am trying to extract beam modal shapes using the example ‘Modal analysis of an elastic structure’ example, cantilever_modal.py, at this address :
https://comet-fenics.readthedocs.io/en/latest/demo/modal_analysis_dynamics/cantilever_modal.py.html

I am using Ubuntu 20, python 3.7 and fenics 2018.1.0 (which seems to be the verison used according to intro page.

I don’t encounter code issue, but the mode shape displayed is wrong.
I am expecting a beam bending mode.

Any suggestion why I am not getting the expected results is welcomed.

Regards,
Matt

Here is the code:

from dolfin import *
import numpy as np
import matplotlib.pyplot as plt

L, B, H = 20., 0.5, 1.

Nx = 200
Ny = int(B/LNx)+1
Nz = int(H/LNx)+1

mesh = BoxMesh(Point(0.,0.,0.),Point(L,B,H), Nx, Ny, Nz)

E, nu = Constant(1e5), Constant(0.)
rho = Constant(1e-3)

Lame coefficient for constitutive relation

mu = E/2./(1+nu)
lmbda = Enu/(1+nu)/(1-2nu)

def eps(v):
return sym(grad(v))
def sigma(v):
dim = v.geometric_dimension()
return 2.0mueps(v) + lmbda*tr(eps(v))*Identity(dim)

V = VectorFunctionSpace(mesh, ‘Lagrange’, degree=1)
u_ = TrialFunction(V)
du = TestFunction(V)

def left(x, on_boundary):
return near(x[0],0.)

bc = DirichletBC(V, Constant((0.,0.,0.)), left)

k_form = inner(sigma(du),eps(u_))*dx
l_form = Constant(1.)*u_[0]*dx
K = PETScMatrix()
b = PETScVector()
assemble_system(k_form, l_form, bc, A_tensor=K, b_tensor=b)

m_form = rho*dot(du,u_)*dx
M = PETScMatrix()
assemble(m_form, tensor=M)

eigensolver = SLEPcEigenSolver(K, M)
eigensolver.parameters[‘problem_type’] = ‘gen_hermitian’
eigensolver.parameters[‘spectral_transform’] = ‘shift-and-invert’
eigensolver.parameters[‘spectral_shift’] = 0.

N_eig = 6 # number of eigenvalues
print(“Computing {} first eigenvalues…”.format(N_eig))
eigensolver.solve(N_eig)

from scipy.optimize import root
from math import cos, cosh
falpha = lambda x: cos(x)cosh(x)+1
alpha = lambda n: root(falpha, (2n+1)*pi/2.)[‘x’][0]

Set up file for exporting results

file_results = XDMFFile(“modal_analysis.xdmf”)
file_results.parameters[“flush_output”] = True
file_results.parameters[“functions_share_mesh”] = True

Extraction

for i in range(N_eig):
# Extract eigenpair
r, c, rx, cx = eigensolver.get_eigenpair(i)

freq_3D = sqrt(r)/2/pi

Initialize function and assign eigenvector

eigenmode = Function(V,name="Eigenvector "+str(i))
eigenmode.vector()[:] = rx

plot(eigenmode)
plt.show()

[Matt141]

Here is the wrong result of the first mode:

[bleyerj]

Hi,
there is no proper 3D plotting support inside dolfin. Use Paraview for instance to plot the eigenmode.

[Matt141]

Thanks you for the answer. I finally managed to get Paraview working, and indeed the 3D plot is better.