Hi all,
I am using mixed-elements for doing modal analysis. I want to extend it to compute the modal participation factor. The geo file is attached here: geometry - Google Drive . The MWE is-
from dolfin import *
import numpy as np
import json
mesh = Mesh()
with XDMFFile("mesh/tetra.xdmf") as infile:
infile.read(mesh)
# Function space, trial and test functions
elem_P = VectorElement('Lagrange', mesh.ufl_cell(), 2)
elem_R = VectorElement('Real', mesh.ufl_cell(), 0)
elem = MixedElement([elem_P, elem_R, elem_R])
V = FunctionSpace(mesh, elem)
u_, U_, C_ = TrialFunctions(V)
du, dU, dC = TestFunctions(V)
# Material properties, constitutive law
E, nu = 200e9, 0.3
rho = Constant(7850)
mu = Constant(E/2./(1+nu))
lmbda = Constant(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)
# Boundary conditions
support1 = CompiledSubDomain("near(x[0], 1.0, tol) && on_boundary", tol=1e-14)
support2 = CompiledSubDomain("near(x[0], 0.0, tol) && on_boundary", tol=1e-14)
support3 = CompiledSubDomain("near(x[1], 1.0, tol) && on_boundary", tol=1e-14)
support4 = CompiledSubDomain("near(x[1], 0.0, tol) && on_boundary", tol=1e-14)
bc_l = DirichletBC(V.sub(0), Constant((0.,0.,0.)), support1)
bc_r = DirichletBC(V.sub(0), Constant((0.,0.,0.)), support2)
bc_t = DirichletBC(V.sub(0), Constant((0.,0.,0.)), support3)
bc_b = DirichletBC(V.sub(0), Constant((0.,0.,0.)), support4)
bc=[bc_l,bc_r,bc_t,bc_b]
# Surface markers
force_face = CompiledSubDomain("(x[0] - 0.41478)*(x[0] - 0.41478) + (x[1] - 0.46787)*(x[1] - 0.46787) < 0.16*0.16")
line_mf = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
line_mf.set_all(0)
force_face.mark(line_mf, 2)
# Parameters of the point mass
hole_radius = 0.15
plate_thk = 0.1
mass = 1000 # kg
area = 2*np.pi*hole_radius*plate_thk
M = Constant(mass)
A = Constant(area)
# Mass and stiffness forms
ds = Measure('ds', mesh, subdomain_data=line_mf)
k_form = inner(sigma(du),eps(u_))*dx - (inner(C_, du-dU) + inner(u_-U_, dC))*ds(2)
m_form = rho*inner(u_,du)*dx + M/A*inner(U_,du)*ds(2)
K = PETScMatrix()
assemble(k_form, tensor=K)
M = PETScMatrix()
assemble(m_form, tensor=M)
for bcs in bc:
bcs.zero(M)
bcs.apply(K)
# Solution
eigensolver = SLEPcEigenSolver(K, M)
eigensolver.parameters['problem_type'] = 'gen_hermitian'
eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
eigensolver.parameters['spectral_shift'] = 0.
n_eig = 5 # number of eigenvalues
print("Computing {} first eigenvalues...".format(N_eig))
eigensolver.solve(n_eig)
# 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
eigenmodes = []
eig_vals = []
eig_v = []
# Extraction
for i in range(N_eig):
# Extract eigenpair
r, c, rx, cx = eigensolver.get_eigenpair(i)
# 3D eigenfrequency
freq_3D = np.sign(r)*np.sqrt(np.abs(r))/2/pi
print("Solid FE ({0:d}): {1:8.5f} [Hz]".format(i+1, freq_3D))
# Initialize function and assign eigenvector
eigenmode = Function(V,name="Eigenvector "+str(i))
eigenmode.vector()[:] = rx
file_results.write(eigenmode.split()[0], 0)
eigenmodes.append(rx)
eig_v.append(eigenmode)
eig_vals.append(freq_3D)
V_new = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
u = Function(V_new, name="Unit displacement")
u.interpolate(Constant((0, 1, 0)))
combined_mass = 0
for i, xi in enumerate(eig_v):
qi = assemble(action(action(m_form, u), xi))
mi = assemble(action(action(m_form, xi), xi))
meff_i = (qi / mi) ** 2
total_mass = assemble(rho * dx(domain=mesh))
print("-" * 50)
print("Mode {}:".format(i + 1))
print(" Modal participation factor: {:}".format(qi))
print(" Modal mass: {:}".format(mi))
print(" Effective mass: {:}".format(meff_i))
The demo for calculation of modal participation factor can be found here