Partly Clamped Hyperbolic Paraboloid

This demo is implemented in the single Python file demo_naghdi-linear-hypar.py.

This demo program solves the linear Naghdi shell equations for a partly-clamped hyperbolic paraboloid (hypar) subjected to a uniform vertical loading. This is a well known bending dominated benchmark problem for testing a FE formulation with respect to membrane locking issues, see [1]. Here, locking is cured using enriched finite elements including cubic bubble shape functions and Partial Selective Reduced Integration (PSRI), see [2, 3].

../../_images/configuration.png

To follow this demo you should know how to:

  • Define a MixedElement and EnrichedElement and a FunctionSpace from it.
  • Write variational forms using the Unified Form Language.
  • Automatically derive Jacobian and residuals using derivative().
  • Apply Dirichlet boundary conditions using DirichletBC and apply().

This demo then illustrates how to:

  • Define and solve a linear Naghdi shell problem with a curved stress-free configuration given as analytical expression in terms of two curvilinear coordinates.

We start with importing the required modules, setting matplotlib as plotting backend, and generically set the integration order to 4 to avoid the automatic setting of FEniCS which would lead to unreasonably high integration orders for complex forms.:

import os, sys

import numpy as np
import matplotlib.pyplot as plt

from dolfin import *
from mshr import *
from ufl import Index

from mpl_toolkits.mplot3d import Axes3D

parameters["form_compiler"]["quadrature_degree"] = 4
output_dir = "output/"
if not os.path.exists(output_dir):
    os.makedirs(output_dir)

We consider an hyperbolic paraboloid shell made of a linear elastic isotropic homogeneous material with Young modulus Y and Poisson ratio nu; mu and lb denote the shear modulus \(\mu = Y/(2 (1 + \nu))\) and the Lamé constant \(\lambda = 2 \mu \nu/(1 - 2 \nu)\). The (uniform) shell thickness is denoted by t.

Y, nu = 2.0e8, 0.3
mu = Y/(2.0*(1.0 + nu))
lb = 2.0*mu*nu/(1.0 - 2.0*nu)
t = Constant(1E-2)

The midplane of the initial (stress-free) configuration \(\mathcal{S}_I\) of the shell is given in the form of an analytical expression

\[y_I:x\in\mathcal{M}\subset R^2\to y_I(x)\in\mathcal{S}_I\subset \mathcal R^3 \qquad \mathcal{S}_I = \{x_1, x_2, x_1^2-x_2^2\}\]

in terms of the curvilinear coordinates \(x\). Hence, we mesh the two-dimensional domain \(\mathcal{M}\equiv [-L/2,L/2]\times [-L/2,L/2]\).:

L = 1.0
P1, P2 = Point(-L/2, -L/2), Point(L/2, L/2)
ndiv = 40
mesh = RectangleMesh(P1, P2, ndiv, ndiv)

We provide the analytical expression of the initial shape as an Expression that we represent on a suitable FunctionSpace (here \(P_2\), but other are choices are possible):

initial_shape = Expression(('x[0]','x[1]','x[0]*x[0] - x[1]*x[1]'), degree = 4)
V_y = FunctionSpace(mesh, VectorElement("P", triangle, degree=2, dim=3))
yI = project(initial_shape, V_y)

We compute the covariant and contravariant component of the metric tensor, aI, aI_contra; the covariant base vectors g0, g1; the contravariant base vectors g0_c, g1_c.

aI = grad(yI).T*grad(yI)
aI_contra, jI = inv(aI), det(aI)
g0, g1 = yI.dx(0), yI.dx(1)
g0_c, g1_c = aI_contra[0,0]*g0 + aI_contra[0,1]*g1, aI_contra[1,0]*g0 + aI_contra[1,1]*g1

Given the midplane, we define the corresponding unit normal as below and project on a suitable function space (here \(P_1\) but other choices are possible):

def normal(y):
    n = cross(y.dx(0), y.dx(1))
    return n/sqrt(inner(n,n))

V_normal = FunctionSpace(mesh, VectorElement("P", triangle, degree = 1, dim = 3))
nI = project(normal(yI), V_normal)

