Flat membrane elements in R3 for elastic problem

Question

Hello everybody,

I just started with FEniCS a couple of weeks ago and would like to use it as a framework for the structural simulation of a ram-air kite. Because the material is very thin, the analysis is based on membrane elements (only resisting in-plane membrane forces) embedded in R3 space.

I implemented a simple test case: a flat rectangle in R3, with body forces acting parallel to its surface (see code). I realised that the volume integration for the body forces is not working as I intended, and that the stress tensor has to be transformed to local coordinates of the element.

Also, according to the FEniCS book an embedded triangle element formulation in R3 has not been implemented. Is that still the case or has this issue been resolved? If not, is there an alternative approach?

from fenics import *
import numpy as np

editor = MeshEditor()
mesh = Mesh()

editor.open(mesh, 2, 3)  # top. and geom. dimension are 2 and 3
editor.init_vertices(4)  # number of vertices
editor.init_cells(2)     # number of cells

editor.add_vertex(0, np.array([0.0, 0.0, 0.0]))
editor.add_vertex(1, np.array([15.0, 0.0, 0.0]))
editor.add_vertex(2, np.array([15.0, 9.0, 16.0]))
editor.add_vertex(3, np.array([0.0, 9.0, 16.0]))

editor.add_cell(0, np.array([0, 1, 2], dtype=np.uintp))
editor.add_cell(1, np.array([0, 2, 3], dtype=np.uintp))
editor.close()

V = VectorFunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
d = u.geometric_dimension()
v = TestFunction(V)
info(V)

# Elasticity parameters
E, nu = 100, 0.3
mu, lmbda = Constant(E / (2 * (1 + nu))), Constant(E * nu / ((1 + nu) * (1 - 2 * nu)))

# Small strain
def epsilon(u):
    return 0.5 * (nabla_grad(u) + nabla_grad(u).T)

# Cauchy stress
def sigma(u):
    return lmbda * nabla_div(u) * Identity(d) + 2 * mu * epsilon(u)

class Boundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[1], 0, tol)

boundary = Boundary()
bc = DirichletBC(V, Constant((0.0, 0.0, 0.0)), boundary)

# Body forces
bVec = Constant((0.0, 3/5, 4/5))
RHS = dot(bVec, v) * dx

# Strain energy
LHS = inner(sigma(u), epsilon(v))*dx

u = Function(V)

solve(LHS == RHS, u, bc)

Best regards,

Paul