Dirichlet boundary condition for modal analysis of an elastic structure

1

Hello everyone, I am new to FEniCS and I am struggling with modal analysis. I followed “modal analysis of an elastic structure” at Modal analysis of an elastic structure — Numerical tours of continuum mechanics using FEniCS master documentation . Which worked pretty well. But after I manupulate the original code and use the aluminum material parameters, the first mesh where boundary is not
keep it’s shape. Look very strange. I am really confused, any hint is greatly appreciated.
Following is my code and parawiev screenshots of the mesh.

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(m_form, tensor=M)

eigensolver = SLEPcEigenSolver(K, M)
eigensolver.parameters['problem_type'] = 'gen_hermitian'
eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
eigensolver.parameters['spectral_shift'] = 0.
N_eig = 30   # 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, i)

1
12242×1115 65.9 KB

2

Isn’t this just the eigenmode resulting from the set Dirichlet DoFs?

3

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)

image
image1498×832 179 KB

1 Like