Mode shapes of Free-Free Beam in FENICS

Linear Algebra
1

Hi guys,

Can anyone tell me how to compute and export mode shapes (D(x,y,z)) and normalized mass (D^TMD) and normalized stiffness (D^TKD) matrices of a Free-Free beam (or any object with no boundary conditions) from FENICS.

I am currently using the following example code

    from dolfin import *
    import numpy as np

    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)

    E, nu = Constant(1e5), Constant(0.)
    rho = Constant(1e-3)
    # 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)
    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['spectral_transform'] = 'shift-and-invert'
    eigensolver.parameters['spectral_shift'] = 0.

    N_eig = 6   # 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("modal_analysis.xdmf")
    file_results.parameters["flush_output"] = True
    file_results.parameters["functions_share_mesh"] = True


    # Extraction
    normal_M = np.zeros(N_eig+2)
    normal_K = np.zeros(N_eig+2)
    lamba = np.zeros(N_eig+2)
    Modes = []
    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+2] = r

        # 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"Solid FE: {freq_beam} [Hz]")

        # Initialize function and assign eigenvector
        eigenmode = Function(V,name="Eigenvector "+str(i))
        eigenmode.vector()[:] = rx
        Modes.append(rx[:])
        normal_M[i+2] = np.dot(rx[:],M*rx[:])
        normal_K[i+2] = np.dot(rx[:],K*rx[:])

Thank you

2

I think I resolved the issue.

I should change

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)

to


# 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) # Change BC to None

m_form = rho*dot(du,u_)*dx
M = PETScMatrix()
assemble(m_form, tensor=M)