Consideration of point mass in modal analysis code

8

Hi @Shubham_Saurabh,

You should look at the Ansys manual to see how point masses are implemented and how the equations are modified. I’m not familiar with modal analysis in Ansys so I don’t know how the frequencies are obtained.

With that being said, in the results you present from Ansys, it looks like the addition of a point mass has a much stronger effect on the frequency of the lowest-frequency mode than any of the other modes. To me, this suggests that the point mass is implemented via an additional vector-valued dof that is constrained to the average displacement of the hole surface. The first-order plate mode has a nonzero average displacement, while the second-order plate mode shapes are such that the average displacement of the hole surface is nearly 0; thus the point mass participates much more strongly in the first mode than in the higher-order modes.

To implement this, define the displacement \mathbf{U} of the point mass m. The point mass is subjected to a force \mathbf{F} exerted upon it by the surface of the hole, so the equation of motion for the point mass is given by

\mathbf{F} = m \ddot{\mathbf{U}} = -\omega^2 m \mathbf{U}

The reaction of the point mass causes the hole surface to be subjected to a uniform traction \mathbf{T} such that

\mathbf{T} = - \frac{\mathbf{F}}{|\partial \Omega_1|} = \frac{\omega^2 m \mathbf{U}}{|\partial \Omega_1|}

where |\partial \Omega_1| is the area of the hole surface. The weak form of the equations of motion for the plate is (where quantities with a hat “\hat{}” are test quantities):

-\int_{\partial \Omega_1} {\mathbf{T} \cdot \hat{\mathbf{u}} \,\mathrm{d}s} + \int_{\Omega} { \mathbf{\sigma} : \hat{\mathbf{\epsilon}} \, \mathrm{d}x} = \omega^2 \int_{\Omega} {\rho \mathbf{u} \cdot \hat{\mathbf{u}} \, \mathrm{d}x} \\ - \frac{m \omega^2}{|\partial \Omega_1|} \int_{\partial \Omega_1} {\mathbf{U} \cdot \hat{\mathbf{u}} \,\mathrm{d}s} + \int_{\Omega} { \mathbf{\sigma} : \hat{\mathbf{\epsilon}} \, \mathrm{d}x} = \omega^2 \int_{\Omega} {\rho \mathbf{u} \cdot \hat{\mathbf{u}} \, \mathrm{d}x} \\ \int_{\Omega} { \mathbf{\sigma} : \hat{\mathbf{\epsilon}} \, \mathrm{d}x} = \omega^2 \Biggl\{ \int_{\Omega} {\rho \mathbf{u} \cdot \hat{\mathbf{u}} \, \mathrm{d}x} + \frac{m}{|\partial \Omega_1|} \int_{\partial \Omega_1} {\mathbf{U} \cdot \hat{\mathbf{u}} \,\mathrm{d}s} \Biggr\}

The displacement of the point mass can be constrained to the average displacement of the hole surface via a Lagrange multiplier approach, i.e. by adding the terms

- \mathbf{C} \cdot \int_{\partial \Omega_1} {(\hat{\mathbf{u}} - \hat{\mathbf{U}}) \, \mathrm{d}s} - \hat{\mathbf{C}} \cdot \int_{\partial \Omega_1} {(\mathbf{u} - \mathbf{U}) \, \mathrm{d}s}

to the functional, where \mathbf{C} is the Lagrange multiplier. Thus we have for the K form:

\int_{\Omega} { \mathbf{\sigma} : \hat{\mathbf{\epsilon}} \, \mathrm{d}x} - \mathbf{C} \cdot \int_{\partial \Omega_1} {(\hat{\mathbf{u}} - \hat{\mathbf{U}}) \, \mathrm{d}s} - \hat{\mathbf{C}} \cdot \int_{\partial \Omega_1} {(\mathbf{u} - \mathbf{U}) \, \mathrm{d}s}

and for the M form:

\int_{\Omega} {\rho \mathbf{u} \cdot \hat{\mathbf{u}} \, \mathrm{d}x} + \frac{m}{|\partial \Omega_1|} \int_{\partial \Omega_1} {\mathbf{U} \cdot \hat{\mathbf{u}} \,\mathrm{d}s}

With this approach, I obtain:

Solid FE (1): 286.95752 [Hz]
Solid FE (2): 1236.11988 [Hz]
Solid FE (3): 1241.04164 [Hz]
Solid FE (4): 1366.32870 [Hz]
Solid FE (5): 1456.47904 [Hz]

For reference, my code 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

# 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)
5 Likes