Hi guys,
I am performing modal analysis of a 2D rectangular beam and I am getting what looks like rigid body modes at non-zero frequencies and repeated modes at higher frequencies.
Can anyone tell me how to resolve this issue?
# To add a new cell, type '# %%'
# To add a new markdown cell, type '# %% [markdown]'
# %%
#
# .. # gedit: set fileencoding=utf8 :
#
# .. 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>
#
# .. _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 dolfin import *
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
Path("plots").mkdir(parents=True, exist_ok=True)
dim = 2
N_eig = 16 # number of eigenvalues
if dim ==2:
L, H = 1.0, 0.02
point0 = Point(-0.4,1.02)
# point0 = Point(-0.5,-0.01)
Nx = 600
Nz = int(H/L*Nx)
print(Nz)
mesh = RectangleMesh(point0,Point(point0[0]+L,point0[1]+H), Nx, Nz)
elif dim == 3:
L, B, H = 1.0, 0.02, 0.02
point0 = Point(-0.4,0.05,1.02)
Nx = 200
Ny = int(B/L*Nx)+1
Nz = int(H/L*Nx)+1
mesh = BoxMesh(point0,Point(point0[0]+L,point0[1]+B,point0[2]+H), Nx, Ny, Nz)
else:
raise ValueError("The dimension is neither 2D or 3D")
mesh
# %%
# 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 = Constant(0.025e9), Constant(0.3)
rho = Constant(500.0)
# E, nu = 1e5, 0.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)
v2d = vertex_to_dof_map(V).reshape((-1, dim))
d2v = dof_to_vertex_map(V)[range(0, len(dof_to_vertex_map(V)), dim)]/dim
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, None, 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['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). A
# spectral transform is therefore 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)::
print("Computing {} first eigenvalues...".format(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]
# Extraction
normal_M = np.zeros(N_eig)
normal_K = np.zeros(N_eig)
lamba = np.zeros(N_eig)
# %%
Modes = []
file_results = XDMFFile("xdmf/elasticity_results1.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
# file_results.write(mesh)
for i in range(N_eig):
# Extract eigenpair
r, c, rx, cx = eigensolver.get_eigenpair(i)
# 3D eigenfrequency
freq_3D = sqrt(r)/2/pi
lamba[i] = r
# freq[i] = sqrt(lamba[i])/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(f"{freq_3D:.3f}")
# Initialize function and assign eigenvector
eigenmode = Function(V,name="Eigenvector "+str(i))
eigenmode.vector()[:] = rx
Modes.append(rx[:].reshape((dim,int(rx[:].shape[0]/dim)),order='F').T)
normal_M[i] = np.dot(rx[:],M*rx[:])
normal_K[i] = np.dot(rx[:],K*rx[:])
# Function can be plotted as usual:
rx_func = Function(V)
rx_func.vector()[:] = rx
file_results.write(rx_func, 0.)
plt.figure(figsize=(20,5))
plot(rx_func)
plt.savefig(f"plots/pmode_{i}_freq_{freq_3D:6f}.png", format="png", dpi=600)
# plt.autoscale()
plt.close()
Modes = np.array(Modes)
freq=np.sqrt(np.abs(lamba))/(2*np.pi)
Modes.shape
# %%
arr = rx_func.vector()[:]
coor = mesh.coordinates()
d2v.shape, v2d.shape
N_eig
# %%
v2d.shape, d2v.shape
# %%
# values = []
# for i, dum in enumerate(coor):
# values.append([arr[v2d[dim*i]],arr[v2d[dim*i+1]]])
# values = np.array(values)
# values
# %%
V.tabulate_dof_coordinates()
# %%
rx_func.vector()[:].reshape((dim,int(rx_func.vector()[:].shape[0]/dim)),order='F').T
# %%
Modes.shape
# %%
Modes_new = np.array([
Modes[:,i,:].reshape(N_eig*dim) for i in range(Modes.shape[1])
])
Modes_new.shape
# %%
mesh.coordinates()[int(d2v[1]),:]
# %%
nodenum = mesh.coordinates().shape[0]
coords = np.zeros((nodenum, dim))
for i in range(mesh.coordinates().shape[0]):
coords[i,:] = mesh.coordinates()[int(d2v[i]),:]
coords
# %%
for i in range(20):
print(f"coord x-{coords[i,0]:4f}, y-{coords[i,1]:4f} \t mode[0] x-{Modes[0,i,0]:4f}, y-{Modes[0,i,1]:4f}")
# %%
if dim==2:
rigid_dof = 3
elif dim==3:
rigifd_dof = 6
Assembled_mode_nodes = np.hstack((mesh.coordinates(),Modes_new[:,rigid_dof*dim:]))
Assembled_mass_stiff = np.transpose([normal_M,normal_K])[rigid_dof:,:]
Modes_new[0,rigid_dof*dim:]
Modes_new[0,:]
Assembled_mode_nodes.shape
# %%
if dim==2:
nodes_2d_to_3d = coords[:,0]
nodes_2d_to_3d = np.vstack((nodes_2d_to_3d, 0.05*np.ones(nodenum)))
nodes_2d_to_3d = np.vstack((nodes_2d_to_3d, coords[:,1]))
for i in range(1,N_eig-rigid_dof+1):
nodes_2d_to_3d = np.vstack((nodes_2d_to_3d, Assembled_mode_nodes[:,2*i]))
nodes_2d_to_3d = np.vstack((nodes_2d_to_3d, np.zeros(nodenum)))
nodes_2d_to_3d = np.vstack((nodes_2d_to_3d, Assembled_mode_nodes[:,2*i+1]))
print(i, 2*i, 2*i+1, int(Assembled_mode_nodes.shape[1]/dim))
print("nodes_2d_to_3d.shape", nodes_2d_to_3d.shape)
nodes_2d_to_3d = nodes_2d_to_3d.T
for i in range(1,N_eig-rigid_dof+1):
print(i)
np.savetxt(f'modeshapes2d_to_3d_mode{i-1}.csv',nodes_2d_to_3d[:,:3*i], delimiter=",",fmt='%.8f')
nodes_2d_to_3d.shape
# %%
np.savetxt(f'mass_stiff{dim}D.csv',Assembled_mass_stiff,delimiter=",",fmt='%.8f')
np.savetxt(f'modeshapes{dim}D.csv',Assembled_mode_nodes, delimiter=",",fmt='%.8f')
# %%
fig, ax = plt.subplots(N_eig,1,figsize=(50,50))
for k, freqi in enumerate(freq):
X = coords[:,0]
Y = coords[:,1]
U = Modes_new[:,2*k]
V1 = Modes_new[:,2*k+1]
color_array = np.sqrt(U**2 + V1**2)
# ax1.set_ylim(-0.02, 0.02)
# ax1.set_title(f'Mode Displacement vector plots for frequency = {freqi:6f} [Hz]')
Q = ax[k].quiver(X, Y, U, V1, color_array, units='inches',color='r')
qk = ax[k].quiverkey(Q, 0.5, 0.5, 0.1, "", labelpos='E',
coordinates='figure')
ax[k].get_xaxis().set_visible(False)
ax[k].set_ylabel(f"{freq[k]:.1f}", rotation=0, fontsize=20, labelpad=20)
plt.tight_layout()
plt.savefig(f"plots/modes.png", format="png", dpi=300)
plt.show()
# %%
N_eig-rigid_dof, nodenum
# %%
lamba_n = np.double(lamba[3:])
lamba_n
