Clamped circular plate bending solution error

1

I am very new to FEniCS and I am currently trying to solve the deflection of a clamped circular plate of radius a under uniform load q . The governing equation is ∇^4u=f , in which f=q/D is constant in my case. I have just adapted the biharmonic demo code from FEniCS website, changing only the mesh type and the forcing function f definition. The code runs perfectly, but the results do not match the analytical solution w=q(a^2−(x^2+y^2))/64D. What am I doing wrong? Also, sorry for the bad formatting of this question, its my first post here.


from dolfin import *

from mshr import *

import math

# Optimization options for the form compiler

parameters["form_compiler"]["cpp_optimize"] = True

parameters["form_compiler"]["optimize"] = True

# A mesh is created, and a quadratic finite element function space::

# Make mesh ghosted for evaluation of DG terms

parameters["ghost_mode"] = "shared_facet"

# Create mesh and define function space

raio = 10e-3  # Radius

q = 9e6       # Uniform pressure

E = 6.9e10    # Young Modulus

h = 2e-3      # Thickness

v = 0.33      # Poisson ratio

D = E*math.pow(h,3.0) / (12*(1-math.pow(v,2.0)))

domain = Circle(Point(0, 0), raio)

mesh = generate_mesh(domain, 64)

V = FunctionSpace(mesh, "CG", 2)

# Define Dirichlet boundary

class DirichletBoundary(SubDomain):

    def inside(self, x, on_boundary):

        return on_boundary

# The Dirichlet boundary condition is created::

# Define boundary condition

u0 = Constant(0.0)

bc = DirichletBC(V, u0, DirichletBoundary())

# Define trial and test functions

u = TrialFunction(V)

v = TestFunction(V)

# Define normal component, mesh size and right-hand side

h = CellDiameter(mesh)

h_avg = (h('+') + h('-'))/2.0

n = FacetNormal(mesh)

f = Constant(q/D)

analitico = Expression("q*pow(pow(a,2)-pow(x[0],2)-pow(x[1],2),2)/(64*D)",q=q,a=raio,D=D,degree=3)

# Penalty parameter

alpha = Constant(8.0)

# The bilinear and linear forms are defined::

# Define bilinear form

a = inner(div(grad(u)), div(grad(v)))*dx \

  - inner(avg(div(grad(u))), jump(grad(v), n))*dS \

  - inner(jump(grad(u), n), avg(div(grad(v))))*dS \

  + alpha/h_avg*inner(jump(grad(u),n), jump(grad(v),n))*dS

# Define linear form

L = f*v*dx

# Solve variational problem

u = Function(V)

solve(a == L, u, bc)

# Error

error_L2 = errornorm(analitico, u, 'L2')

print('error_L2 =', error_L2)
2

Please properly format your code using 3x` Encapsulation

1 Like
3

Running your code, I get the output

error_L2 = 5.465542724341058e-07

which is a very good match, considering that you are approximating a fourth order polynomial with a set of second order polynomials.

I am not sure what error you expect to get?