Hi,
I am working on the problem of finding the natural frequencies and mode shapes of a semi-circle which is clamped around the circumferential edge and free along the straight radial edge. I have managed to implement the fourth order bending equation $$ \nabla^{4} w = \lambda w $$ for a circle with clamped edge boundary conditions using a mixed method.
However when I change the geometry to a semi circle I need to apply the zero moment following boundary conditions at the free straight radial edge:
$$ M_{\theta} = \nu \dfrac{\partial^2 w}{\partial r^2} + \dfrac{1}{r}\dfrac{\partial w}{\partial r} + \dfrac{1}{r^2}\dfrac{\partial^2 w}{\partial \theta^2} = 0 $$
$$ V_{\theta} = \dfrac{1}{r}\dfrac{\partial}{\partial \theta} \left( \dfrac{\partial^2 w}{\partial r^2} + \dfrac{1}{r}\dfrac{\partial w}{\partial r} + \dfrac{1}{r^2}\dfrac{\partial^2 w}{\partial \theta^2} \right) + \left( 1 - \nu \right) \dfrac{\partial}{\partial r} \left[ \dfrac{\partial}{\partial r} \left( \dfrac{1}{r} \dfrac{\partial w}{\partial \theta} \right) \right] $$
Where \nu is the poisson's ratio of the material. Is there any way to fix certain differential equations to be satisfied on the boundary in FEnics. I have tried putting these equations into a format where they are dirichlet and neumman boundary conditions on $$ \nabla^2 w $$ and come up with the following
$$ \nabla^2 w = \left( 1 - \nu \right) \dfrac{\partial^2 w}{\partial r^2} $$
and
$$ \nabla \left( \nabla^2 w \right) \cdot n = \left( \nu - 1 \right)\dfrac{\partial^2}{\partial r^2} \left( \dfrac{1}{r} \dfrac{\partial w}{\partial \theta} \right) $$
But this didn't help much as I don't know if there is a way of implementing boundary conditions which depend on the solution and also what the syntax would be for \dfrac{\partial^2 w}{\partial r^2}.
Here is the python code I am using without any boundary condition applied to the free edge.
from dolfin import *
from mshr import *
import numpy
import math
import matplotlib.pyplot as plt
import csv
# Construct Mesh
radius = 1
r1 = Circle(Point(0.0, 0.0), radius)
r2 = Rectangle(dolfin.Point(-2.0*radius, -radius), dolfin.Point(0.0, radius))
# Define domain and resolution
domain = r1-r2
res = 10
# Generate mesh
mesh = generate_mesh(domain, res)
# Define function spaces and mixed (product) space
Space1 = FunctionSpace(mesh, "Lagrange", 1)
Space2 = FunctionSpace(mesh, "Lagrange", 1)
W = Space1 * Space2
def fixed_boundary(x,on_boundary):
r = math.sqrt(x[0] * x[0] + x[1] * x[1])
return on_boundary and (near(r,radius,0.01))
def flat_boundary_free(x,on_boundary):
return on_boundary and near(x[0],0,0.01) and abs(x[1])<1
bc2 = DirichletBC(W.sub(1),Constant(0.0),fixed_boundary)
# Define trial and test functions
(sigma, u) = TrialFunctions(W)
(tau, v) = TestFunctions(W)
# Define variational form
a = (inner(grad(u), grad(tau))+inner(grad(sigma),grad(v))+(sigma*tau))*dx
m = -u*v*dx
# Assemble stiffness form
A = PETScMatrix()
assemble(a, tensor=A)
M = PETScMatrix()
assemble(m, tensor=M)
bc2.apply(A)
bc2.apply(M)
# Create eigensolver
eigensolver = SLEPcEigenSolver(A,M)
eigensolver.parameters['spectrum'] = 'smallest magnitude'
eigensolver.parameters['solver'] = 'lapack'
# Compute all eigenvalues of A x = \lambda x
waitMessage = "Computing eigenvalues. Res ="+str(res)
print(waitMessage)
eigensolver.solve()
n = eigensolver.get_number_converged()
x = []
i = 0
#extract next eigenpair
print("Looping to find first 9 non-spurious eigenvalues")
while True:
r, c, rx, cx = eigensolver.get_eigenpair(i)
if abs(r-1) > 0.03 and abs(r)>0.01:
if len(x) < 10:
x.append(i)
print(i)
print 'Sqrt eigenvalue:', math.sqrt(r)
if len(x) > 8:
break
i = i+1
for j in range(0,9):
r, c, rx, cx = eigensolver.get_eigenpair(x[j])
w = Function(W)
w.vector()[:] = rx
plot(w[1])
interactive()
This code produces similar results to analytical solutions when a fully circular plate is used.
This seems to produce sensible plots for the semi circle but I need to ensure the free edge boundary condition is satisfied and I am not sure how to implement this in FEnics. Any help would be greatly appreciated. Thanks