We can visualize the shell shape and its normal with this utility function

def plot_shell(y,n=None):
    y_0, y_1, y_2 = y.split(deepcopy=True)
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    ax.plot_trisurf(y_0.compute_vertex_values(),
                    y_1.compute_vertex_values(),
                    y_2.compute_vertex_values(),
                    triangles=y.function_space().mesh().cells(),
                    linewidth=1, antialiased=True, shade = False)
    if n:
        n_0, n_1, n_2 = n.split(deepcopy=True)
        ax.quiver(y_0.compute_vertex_values(),
              y_1.compute_vertex_values(),
              y_2.compute_vertex_values(),
              n_0.compute_vertex_values(),
              n_1.compute_vertex_values(),
              n_2.compute_vertex_values(),
              length = .2, color = "r")
    ax.view_init(elev=50, azim=30)
    ax.set_xlim(-0.5, 0.5)
    ax.set_ylim(-0.5, 0.5)
    ax.set_zlim(-0.5, 0.5)
    ax.set_xlabel(r"$x_0$")
    ax.set_ylabel(r"$x_1$")
    ax.set_zlabel(r"$x_2$")
    return ax

plot_shell(yI)
plt.savefig("output/initial_configuration.png")
../../_images/initial_configuration.png

In our 5-parameter Naghdi shell model the configuration of the shell is assigned by

  • the 3-component vector field u_ representing the (small) displacement with respect to the initial configuration yI
  • the 2-component vector field theta_ representing the (small) rotation of fibers orthogonal to the middle surface.

Following [2, 3], we use a P2+bubble element for y_ and a P2 element for theta_, and collect them in the state vector z_=[u_, theta_]. We further define Function, TestFunction, and TrialFucntion and their different split views, which are useful for expressing the variational formulation.

P2 = FiniteElement("P", triangle, degree = 2)
bubble = FiniteElement("B", triangle, degree = 3)

Z = FunctionSpace(mesh, MixedElement(3*[P2 + bubble ] + 2*[P2]))
z_ = Function(Z)
z, zt = TrialFunction(Z), TestFunction(Z)

u0_, u1_, u2_, th0_, th1_ = split(z_)
u0t, u1t, u2t, th0t, th1t = split(zt)
u0, u1, u2, th0, th1 = split(z)

We define the displacement vector and the rotation vector, with this latter tangent to the middle surface, \(\theta = \theta_\sigma \, g^\sigma\), \(\sigma = 0, 1\)

u_, u, ut = as_vector([u0_, u1_, u2_]), as_vector([u0, u1, u2]), as_vector([u0t, u1t, u2t])
theta_, theta, thetat = th0_*g0_c + th1_*g1_c, th0*g0_c + th1*g1_c, th0t*g0_c + th1t*g1_c

The extensional, e_naghdi, bending, k_naghdi, and shearing, g_naghdi, strains in the linear Naghdi model are defined by

e_naghdi = lambda v: 0.5*(grad(yI).T*grad(v) + grad(v).T*grad(yI))
k_naghdi = lambda v, t: -0.5*(grad(yI).T*grad(t) + grad(t).T*grad(yI)) - 0.5*(grad(nI).T*grad(v) + grad(v).T*grad(nI))
g_naghdi = lambda v, t: grad(yI).T*t + grad(v).T*nI

Using curvilinear coordinates, the constitutive equations are written in terms of the matrix A_hooke below, representing the contravariant components of the constitutive tensor for isotropic elasticity in plane stress, see e.g. [4]. We use the index notation offered by UFL to express operations between tensors

i, j, k, l = Index(), Index(), Index(), Index()
A_hooke = as_tensor((((2.0*lb*mu)/(lb + 2.0*mu))*aI_contra[i,j]*aI_contra[k,l] + 1.0*mu*(aI_contra[i,k]*aI_contra[j,l] + aI_contra[i,l]*aI_contra[j,k])),[i,j,k,l])

The membrane stress and bending moment tensors, N and M, and shear stress vector, T, are

