cantilever_modal.py 7.16 KB
 Jeremy BLEYER committed May 28, 2018 1 2 3 4 5 6 7 8 9 10 11 12 # # .. _ModalAnalysis: # # ========================================== # Modal analysis of an elastic structure # ========================================== # # ------------- # Introduction # ------------- # # This program performs a dynamic modal analysis of an elastic cantilever beam  Jeremy BLEYER committed Jul 06, 2018 13 14 # represented by a 3D solid continuum. The eigenmodes are computed using the # **SLEPcEigensolver** and compared against an analytical solution of beam theory.  Jeremy BLEYER committed May 28, 2018 15 16 17 18 19 20 21 22 # 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 committed Jul 06, 2018 23 24 # 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 committed May 28, 2018 25 26 27 28 29 # # --------------- # Implementation # --------------- #  Jeremy BLEYER committed Jul 06, 2018 30 # After importing the relevant modules, the geometry of a beam of length :math:L=20  Jeremy BLEYER committed May 28, 2018 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 # 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 committed Jul 06, 2018 46 # take :math:\nu=0) and we also introduce the material density :math:\rho for  Jeremy BLEYER committed May 28, 2018 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 # 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 committed Jul 06, 2018 62 # Standard FunctionSpace is defined and boundary conditions correspond to a  Jeremy BLEYER committed May 28, 2018 63 64 65 66 67 68 69 70 71 72 73 74 75 # 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 committed Jul 06, 2018 76 # The system stiffness matrix :math:[K] and mass matrix :math:[M] are  Jeremy BLEYER committed May 28, 2018 77 78 79 80 81 82 83 84 85 86 87 88 # 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 committed Jul 06, 2018 89 # Matrices :math:[K] and :math:[M] are first defined as PETSc Matrix and  Jeremy BLEYER committed May 28, 2018 90 91 92 # 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 committed Jul 06, 2018 93 # vector have been introduced to call this function).  Jeremy BLEYER committed May 28, 2018 94 95 # #  Jeremy BLEYER committed Jul 06, 2018 96 # Modal dynamic analysis consists in solving the following generalized  Jeremy BLEYER committed May 28, 2018 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 # 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 committed Jul 06, 2018 124   Jeremy BLEYER committed May 28, 2018 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 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 committed Jul 06, 2018 144   Jeremy BLEYER committed May 28, 2018 145 146  # 3D eigenfrequency freq_3D = sqrt(r)/2/pi  Jeremy BLEYER committed Jul 06, 2018 147   Jeremy BLEYER committed May 28, 2018 148 149 150 151 152 153  # 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 committed Jul 06, 2018 154   Jeremy BLEYER committed May 28, 2018 155 156 157 158  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 committed Jul 06, 2018 159  eigenmode.vector()[:] = rx  Jeremy BLEYER committed May 28, 2018 160 161 162  # 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 committed Jul 06, 2018 163 # where :math:S=BH is the beam section, :math:I the bending inertia and  Jeremy BLEYER committed May 28, 2018 164 165 166 167 168 169 170 171 172 173 174 175 # :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 committed Jul 06, 2018 176 #  Jeremy BLEYER committed May 28, 2018 177 # With Nx=400, we obtain the following comparison between the FE eigenfrequencies  Jeremy BLEYER committed Jul 06, 2018 178 # and the beam theory eigenfrequencies :  Jeremy BLEYER committed May 28, 2018 179 180 181 182 183 184 185 186 # # # ===== ============= ================= # Mode Eigenfrequencies # ----- -------------------------------- # # Solid FE [Hz] Beam theory [Hz] # ===== ============= ================= # 1 2.04991 2.01925  Jeremy BLEYER committed Jul 06, 2018 187 # 2 4.04854 4.03850  Jeremy BLEYER committed May 28, 2018 188 # 3 12.81504 12.65443  Jeremy BLEYER committed Jul 06, 2018 189 # 4 25.12717 25.30886  Jeremy BLEYER committed May 28, 2018 190 # 5 35.74168 35.43277  Jeremy BLEYER committed Jul 06, 2018 191 # 6 66.94816 70.86554  Jeremy BLEYER committed May 28, 2018 192 193 194 # ===== ============= ================= # #