cantilever_modal.py 7.55 KB
Newer Older
1 2 3 4
raw:: html

<a rel="license" href="http://creativecommons.org/licenses/by-sa/4.0/"><p align="center"><img alt="Creative Commons License" style="border-width:0" src="https://i.creativecommons.org/l/by-sa/4.0/88x31.png"/></a><br />This work is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by-sa/4.0/">Creative Commons Attribution-ShareAlike 4.0 International License</a></p>

Jeremy BLEYER's avatar
Jeremy BLEYER committed
5 6 7 8 9 10 11 12 13 14 15
# .. _ModalAnalysis:
#
# ==========================================
# Modal analysis of an elastic structure
# ==========================================
#
# -------------
# Introduction
# -------------
#
# This program performs a dynamic modal analysis of an elastic cantilever beam
Jeremy BLEYER's avatar
Jeremy BLEYER committed
16 17
# represented by a 3D solid continuum. The eigenmodes are computed using the
# **SLEPcEigensolver** and compared against an analytical solution of beam theory.
Jeremy BLEYER's avatar
Jeremy BLEYER committed
18 19 20 21 22 23 24 25
# 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 %
#
Jeremy BLEYER's avatar
Jeremy BLEYER committed
26 27
# 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.
Jeremy BLEYER's avatar
Jeremy BLEYER committed
28 29 30 31 32
#
# ---------------
# Implementation
# ---------------
#
Jeremy BLEYER's avatar
Jeremy BLEYER committed
33
# After importing the relevant modules, the geometry of a beam of length :math:`L=20`
Jeremy BLEYER's avatar
Jeremy BLEYER committed
34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
# 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
Jeremy BLEYER's avatar
Jeremy BLEYER committed
49
# take :math:`\nu=0`) and we also introduce the material density :math:`\rho` for
Jeremy BLEYER's avatar
Jeremy BLEYER committed
50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
# 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)

Jeremy BLEYER's avatar
Jeremy BLEYER committed
65
# Standard FunctionSpace is defined and boundary conditions correspond to a
Jeremy BLEYER's avatar
Jeremy BLEYER committed
66 67 68 69 70 71 72 73 74 75 76 77 78
# 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)


Jeremy BLEYER's avatar
Jeremy BLEYER committed
79
# The system stiffness matrix :math:`[K]` and mass matrix :math:`[M]` are
Jeremy BLEYER's avatar
Jeremy BLEYER committed
80 81 82 83 84 85 86 87 88 89 90 91
# 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)

Jeremy BLEYER's avatar
Jeremy BLEYER committed
92
# Matrices :math:`[K]` and :math:`[M]` are first defined as PETSc Matrix and
Jeremy BLEYER's avatar
Jeremy BLEYER committed
93 94 95
# 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
Jeremy BLEYER's avatar
Jeremy BLEYER committed
96
# vector have been introduced to call this function).
Jeremy BLEYER's avatar
Jeremy BLEYER committed
97 98
#
#
Jeremy BLEYER's avatar
Jeremy BLEYER committed
99
# Modal dynamic analysis consists in solving the following generalized
Jeremy BLEYER's avatar
Jeremy BLEYER committed
100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126
# 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)::
Jeremy BLEYER's avatar
Jeremy BLEYER committed
127

Jeremy BLEYER's avatar
Jeremy BLEYER committed
128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146
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)
Jeremy BLEYER's avatar
Jeremy BLEYER committed
147

Jeremy BLEYER's avatar
Jeremy BLEYER committed
148 149
    # 3D eigenfrequency
    freq_3D = sqrt(r)/2/pi
Jeremy BLEYER's avatar
Jeremy BLEYER committed
150

Jeremy BLEYER's avatar
Jeremy BLEYER committed
151 152 153 154 155 156
    # 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
Jeremy BLEYER's avatar
Jeremy BLEYER committed
157

Jeremy BLEYER's avatar
Jeremy BLEYER committed
158 159 160 161
    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))
Jeremy BLEYER's avatar
Jeremy BLEYER committed
162
    eigenmode.vector()[:] = rx
Jeremy BLEYER's avatar
Jeremy BLEYER committed
163 164 165

# 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}}`
Jeremy BLEYER's avatar
Jeremy BLEYER committed
166
# where :math:`S=BH` is the beam section, :math:`I` the bending inertia and
Jeremy BLEYER's avatar
Jeremy BLEYER committed
167 168 169 170 171 172 173 174 175 176 177 178
# :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`.
Jeremy BLEYER's avatar
Jeremy BLEYER committed
179
#
Jeremy BLEYER's avatar
Jeremy BLEYER committed
180
# With ``Nx=400``, we obtain the following comparison between the FE eigenfrequencies
Jeremy BLEYER's avatar
Jeremy BLEYER committed
181
# and the beam theory eigenfrequencies :
Jeremy BLEYER's avatar
Jeremy BLEYER committed
182 183 184 185 186 187 188 189
#
#
# =====  =============  =================
# Mode      Eigenfrequencies
# -----  --------------------------------
#  #     Solid FE [Hz]   Beam theory [Hz]
# =====  =============  =================
#   1      2.04991           2.01925
Jeremy BLEYER's avatar
Jeremy BLEYER committed
190
#   2      4.04854           4.03850
Jeremy BLEYER's avatar
Jeremy BLEYER committed
191
#   3      12.81504         12.65443
Jeremy BLEYER's avatar
Jeremy BLEYER committed
192
#   4      25.12717         25.30886
Jeremy BLEYER's avatar
Jeremy BLEYER committed
193
#   5      35.74168         35.43277
Jeremy BLEYER's avatar
Jeremy BLEYER committed
194
#   6      66.94816         70.86554
Jeremy BLEYER's avatar
Jeremy BLEYER committed
195 196 197
# =====  =============  =================
#
#