Rotational Constraint at the center of a beam

Hi Guys,

I am trying to compute mode shapes for a 2D wedge (in XY plane) with U_x=0 and UR_z=0 (rotation about the z-axis) constraints at the center.

When I tried for the left and right clamped constraints it was working correctly. But for the center constraint, I am getting incorrect results

from dolfin import *
from ufl import nabla_div, nabla_grad
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
import json
from mshr import Polygon, generate_mesh
from scipy.optimize import root
from math import cos, cosh
import sys
import argparse
from functools import partial
import pandas as pd
import meshio as mio
from stl import mesh as mstl

#======================================================================================
# General Functions
#======================================================================================

def int_zfill(index, total_len): 
    from math import floor, log10
    return str(index).zfill(1+floor(log10(total_len)))

print(int_zfill(2,271))
#======================================================================================
# File locations
#======================================================================================


main_dir = Path(".")
out_dir = main_dir / "CaseElasWedge_out"

properties_file = main_dir / "properties.json"
mesh_file = main_dir / "mesh.xml"
plot_dir = main_dir / "plots"
fdebug_file = out_dir / "FloatingAce_mkbound_0.csv"
fmotion_file = out_dir / "floatinginfo/FloatingMotion_mk11.csv"
fforce_file = out_dir / "forces/_FloatingForce.csv"
stress_dir = out_dir / "stress"

nonlin = False
Path.cwd()

mesh_file

#======================================================================================
# Meshing
#======================================================================================

# Paths to the material property folders

plot_dir.mkdir(parents=True, exist_ok=True)

# importing properties
with properties_file.open(mode="r") as pf:
    properties = json.load(pf)

wedge_length = 0.3
dim = 2
N_eig = 16   # number of eigenvalues

wedge_thickness = properties["thickness"]
youngs_modulus = Constant(properties["young's modulus"])
poisson_ratio = Constant(properties["poisson's ratio"])
density = Constant(properties["density"])
beta_deg = properties["beta"]
beta_rad = np.deg2rad(beta_deg)
cosb = np.cos(beta_rad)
sinb = np.sin(beta_rad)
gravity = 9.81
    
domain = Polygon([
    Point(-wedge_length * cosb, wedge_length * sinb),
    Point(0.0, 0.0),
    Point(wedge_length * cosb, wedge_length * sinb),
    Point(wedge_length * cosb, wedge_length * sinb + wedge_thickness * cosb),
    Point(0.0, wedge_thickness),
    Point(-wedge_length * cosb, wedge_length * sinb + wedge_thickness * cosb),
])

mesh = generate_mesh(domain,100)
mesh

strain_gauges = [
    Point(0.03 * cosb, 0.03 * sinb + wedge_thickness / cosb),
    Point(0.12 * cosb, 0.12 * sinb + wedge_thickness / cosb)
]
strain_gauges

File("mesh.pvd") << mesh
File("mesh.xml") << mesh

msh = mio.read("mesh000000.vtu")
mio.write("mesh_xy.stl",msh)

stl_mesh = mstl.Mesh.from_file('mesh_xy.stl')
stl_mesh.rotate([1.0, 0.0, 0.0], np.radians(270))
stl_mesh.save('mesh_xz.stl')
[
    (0.0, 0.0),
    (wedge_length * cosb, wedge_length * sinb),
    (wedge_length * cosb, wedge_length * sinb + wedge_thickness * cosb),
    (0.0, wedge_thickness),
    (-wedge_length * cosb, wedge_length * sinb + wedge_thickness * cosb),
    (-wedge_length * cosb, wedge_length * sinb),
    properties
]

#======================================================================================
# Strain and Stress Functions
#======================================================================================

V = VectorFunctionSpace(mesh, 'Lagrange', degree=1)
v2d = vertex_to_dof_map(V).reshape((-1, dim))
d2v = dof_to_vertex_map(V)[range(0, len(dof_to_vertex_map(V)), dim)]/dim
u_ = TrialFunction(V)
du = TestFunction(V)

