Dynamic modal analysis example

doc/index.rst

@@ -13,6 +13,7 @@ Contents:
 
    intro
    2D_elasticity
+   demo/modal_analysis_dynamics/cantilever_modal.py.rst
    demo/reissner_mindlin/reissner_mindlin.rst
 
 

doc/intro.rst

@@ -48,6 +48,13 @@ in Solid and Structural Mechanics at `Laboratoire Navier <http://navier.enpc.fr>
 a joint research unit of `Ecole Nationale des Ponts et Chaussées <http://www.enpc.fr>`_, 
 `IFSTTAR <http://www.ifsttar.fr>`_ and `CNRS <http://www.cnrs.fr>`_ (UMR 8205).
 
+.. image:: _static/logo_Navier.png
+   :scale: 8 %
+   :align: left
+.. image:: _static/logo_tutelles.png
+   :scale: 20 %
+   :align: right
+
 
 
 

doc/rstprocess.py

@@ -72,7 +72,7 @@ def process():
                         if not ret == 0:
                             raise RuntimeError("Unable to convert rst file to a .py ({})".format(f))
 
-                png_files = [f for f in os.listdir(version_path) if os.path.splitext(f)[1] == ".png" ]
+                png_files = [f for f in os.listdir(version_path) if os.path.splitext(f)[1] in [".png",".gif"] ]
                 for f in png_files:
                     source = os.path.join(version_path, f)
                     print("Copying {} to {}".format(source, demo_dir)) 

examples/modal_analysis_dynamics/cantilever_modal.py.rst

+
+.. _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 
+=====  =============  =================
+
+