Boundary conditions for circular domain in eigenvalue problem


61
views
0
5 weeks ago by
Daniel  
I am trying to compute the eigenmode of a circular structure with Dirichlet boundary conditions on the circumference of the circular domain. When I run the program without bcs I get the expected result like this mode shape



However, as soon as I impose boundary conditions to the problem using a function for defining the boundary I get strange results like this one



So I suppose that I do something wrong when defining the boundary in my script.

import dolfin
import fenics as fe
import mshr

# Numerical parameters
n_eigen = 30

# Mesh paramters
n_mesh = 10

# Geometry parameters
r = 250
h = 0.01

# Material parameters
E = 170e9
nu = 0.28
rho = 2330.0

# Lame parameters
mu = E / (2 * (1 + nu))
lmbda = E * nu / ((1 + nu) * (1 - 2 * nu))

# Definitions of stress and strain
def epsilon(u):
    return fe.sym(fe.grad(u))

def sigma(u):
    n_dim = u.geometric_dimension()
    return 2 * mu * epsilon(u) + lmbda * fe.tr(epsilon(u)) * fe.Identity(n_dim)

mesh = fe.Mesh('pygmsh_circle.xml')

# # Generate geometry
# domain = mshr.Cylinder(fe.Point(0, 0, h), fe.Point(0, 0, 0), r, r)
# mesh = mshr.generate_mesh(domain, n_mesh)

f = fe.File('mesh.pvd')
f << mesh

# Create test and trial functions
V = fe.VectorFunctionSpace(mesh, 'Lagrange', degree=1)
u = fe.TrialFunction(V)
v = fe.TestFunction(V)

# Define one clamped end as boundary condition
def boundary(x, on_boundary):
    if fe.near(x[0]**2 + x[1]**2, r**2, 10) and on_boundary:
        # print(x[0]**2 + x[1]**2)
        print(x[0])
        print(x[1])
        # print(x[2])
        print('******************')
        return True
    else:
        return False
bc = fe.DirichletBC(V, fe.Constant((0.0, 0.0, 0.0)), boundary)

# Assemble stiffness matrix
k_form = fe.inner(sigma(u), epsilon(v)) * fe.dx
l_form = fe.Constant(1) * u[0] * fe.dx
K = fe.PETScMatrix()
b_tmp = fe.PETScVector()
fe.assemble_system(k_form, l_form, bc, A_tensor=K, b_tensor=b_tmp)
#fe.assemble_system(k_form, l_form, None, A_tensor=K, b_tensor=b_tmp)

# Assemble mass matrix
m_form = rho * fe.dot(u, v) * fe.dx
M = fe.PETScMatrix()
fe.assemble(m_form, tensor=M)

# Setup SLEPc eigensolver
eigensolver = fe.SLEPcEigenSolver(K, M)
eigensolver.parameters['problem_type'] = 'gen_hermitian'
eigensolver.parameters['spectrum'] = 'smallest real'
eigensolver.parameters['spectral_transform'] = 'shift-and-invert'
eigensolver.parameters['spectral_shift'] = 0.0

# Compute eigenvalues
eigensolver.solve(n_eigen)

 # Eigenvalue extraction loop
f = fe.File('modes.pvd')
for i_eigen in range(n_eigen):   
    # Pick on eigenvalue
    r, c, rx, cx = eigensolver.get_eigenpair(i_eigen)

    u_ex = fe.Function(V)
    u_ex.vector()[:] = rx
    f << u_ex​

The script for generating the mesh looks like that:

import pygmsh
import meshio

r = 250
h = 10
lcar = 10

geom = pygmsh.opencascade.Geometry()
c1 = geom.add_cylinder([0, 0.0, 0.0], [0, 0, h], r, char_length=lcar)

points, cells, _, _, _ = pygmsh.generate_mesh(geom, geo_filename="tmp.geo")
m = meshio.Mesh(points, cells)
meshio.write("pygmsh_circle.vtu", m)
meshio.write("pygmsh_circle.xml", m, file_format='dolfin-xml')

I suppose it must be a stupid mistake in

def boundary(x, on_boundary):
    if fe.near(x[0]**2 + x[1]**2, r**2, 10) and on_boundary:
        # print(x[0]**2 + x[1]**2)
        print(x[0])
        print(x[1])
        # print(x[2])
        print('******************')
        return True
    else:
        return False

where I also added some print statements for debugging. Can someone help me with this problem or point me into the right direction? Thanks in advance
Community: FEniCS Project
I edited the original post because I realized the problem is not suitable for a 2d treatment.
written 5 weeks ago by Daniel  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »