Axisymmetric formulation for hyperelastic structures of revolution

variational formulation
1

Hi,

I have been trying to implement a hyperelastic material model defined here using an axisymmetric formulation for a cylinder (a 2 dimensional slice of). I have gone through what I can to account for missing terms in the 2D representation, changing

Pi = psi*df.dx

to

Pi = psi*2*np.pi*X[1]*df.dx

for example (revolving around X[0]).

FEniCS returns a reasonable looking solution for the displacement field, but I’m struggling to work out if I need to change the lines

residual = df.derivative(Pi, u, v)
Jacobian = df.derivative(residual, u, du)

to account for missing terms in the 2D representation. I’m not sure how I could change them very easily if I do need to. I have left my reduced code below and I’m using FEniCS version 2019.2.0.dev0

Thanks,
Kit

import dolfin as df, numpy as np

# Define mesh and displacement of RHS
L, R, density = 10, 3, 10
mesh = df.RectangleMesh(df.Point(0, 0), df.Point(L, R), L*density, R*density)
displacement = 2

# Define function space
V = df.VectorFunctionSpace(mesh, 'CG', 1)
du = df.TrialFunction(V)
v = df.TestFunction(V)
u = df.Function(V)

# Define material parameters
k0 = 1e3
mu0 = 50
k1 = 1e3
k2 = 2
a04 = df.Constant((1.0, 0))

# Kinematics
X = df.SpatialCoordinate(mesh) # X[0] = Z, X[1] = R
dimension = mesh.geometry().dim()
I         = df.Identity(dimension)           # Identity tensor
F         = I + df.grad(u)                   # Deformation gradient
C         = F.T * F                          # Right Cauchy-Green tensor
B         = F * F.T                          # Left Cauchy-Green tensor
J_naive   = df.det(F)                        # Naive determinant of F
J         = (1 + u[1]/X[1]) * J_naive        # Corrected determinant of F
I1_naive  = df.tr(C)                         # Naive first invariant of C
I1        = I1_naive + (1 + u[1]/X[1]) ** 2  # Corrected first invariant of C
I1_ico    = J ** (-2/3) * I1                 # Isochoric first invariant of C
I4        = df.dot(a04, C * a04)             # Anisotropic invariant

# Strain energy density
psi_vol = 1/2 * k0 * (J - 1) ** 2
psi_ico = 1/2 * mu0 * (I1_ico - 3)
psi_ani = (k1 / (2 * k2)) * (df.exp(k2 * (I4 - 1) ** 2) - 1)
psi = psi_vol + psi_ico + psi_ani

# Total potential energy
Pi = psi*2*np.pi*X[1]*df.dx

# Define boundaries: x[0] = Z, x[1] = R
x = np.unique(mesh.coordinates()[:,0])
y = np.unique(mesh.coordinates()[:,1])
class Left(df.SubDomain):
    def inside(self, x, on_boundary):
        return df.near(x[0], min(x)) and on_boundary
class Right(df.SubDomain):
    def inside(self, x, on_boundary):
        return df.near(x[0], max(x)) and on_boundary
class Axis(df.SubDomain):
    def inside(self, x, on_boundary):
        return df.near(x[1], min(y)) and on_boundary
class CurvedSurface(df.SubDomain):
    def inside(self, x, on_boundary):
        return df.near(x[1], max(y)) and on_boundary
    
# Mark facets
facets = df.MeshFunction("size_t", mesh, 1)
facets.set_all(0)
Left().mark(facets, 1)
Right().mark(facets, 2)
Axis().mark(facets, 3)
CurvedSurface().mark(facets, 4)

# Clamped at left, sliding at right, sliding on axis
bcs = [df.DirichletBC(V, df.Constant((0., 0.)), facets, 1),
        df.DirichletBC(V.sub(0), df.Constant(displacement), facets, 2),
        df.DirichletBC(V.sub(1), df.Constant(0.), facets, 3)]

# Displacement field solution
residual = df.derivative(Pi, u, v)
Jacobian = df.derivative(residual, u, du)
df.solve(residual == 0, u, bcs, J = Jacobian, 
solver_parameters={"newton_solver":{"linear_solver":"lu"}})

df.plot(u, mode='displacement')