As @nate pointed out, this mode is associated with the DOFs on the surface where the Dirichlet BC is applied. FEniCS does not remove these rows and columns from K and M, and instead zeroes the corresponding rows and columns and makes an entry on the diagonal that depends on the number of connected cells. For the set of properties chosen in the tutorial, the frequency of these modes is higher than the beam modes, so the eigensolver recovers the beam modes. Using your properties for aluminum, these spurious modes become lower in frequency than the beam modes.
To avoid seeing these modes, you can set the diagonal entries of K and M such that the spurious modes have a higher frequency than the frequency range of interest, as suggested here. One possible choice is below, where I calculate the highest eigenfrequency of the system and set the diagonal entries so that the spurious modes have a frequency equal to the highest frequency. You can see that the beam modes are now correctly recovered.
from dolfin import *
import numpy as np
L, B, H = 1, 0.1, 0.04
Nx = 40
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)
E, nu = Constant(69.0E9), Constant(0.3)
rho = Constant(2712.0)
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)
V = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
u_ = TrialFunction(V)
du = TestFunction(V)
def left(x, on_boundary):
return near(x[0],0.)
bc = DirichletBC(V, Constant((0.,0.,0.)), left)
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_system(m_form, l_form, bc, A_tensor=M, b_tensor=b)
# Calculate the largest eigenvalue
eigensolver = SLEPcEigenSolver(K, M)
eigensolver.parameters['problem_type'] = 'gen_hermitian'
eigensolver.solve(1)
eig_max, _, _, _ = eigensolver.get_eigenpair(0)
# Set diagonal entries of K and M so that the frequency of Dirichlet BC
# eigenmodes is equal to the highest eigenfrequency
bc_idx = np.array([i for i in bc.get_boundary_values().keys()], dtype=np.int32)
K.get_diagonal(b)
b[bc_idx] = np.sqrt(eig_max)
K.set_diagonal(b)
M.get_diagonal(b)
b[bc_idx] = 1/np.sqrt(eig_max)
M.set_diagonal(b)
# Solve the beam eigenmodes
eigensolver = SLEPcEigenSolver(K, M)
eigensolver.parameters['problem_type'] = 'gen_hermitian'
eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
eigensolver.parameters['spectral_shift'] = 0.
N_eig = 10 # number of eigenvalues
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]
# Set up file for exporting results
file_results = XDMFFile("mod_syn_evp-commet.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
eigenmodes = []
# 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(float(E)*I_bend/(float(rho)*B*H*L**4))/2/pi
eigenmode = Function(V,name="Eigenvector "+str(i))
eigenmode.vector()[:] = rx
eigenmodes.append(eigenmode)
file_results.write(eigenmode, 0)
