I am importing a bdf file from commercial package(Ansys) to Gmsh and converting it to XDMF using meshio and importing that XDMF file into fenics. In commercial package(Ansys), I have considered point masses at certain locations. But I am not able to incorporate these point masses in fenics implementation. Can you please guide me regarding this?
Hi,
you can use PointSource to add additional contribution to the mass matrix
Thank you. It worked for point mass at any node. But my problem is slightly different. Consider a 3D plate with holes. In these holes, there were some other geometries. That I am not considering for this study. I am using point mass to represent those geometries. How to account for masses of those geometries.
I assume that you should somehow distribute the mass on the corresponding hole boundary, just adding to the mass form something like:
M/S*dot(u,v)*ds(1)
where M is the corresponding mass and S the are of the boundary of index 1.
1 Like
With application for above syntax, I am getting the same natural frequencies with or without point mass. I am attaching a MWE for the same. The imported geometry is a square plate with hole. 1000 Kg mass is being applied to the internal face of this hole.
from dolfin import *
import numpy as np
import json
import matplotlib.pyplot as plt
import scipy as sp
%matplotlib inline
mesh = Mesh()
with XDMFFile("mesh/tetra.xdmf") as infile:
infile.read(mesh)
V = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
E, nu = (71e9), (0.33) #N/m**2
rho=2820 #Kg/m**3
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)
u_ = TrialFunction(V)
du = TestFunction(V)
support1 = CompiledSubDomain("x[0] == 1")
support2 = CompiledSubDomain("x[0] == 0")
support3 = CompiledSubDomain("x[1] == 1")
support4 = CompiledSubDomain("x[1] == 0")
bc=[bc_l,bc_r,bc_t,bc_b]
force_face = CompiledSubDomain("(x[0] - 0.41478)*(x[0] - 0.41478) + (x[1] - 0.46787)*(x[1] - 0.46787) + (x[2] - 0.05)*(x[2] - 0.05) < 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)
k_form = inner(sigma(du),eps(u_))*dx
m_form = rho*dot(du,u_)*dx + 1000/(2*3.14*0.16*0.01)*dot(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)
eigensolver = SLEPcEigenSolver(K, M)
eigensolver.parameters['problem_type'] = 'gen_hermitian'
eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
eigensolver.parameters['spectral_shift'] = 0.
N_eig = 40 # 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
eig=[]
eig_v = []
eig_l = []
for i in range(N_eig):
# Extract eigenpair
r, c, rx, cx = eigensolver.get_eigenpair(i)
# 3D eigenfrequency
freq_3D = np.sqrt(r)/2/pi
print("Mode {}:".format(i + 1))
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
file_results.write(eigenmode, 0)
eig_l.append(freq_3D)
eig_v.append(eigenmode)
eig.append(rx)
file_results.close()
In both cases, with or without external mass, the natural frequencies are coming as:
Solid FE: 566.11561 [Hz]
Solid FE: 896.82914 [Hz]
Solid FE: 1304.03435 [Hz]
Solid FE: 1548.16488 [Hz]
You need to define a Measure for the surface mass contribution in order to recognize the markers, i.e.:
ds = Measure('ds', mesh, subdomain_data=line_mf)
I attempted to recreate your mesh using a 1\times 1\times 0.1 plate with a cylindrical hole of radius 0.16 located at (0.41478,0.46787), for which I obtain
Solid FE (1): 864.85799 [Hz]
Solid FE (2): 1386.40250 [Hz]
Solid FE (3): 1490.76239 [Hz]
Solid FE (4): 2034.41484 [Hz]
Solid FE (5): 2362.79009 [Hz]
before I add the Measure, and
Solid FE (1): 101.99018 [Hz]
Solid FE (2): 168.69723 [Hz]
Solid FE (3): 233.41455 [Hz]
Solid FE (4): 355.96953 [Hz]
Solid FE (5): 377.74612 [Hz]
after I add the Measure.
I was unable to run your original code (in dolfin 2019.1.0) due to missing definitions of bc_l, bc_r, bc_t, and bc_b. In addition, assuming the hole is cylindrical, you should use
force_face = CompiledSubDomain("(x[0] - 0.41478)*(x[0] - 0.41478)"
+ "+ (x[1] - 0.46787)*(x[1] - 0.46787) < 0.16*0.16")
to mark the surface of the hole, i.e. removing the terms with x[2]. For reference, my final code is below:
from dolfin import *
import numpy as np
import json
import matplotlib.pyplot as plt
import scipy as sp
mesh = Mesh()
with XDMFFile("mesh/tetra.xdmf") as infile:
infile.read(mesh)
V = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
E, nu = (71e9), (0.33) #N/m**2
rho=2820 #Kg/m**3
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)
u_ = TrialFunction(V)
du = TestFunction(V)
support1 = CompiledSubDomain("""abs(x[0]) < DOLFIN_EPS ||
abs(x[0] - 1) < DOLFIN_EPS ||
abs(x[1]) < DOLFIN_EPS ||
abs(x[1] - 1) < DOLFIN_EPS""")
bc=[DirichletBC(V, (0,0,0), support1)]
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)
k_form = inner(sigma(du),eps(u_))*dx
ds = Measure('ds', mesh, subdomain_data=line_mf)
m_form = rho*dot(du,u_)*dx + 1000/(2*3.14*0.16*0.01)*dot(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)
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
eig=[]
eig_v = []
eig_l = []
for i in range(N_eig):
# Extract eigenpair
r, c, rx, cx = eigensolver.get_eigenpair(i)
# 3D eigenfrequency
freq_3D = np.sqrt(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, 0)
eig_l.append(freq_3D)
eig_v.append(eigenmode)
eig.append(rx)
file_results.close()
2 Likes
Thank you for the clarification. I tried this with plate having hole. The first fundamental frequency is close to the frequency obtained from commercial package but other frequencies are not matching. Is there any way with which we can refine the formulation.
For the reference, the first five frequencies obtained from the code are:
Solid FE (1): 282.81 [Hz]
Solid FE (2): 500.242 [Hz]
Solid FE (3): 509.604 [Hz]
Solid FE (4): 927.270 [Hz]
Solid FE (5): 937.056 [Hz]
and the first five frequencies obtained from the commercial package are:
Solid FE (1): 287.55 [Hz]
Solid FE (2): 1239.7 [Hz]
Solid FE (3): 1244.5 [Hz]
Solid FE (4): 1366.7 [Hz]
Solid FE (5): 1453.9 [Hz]
The problem statement is: In a plate of 1 m*1 m with 0.1m thickness there is a hole of radius 0.15m inside it, fixed boundary conditions applied at all fours edges and point mass of 1000 kg is applied at the centre of hole. The link for geometry file has been attached at the end.
The final code that I am using is:
from dolfin import *
import numpy as np
import json
import matplotlib.pyplot as plt
import scipy as sp
mesh = Mesh()
with XDMFFile("mesh/tetra.xdmf") as infile:
infile.read(mesh)
V = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
E, nu = 2e11, 0.3
rho = 7850
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)
u_ = TrialFunction(V)
du = TestFunction(V)
support1 = CompiledSubDomain("x[0] == 1")
support2 = CompiledSubDomain("near(x[0], 0.0, tol) && on_boundary", tol=1e-14)
support3 = CompiledSubDomain("x[1] == 1")
support4 = CompiledSubDomain("x[1] == 0")
support=[support1, support2, support3, support4]
bc_l = DirichletBC(V, Constant((0.,0.,0.)), support1)
bc_r = DirichletBC(V, Constant((0.,0.,0.)), support2)
bc_t = DirichletBC(V, Constant((0.,0.,0.)), support3)
bc_b = DirichletBC(V, Constant((0.,0.,0.)), support4)
bc=[bc_l,bc_r,bc_t,bc_b]
force_face = CompiledSubDomain("(x[0] - 0.41478)*(x[0] - 0.41478) + (x[1] - 0.46787)*(x[1] - 0.46787) + (x[2] - 0.05)*(x[2] - 0.05) < 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)
k_form = inner(sigma(du),eps(u_))*dx
ds = Measure('ds', mesh, subdomain_data=line_mf)
m_form = rho*dot(du,u_)*dx + 1000/(2*np.pi*0.15*0.1)*dot(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)
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
eig=[]
eig_v = []
eig_l = []
for i in range(N_eig):
# Extract eigenpair
r, c, rx, cx = eigensolver.get_eigenpair(i)
# 3D eigenfrequency
freq_3D = np.sqrt(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, 0)
eig_l.append(freq_3D)
eig_v.append(eigenmode)
eig.append(rx)
The geo file is attached here: geometry - Google Drive
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
Thanks a lot for the detailed description.