Navier plate in Python with FEniCS

navier_plate_fenics.py

# coding: utf-8
"""Simply supported rectangular plate on its four edges with FEniCS

Author: Juan Luis Cano Rodríguez <[email protected]>

References
----------
* Timoshenko, Stephen, and S. Woinowsky-Krieger. "Theory of Plates and Shells".
  New York: McGraw-Hill, 1959.

"""
import sys
from dolfin import *

num_elem = 32
if len(sys.argv) > 1:
    num_elem = int(sys.argv[1])

## Mixed approach

# --- Problem data
a = 1.0  # m
b = 1.0  # m
h = 0.050  # m

E = 69e9  # Pa
nu = 0.35

P = -10e3  # N
# ---

# Square mesh
mesh = UnitSquareMesh(num_elem, num_elem)

# Lagrange, 2nd order
V = FunctionSpace(mesh, "CG", 3)

u0 = Constant("0")

def u0_boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u0, u0_boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

D = E * h**3 / (12 * (1 - nu**2))

point = Point(0.5, 0.5)
f = PointSource(V, point, P)

# Weak form of the equation
a = Constant(D) * inner(nabla_grad(u), nabla_grad(v)) * dx
L = Constant(0.0) * v * dx

A, b = assemble_system(a, L, bc)
f.apply(b)

u = Function(V)
U = u.vector()
solve(A, U, b)

w = TrialFunction(V)
a = inner(nabla_grad(w), nabla_grad(v)) * dx
L = u*v*dx

w = Function(V)
solve(a == L, w, bc)

# Plot solution
plot(w, interactive=True)

# Print maximum displacement
max_w = max(abs(w.vector().array()))
print "Maximum displacement = %14.12f mm" % (max_w * 1e3)