diff --git a/doc/_static/logo_Navier.png b/doc/_static/logo_Navier.png new file mode 100644 index 0000000000000000000000000000000000000000..960089fb7881cdc97fda8194f96525be20177cbf Binary files /dev/null and b/doc/_static/logo_Navier.png differ diff --git a/doc/_static/logo_tutelles.png b/doc/_static/logo_tutelles.png new file mode 100644 index 0000000000000000000000000000000000000000..407535e38a9c6f7acee9e037eb13d08aa3a5d0d2 Binary files /dev/null and b/doc/_static/logo_tutelles.png differ diff --git a/doc/index.rst b/doc/index.rst index a6826774e1bfa90dbd6bd27d75e34ccd1ebbcacb..3018a52e6754eacfcefe3eba03b3d22356b23874 100644 --- a/doc/index.rst +++ b/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 diff --git a/doc/intro.rst b/doc/intro.rst index 199d0bc12e8764b93e03c74fb13754b1049f1233..a04b8287c5e538d353e647dc29aa185812022648 100644 --- a/doc/intro.rst +++ b/doc/intro.rst @@ -48,6 +48,13 @@ in Solid and Structural Mechanics at `Laboratoire Navier a joint research unit of `Ecole Nationale des Ponts et Chaussées `_, `IFSTTAR `_ and `CNRS `_ (UMR 8205). +.. image:: _static/logo_Navier.png + :scale: 8 % + :align: left +.. image:: _static/logo_tutelles.png + :scale: 20 % + :align: right + diff --git a/doc/rstprocess.py b/doc/rstprocess.py index 7364782483967f5e124a206cd34880d09597db64..effbca79ea620c9c995f338bed8474867e840b57 100644 --- a/doc/rstprocess.py +++ b/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)) diff --git a/examples/modal_analysis_dynamics/cantilever_modal.py.rst b/examples/modal_analysis_dynamics/cantilever_modal.py.rst new file mode 100644 index 0000000000000000000000000000000000000000..4aed8722c2476332148d423a54715626281adfc8 --- /dev/null +++ b/examples/modal_analysis_dynamics/cantilever_modal.py.rst @@ -0,0 +1,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 +===== ============= ================= + + diff --git a/examples/modal_analysis_dynamics/vibration_modes.gif b/examples/modal_analysis_dynamics/vibration_modes.gif new file mode 100644 index 0000000000000000000000000000000000000000..1ae977159786e93bf97338ee01ac728c9f7a961e Binary files /dev/null and b/examples/modal_analysis_dynamics/vibration_modes.gif differ