Problem with inflation of neo-hookean membrane

Hi all,
I have a problem simulating the inflation of a hyper-elastic membrane to increasing pressures.
Since I’m new in Fenics, I started from the demo code (https://fenicsproject.org/docs/dolfin/1.4.0/python/demo/documented/hyperelasticity/python/documentation.html) and adapt it as follows:

from dolfin import *
import numpy as np 
import pdb

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

# Create mesh and define function space
mesh = BoxMesh(Point(0.0, 0.0, 0.0),Point(1.0, 0.05, 1.0), 20, 10, 20)
#UnitCubeMesh(24, 16, 16)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

N = FacetNormal(mesh)

# Mark boundary subdomians
left =  CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)

top =  CompiledSubDomain("near(x[1], side) && on_boundary", side = 0.05)
bottom = CompiledSubDomain("near(x[1], side) && on_boundary", side = 0.0)

front =  CompiledSubDomain("near(x[2], side) && on_boundary", side = 1.0)
back = CompiledSubDomain("near(x[2], side) && on_boundary", side = 0.0)
# Define Dirichlet boundary (x = 0 or x = 1)
c = Expression(("0.0", "0.0", "0.0"),degree=3)


bcb = DirichletBC(V, c, bottom)
bct = DirichletBC(V, c, top)
bcl = DirichletBC(V, c, left)
bcr = DirichletBC(V, c, right)
bcf = DirichletBC(V, c, front)
bcba = DirichletBC(V, c, back)
#bcr = DirichletBC(V, r, right)
bcs = [bcl, bcr,bcf,bcba]

# Define functions
du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement from previous iteration
B  = Constant((0.0, 0.0, 0.0))  # Body force per unit volume
T  = Constant((0.0,  0.0, 0.0))  # Traction force on the boundary

# Kinematics
d = u.geometric_dimension()
I = Identity(d)             # Identity tensor

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

file = File("results/displacement.pvd");

dt=0.05
p_f=1
p_0=0.01
NN=(p_f-p_0)/dt
pressures=np.linspace(p_0,p_f,int(NN))
for p in pressures:
	print()
	print('pressure: ',p)
    	T  = Constant((0.0, p, 0.0))  # Body force per unit volume
    	F = I + grad(u)             # Deformation gradient
    	C = F.T*F                   # Right Cauchy-Green tensor
    	Finv=inv(F)
    	
    	# Invariants of deformation tensors
    	Ic = tr(C)
    	J  = det(F)

    	n = J*Finv.T*N

    	# Stored strain energy density (compressible neo-Hookean model)
    	psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2
    	pr=Constant((p))
    	# Total potential energy
    	Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds - dot(u,pr*n)*ds
            #Pi = psi*dx - dot(B, u)*dx - dot(u,pr*n)*ds
    	# Compute first variation of Pi (directional derivative about u in the direction of v)
    	F = derivative(Pi, u, v)

    	# Compute Jacobian of F
    	J = derivative(F, u, du)
    	# Solve variational problem
    	solve(F == 0, u, bcs, J=J,
          	form_compiler_parameters=ffc_options,
          	   solver_parameters={
    						"newton_solver":{\
    							"relaxation_parameter":0.5,
    							"relative_tolerance":1e-6 ,
    							"maximum_iterations":100, 
    						} })
    	#pdb.set_trace()


    	file << (u,p);

The idea of the for loop is that by increasing the pressure progressively, a partial solution can be computed and then used as an initial point for the next step.

I am not sure how to implement the internal pressure forces. I tried two approaches:

  • Changing the value to T = Constant((0.0, p, 0.0)), where p is the pressure value. This method results in applying forces to each cell at the boundary surface:

  • As pressure can be seen as a force that is applied normal to the surface, I used the following (following implementation of neuman bc from Fenics mechancis based on Nanson’s formula ):
    Pi = psi*dx - dot(B, u)dx - dot(u,prn)*ds. This results in a weirdly deformed membrane

I would appreciate any hints on the implementation of an inflatable membrane.
Sorry if the question/solution is obvious but my understanding of finite element is limited.
Thanks

The approach using Nanson formula would be the way to go. Your implementation seems OK to me except for the boundary measure ds. I would expect using ds(1) or something where 1 would be the portion of the boundary on which pressure applies. Here you exert pressure on the whole boundary which certainly is the cause of your “weirdly deformed membrane”

Thanks for the quick reply. That works. Here is the code in case someone is having the same issue:

class bottom_b(SubDomain):
    def inside(self,x,on_boundary):
        return on_boundary and near(x[1],0.0) 

boundaries = MeshFunction("size_t", mesh,mesh.topology().dim()-1, 0)
boundaries.set_all(0)
bottom_bound=bottom_b()
bottom_bound.mark(boundaries, 1)
ds_region = Measure('ds')[boundaries]

Pi = psi*dx - dot(B, u)*dx - dot(u,pr*n)*ds_region(1)