cantilever_modal.py.rst 7.07 KB
Newer Older
Jeremy BLEYER's avatar
Jeremy BLEYER committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 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 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194

.. _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/omega

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 to two frequencies with bending along the weak axis (:math:`I=I_{\text{weak}} = HB^3/12`)
and 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 
=====  =============  =================