N = as_tensor((A_hooke[i,j,k,l]*e_naghdi(u_)[k,l]),[i, j])
M = as_tensor(((1./12.0)*A_hooke[i,j,k,l]*k_naghdi(u_, theta_)[k,l]),[i, j])
T = as_tensor((mu*aI_contra[i,j]*g_naghdi(u_, theta_)[j]), [i])

Hence, the contributions to the elastic energy densities due to membrane, psi_m, bending, psi_b, and shear, psi_s, are

psi_m = .5*inner(N, e_naghdi(u_))
psi_b = .5*inner(M, k_naghdi(u_, theta_))
psi_s = .5*inner(T, g_naghdi(u_, theta_))

Shear and membrane locking are treated using the PSRI proposed by Arnold and Brezzi, see [2, 3]. In this approach, shear and membrane energies are splitted as a sum of two weighted contributions, one of which is computed with a reduced integration. Thus, shear and membrane energies have the form

\[(i=m,s) \qquad \alpha \int_\mathcal{M} \psi_i \, \sqrt{j_I} \mathrm{d}x + (\kappa - \alpha) \int_\mathcal{M} \psi_i \, \sqrt{j_I} \mathrm{d}x_h, \qquad \text{where} \, \kappa \propto t^{-2}\]

While [2, 3] suggest a 1-point reduced integration, we observed that this leads to spurious modes in the present case. We use then \(2\times 2\)-points Gauss integration for the portion \(\kappa - \alpha\) of the energy, whilst the rest is integrated with a \(4\times 4\) scheme. As suggested in [3], we adopt an optimized weighting factor \(\alpha = 1\)

dx_h = dx(metadata={'quadrature_degree': 2})
alpha = 1.0
kappa = 1.0/t**2
shear_energy = alpha*psi_s*sqrt(jI)*dx + (kappa - alpha)*psi_s*sqrt(jI)*dx_h
membrane_energy = alpha*psi_m*sqrt(jI)*dx + (kappa - alpha)*psi_m*sqrt(jI)*dx_h
bending_energy =  psi_b*sqrt(jI)*dx

Then, the elastic energy is

elastic_energy = (t**3)*(bending_energy + membrane_energy + shear_energy)

The shell is subjected to a constant vertical load. Thus, the external work is

body_force = 8.*t
f = Constant(body_force)
external_work = f*u2_*sqrt(jI)*dx

We now compute the total potential energy with its first and second derivatives

Pi_total = elastic_energy - external_work
residual = derivative(Pi_total, z_, zt)
hessian = derivative(residual, z_, z)

The boundary conditions prescribe a full clamping on the \(x_0 = -L/2\) boundary, while the other sides are left free

left_boundary = lambda x, on_boundary: abs(x[0] + L/2) <= DOLFIN_EPS and on_boundary
clamp = DirichletBC(Z, project(Expression(("0.0", "0.0", "0.0", "0.0", "0.0"), degree = 1), Z), left_boundary)
bcs = [clamp]

We now solve the linear system of equations

output_dir = "output/"
A, b = assemble_system(hessian, residual, bcs=bcs)
solver = PETScLUSolver("mumps")
solver.solve(A, z_.vector(), b)
u0_h, u1_h, u2_h, th0_h, th1_h = z_.split(deepcopy=True)

Finally, we can plot the final configuration of the shell

scale_factor = 1e4
plot_shell(project(scale_factor*u_ + yI, V_y))
plt.savefig("output/finalconfiguration.png")

References

[1] K. J. Bathe, A. Iosilevich, and D. Chapelle. An evaluation of the MITC shell elements. Computers & Structures. 2000;75(1):1-30.

[2] D. Arnold and F.Brezzi, Mathematics of Computation, 66(217): 1-14, 1997. https://www.ima.umn.edu/~arnold//papers/shellelt.pdf

[3] D. Arnold and F.Brezzi, The partial selective reduced integration method and applications to shell problems. Computers & Structures. 64.1-4 (1997): 879-880.

[4] P. G. Ciarlet. An introduction to differential geometry with applications to elasticity. Journal of Elasticity, 78-79(1-3):1–215, 2005.