Modal analysis of a barge made of shell elements

Linear Algebra
1

Hi guys,
I am trying to perform a modal analysis of a cuboid barge (with stiffeners) made of shell elements. Here is the code

barge
barge736×634 36.8 KB

from dolfin import *
from vedo.dolfin import plot
import numpy as np

mesh = Mesh("../barge.xml")
plot(mesh)

E, nu = Constant(1e5), Constant(0.)
rho = Constant(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)

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)

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)
bc.zero(M)

eigensolver = SLEPcEigenSolver(K, M)
eigensolver.parameters['problem_type'] = 'gen_hermitian'
eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
eigensolver.parameters['spectral_shift'] = 0.
eigensolver.parameters["spectrum"] = "largest magnitude"
eigensolver.parameters['solver'] = "power"
print(eigensolver.parameters.str(True))

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
for i in range(N_eig):
    # Extract eigenpair
    r, c, rx, cx = eigensolver.get_eigenpair(i)

    # 3D eigenfrequency
    freq_3D = sqrt(r)/2/pi

    print("Solid FE: {0:8.5f} [Hz]   ".format(freq_3D))

    # Initialize function and assign eigenvector
    eigenmode = Function(V,name="Eigenvector "+str(i))
    eigenmode.vector()[:] = rx

But I am getting the following error

Computing 6 first eigenvalues...
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-9-62814374ba19> in <module>
     17 for i in range(N_eig):
     18     # Extract eigenpair
---> 19     r, c, rx, cx = eigensolver.get_eigenpair(i)
     20 
     21     # 3D eigenfrequency

RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to extract eigenpair from SLEPc eigenvalue solver.
*** Reason:  Requested eigenpair (0) has not been computed.
*** Where:   This error was encountered inside SLEPcEigenSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  c5b9b269f4a6455a739109e3a66e036b5b8412f5
*** -------------------------------------------------------------------------

Can anyone tell how to resolve this issue ?

2

Hi,
I don’t understand. Your mesh is made of tetrahedra or of 2D triangular shell elements ? In the latter case, you need to implement a shell model, the above code implements 3D elasticity which will not be suited to your needs.
However, implementing a locking-free shell model with a non-manifold mesh is not so straightforward with FEniCS.

3

Thanks for the reply.

I am currently trying to implement it in the fenics-shell. But even there I am getting error.

Can you tell me how to perform a modal analysis for the above barge?

Thank you