raw:: html

This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License

# .. _ModalAnalysis:
#
# ==========================================
# Modal analysis of an elastic structure
# ==========================================
#
# -------------
# Introduction
# -------------
#
# This program performs a dynamic modal analysis of an elastic cantilever beam
# represented by a 3D solid continuum. The eigenmodes are computed using the
# **SLEPcEigensolver** and compared against an analytical solution of beam theory.
# The corresponding file can be obtained from :download:`cantilever_modal.py`.
#
#
# The first four eigenmodes of this demo will look as follows:
#
# .. image:: vibration_modes.gif
# :scale: 80 %
#
# The first two fundamental modes are on top with bending along the weak axis (left) and along
# the strong axis (right), the next two modes are at the bottom.
#
# ---------------
# Implementation
# ---------------
#
# After importing the relevant modules, the geometry of a beam of length :math:`L=20`
# and rectangular section of size :math:`B\times H` with :math:`B=0.5, H=1` is first defined::
from fenics import *
import numpy as np
L, B, H = 20., 0.5, 1.
Nx = 200
Ny = int(B/L*Nx)+1
Nz = int(H/L*Nx)+1
mesh = BoxMesh(Point(0.,0.,0.),Point(L,B,H), Nx, Ny, Nz)
# Material parameters and elastic constitutive relations are classical (here we
# take :math:`\nu=0`) and we also introduce the material density :math:`\rho` for
# later definition of the mass matrix::
E, nu = 1e5, 0.
rho = 1e-3
# Lame coefficient for constitutive relation
mu = E/2./(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
def eps(v):
return sym(grad(v))
def sigma(v):
dim = v.geometric_dimension()
return 2.0*mu*eps(v) + lmbda*tr(eps(v))*Identity(dim)
# Standard FunctionSpace is defined and boundary conditions correspond to a
# fully clamped support at :math:`x=0`::
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)
# The system stiffness matrix :math:`[K]` and mass matrix :math:`[M]` are
# respectively obtained from assembling the corresponding variational forms::
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)
# Matrices :math:`[K]` and :math:`[M]` are first defined as PETSc Matrix and
# forms are assembled into it to ensure that they have the right type.
# Note that boundary conditions have been applied to the stiffness matrix using
# ``assemble_system`` so as to preserve symmetry (a dummy ``l_form`` and right-hand side
# vector have been introduced to call this function).
#
#
# Modal dynamic analysis consists in solving the following generalized
# eigenvalue problem :math:`[K]\{U\}=\lambda[M]\{U\}` where the eigenvalue
# is related to the eigenfrequency :math:`\lambda=\omega^2`. This problem
# can be solved using the ``SLEPcEigenSolver``. ::
eigensolver = SLEPcEigenSolver(K, M)
eigensolver.parameters['problem_type'] = 'gen_hermitian'
eigensolver.parameters["spectrum"] = "smallest real"
eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
eigensolver.parameters['spectral_shift'] = 0.
# The problem type is specified to be a generalized eigenvalue problem with
# Hermitian matrices. By default, SLEPc computes the largest eigenvalues, here
# we instead look for the smallest eigenvalues (they should all be real). To
# improve convergence of the eigensolver for finding the smallest eigenvalues
# (by default it computes the largest ones), a spectral transform is performed
# using the keyword ``shift-invert`` i.e. the original problem is transformed into
# an equivalent problem with eigenvalues given by :math:`\dfrac{1}{\lambda - \sigma}`
# instead of :math:`\lambda` where :math:`\sigma` is the value of the spectral shift.
# It is therefore much easier to compute eigenvalues close to :math:`\sigma` i.e.
# close to :math:`\sigma = 0` in the present case. Eigenvalues are then
# transformed back by SLEPc to their original value :math:`\lambda`.
#
#
# We now ask SLEPc to extract the first 6 eigenvalues by calling its solve function
# and extract the corresponding eigenpair (first two arguments of ``get_eigenpair``
# correspond to the real and complex part of the eigenvalue, the last two to the
# real and complex part of the eigenvector)::
N_eig = 6 # number of eigenvalues
print "Computing %i first eigenvalues..." % N_eig
eigensolver.solve(N_eig)
# Exact solution computation
from scipy.optimize import root
from math import cos, cosh
falpha = lambda x: cos(x)*cosh(x)+1
alpha = lambda n: root(falpha, (2*n+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)
# 3D eigenfrequency
freq_3D = sqrt(r)/2/pi
# Beam eigenfrequency
if i % 2 == 0: # exact solution should correspond to weak axis bending
I_bend = H*B**3/12.
else: #exact solution should correspond to strong axis bending
I_bend = B*H**3/12.
freq_beam = alpha(i/2)**2*sqrt(E*I_bend/(rho*B*H*L**4))/2/pi
print("Solid FE: {0:8.5f} [Hz] Beam theory: {1:8.5f} [Hz]".format(freq_3D, freq_beam))
# Initialize function and assign eigenvector (renormalize by stiffness matrix)
eigenmode = Function(V,name="Eigenvector "+str(i))
eigenmode.vector()[:] = rx
# The beam analytical solution is obtained using the eigenfrequencies of a clamped
# beam in bending given by :math:`\omega_n = \alpha_n^2\sqrt{\dfrac{EI}{\rho S L^4}}`
# where :math:`S=BH` is the beam section, :math:`I` the bending inertia and
# :math:`\alpha_n` is the solution of the following nonlinear equation:
#
# .. math::
# \cos(\alpha)\cosh(\alpha)+1 = 0
#
# the solution of which can be well approximated by :math:`(2n+1)\pi/2` for :math:`n\geq 3`.
# Since the beam possesses two bending axis, each solution to the previous equation is
# associated with two frequencies, one with bending along the weak axis (:math:`I=I_{\text{weak}} = HB^3/12`)
# and the other along the strong axis (:math:`I=I_{\text{strong}} = BH^3/12`). Since :math:`I_{\text{strong}} = 4I_{\text{weak}}`
# for the considered numerical values, the strong axis bending frequency will be twice that corresponsing
# to bending along the weak axis. The solution :math:`\alpha_n` are computed using the
# ``scipy.optimize.root`` function with initial guess given by :math:`(2n+1)\pi/2`.
#
# With ``Nx=400``, we obtain the following comparison between the FE eigenfrequencies
# and the beam theory eigenfrequencies :
#
#
# ===== ============= =================
# Mode Eigenfrequencies
# ----- --------------------------------
# # Solid FE [Hz] Beam theory [Hz]
# ===== ============= =================
# 1 2.04991 2.01925
# 2 4.04854 4.03850
# 3 12.81504 12.65443
# 4 25.12717 25.30886
# 5 35.74168 35.43277
# 6 66.94816 70.86554
# ===== ============= =================
#
#