Rotational springs at the boundary for rectangualr plate

I am trying to obtain eigenvalues for mindlin rectangular plate with rotational springs (rotational stiffness) in two opposite boundaries and free edge condition on the other two edges. The edges where rotational stiffness is apllied, all other degrees of freedom are contrained. Only rotation is allowed along the length of the plate but with rotational springs, so that, frequency of the plate is controlled. I wrote the code, however, even after applying the moment boundary conditions (product of rotational stiffness and theta, the eigen frequencies are not showing any difference. Here, is the snippet of my code. Any suggestions???

Spring stiffness applied on the left and the right boundary
class lefty(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0) and on_boundary

class righty(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.074) and on_boundary

exterior facets MeshFunction
left_traction = lefty()
exterior_facet_domain1 = MeshFunction(“size_t”, mesh, mesh.geometry().dim())
exterior_facet_domain1.set_all(0)
left_traction.mark(exterior_facet_domain1, 1)
ds = Measure(“ds”) (subdomain_data = exterior_facet_domain1)

M_left = Expression((‘M_l’), M_l = 100000, degree = 0)

right_traction = righty()
exterior_facet_domain2 = MeshFunction(“size_t”, mesh, mesh.geometry().dim())
exterior_facet_domain2.set_all(0)
right_traction.mark(exterior_facet_domain2, 1)
ds = Measure(“ds”) (subdomain_data = exterior_facet_domain2)

M_right = Expression((‘M_r’), M_r = 100000, degree = 0)

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

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

Simply supported boundary conditions.
bc1 = DirichletBC(V.sub(0), Constant(0.0), left) # V.sub(0) means applying restriction along the x direction

bc2 = DirichletBC(V.sub(0), Constant(0.0), right)

bc = [bc1, bc2]

Some useful functions for implementing generalized constitutive relations are now defined:

def strain2voigt(eps):
return as_vector([eps[0, 0], eps[1, 1], 2*eps[0, 1]])

def voigt2stress(S):
return as_tensor([[S[0], S[2]], [S[2], S[1]]])

def curv(u):
(w, theta) = split(u)
return sym(grad(theta))

def shear_strain(u):
(w, theta) = split(u)
return theta-grad(w)

def bending_moment(u):
DD = as_tensor([[D, nuD, 0], [nuD, D, 0],[0, 0, D*(1-nu)/2.]])
return voigt2stress(dot(DD,strain2voigt(curv(u))))

def shear_force(u):
return F*shear_strain(u)

u = Function(V)
u_ = TestFunction(V)
du = TrialFunction(V)

(w_, theta_) = split(u_)
(dw, dtheta) = split(du)

dx_shear = dx(metadata={“quadrature_degree”: 2*deg-2})

#-----------------------------------------------------------------------------------------------------

W_ext1 = M_left * theta_[0]*ds(1)

W_ext2 = M_right * theta_[0]*ds(1)

#-----------------------------------------------------------------------------------------------------

l_form = Constant(1.)*u_[0]*dx + W_ext1 + W_ext2
k_form = inner(bending_moment(u_), curv(du))*dx + dot(shear_force(u_), shear_strain(du))*dx_shear

K = PETScMatrix()
b = PETScVector()

assemble_system(k_form, l_form, bc, A_tensor=K, b_tensor = b)

m_form = rho * thick * dot(dw, w_) * dx + (rho/12.) * (thick ** 3) * dot(dtheta, theta_) * dx + (rho / 12) * (thick ** 3) * dot(dtheta, theta_) dx
M = PETScMatrix()
assemble(m_form, tensor=M)

Hi,
sorry but your code is not properly formatted with ``` so its hard to read but it seems to me that you are not accounting for the rotational spring stiffness in the bilinear k_form. There you should add the rotational spring contribution using something like:

k_form += K_spring*theta_[0]*dtheta[0]*ds(1)

where K_spring is the spring stiffness.
Here you just added some Neumann terms in the linear form which is not even used in a modal analysis.

1 Like

Hi Jeremy, thanks for your response. I introduced the spring stiffness term in the bilinear form as suggested. However, the eigenfrequencies obtained by keeping the k_spring value very high and maintaining the simply supported case for translation dof restriction on the same edges (while other two edges free), I am not able to achieve an equivalent clamped case. What would you suggest??

Are you sure that ds(1) correctly corresponds to the two edges ? It does not seem to me as you seem to have defined twice ds with two subdomain_data, hence only the second is taken into account

1 Like

Hey Jeremy, thanks for your valuable support. The program worked and is pointing to a clamped case with very high value of the rotational spring stiffness. Thanks, once again!!!