# Lame coefficient for constitutive relation
mu = Constant( youngs_modulus / (2.0*(1+poisson_ratio)) )
lmbda = Constant(youngs_modulus * poisson_ratio/((1+poisson_ratio)*(1-2*poisson_ratio)))

def eps(v):
    return 0.5*(nabla_grad(v) + nabla_grad(v).T)

def green_eps(v):
    I = Identity(dim)
    F = I + nabla_grad(v)
    return 0.5*(F.T*F - I)

def sigma(eps,v):
    return 2.0*mu*eps(v) + lmbda*tr(eps(v))*Identity(dim)

strain = eps
stress = partial(sigma,strain)


mesh.coordinates()

#======================================================================================
# Boundary Conditions
#======================================================================================

def left(x, on_boundary):
    return on_boundary and near(x[0], -wedge_length * cosb)

def right(x, on_boundary):
    return on_boundary and near(x[0], wedge_length * cosb)

def center(x, on_boundary):
    return on_boundary and near(x[0], 0)

# https://fenicsproject.org/qa/9617/2d-boundary-condition/
# bc = DirichletBC(V, Constant((0,0)), left)
# bc = DirichletBC(V, Constant((0,0)), right)
bc = DirichletBC(V.sub(0), Constant(0), center)

# bc = None

#======================================================================================
# Compute solution
#======================================================================================

# u = TrialFunction(V)
# d = u.geometric_dimension() # space dimension
# v = TestFunction(V)
# f = Constant((0,-density*gravity))
# T = Constant((0,0))
# a = inner(stress(u), strain(v))*dx
# L = dot(f, v)*dx + dot(T, v)*ds

# u = Function(V, name="displacement")
# solve(a == L, u, bc)
# with XDMFFile(f"{main_dir}/deflection.xdmf") as file_results:
#     file_results.parameters["flush_output"] = True
#     file_results.parameters["functions_share_mesh"] = True
#     file_results.write(u, 0.)
# plot(u)

#======================================================================================
# Compute Mode Shapes and frequencies
#======================================================================================
k_form = inner(stress(du),strain(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 = density*dot(du,u_)*dx
M = PETScMatrix()
assemble(m_form, tensor=M)
# https://fenicsproject.discourse.group/t/modal-analysis-of-a-beam-with-plain-strain/3952
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.0

print("Computing {} first eigenvalues...".format(N_eig))
eigensolver.solve(N_eig)
falpha = lambda x: cos(x)*cosh(x)+1
alpha = lambda n: root(falpha, (2*n+1)*pi/2.)['x'][0]
normal_M =  []
normal_K =  []
lamba = []
Modes = []
Modes_vectors = []
file_results = XDMFFile(f"{main_dir}/modes.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True

format_len = np.ceil(np.log10(N_eig)).astype(int)
print("id - freq(Hz) - Ang Freq(rad/s) - Lambda")
for i in range(N_eig):
    # Extract eigenpair
    r, c, rx, cx = eigensolver.get_eigenpair(i)

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

    print(f"{i+1} - {freq_3D:.3f} - {2*np.pi*freq_3D:.3f} - {r:.3f}", )

    # Initialize function and assign eigenvector
    eigenmode = Function(V,name=f"Eigenvector {int_zfill(i+1, N_eig)}")
    eigenmode.vector()[:] = rx
    Modes.append(rx[:].reshape((dim,int(rx[:].shape[0]/dim)),order='F').T)
    Modes_vectors.append(eigenmode)
    normal_M.append(np.dot(rx[:],M*rx[:]))
    normal_K.append(np.dot(rx[:],K*rx[:]))


    # Function can be plotted as usual:
    rx_func = Function(V)
    rx_func.vector()[:] = rx
    file_results.write(eigenmode, 0.)

Can anyone tell me, how to apply U_x=0 and UR_z=0 constraint at the center of the wedge?