Error occuring sometimes when solving the hyperelastic problem

Linear Algebra
1

Hello everyone,

I’m trying to solve the hyper-elastic problem with a neo-Hookean stored energy model, for a beam of size 10x1x1 under gravity. I use the classical example of elasticity for fenics (found here : https://fenicsproject.org/docs/dolfin/1.4.0/python/demo/documented/hyperelasticity/python/documentation.html) and modify it just a little bit to suit better my problem. I’ll post my code at the end of this message.

My issue is that the solver does not always converges, but depending on the mesh. For example if I refine my beam (initially formed with 8 points, at each extremity) 3 times and add a volume inside (I get 440 points then), the solver diverges when I use rho = 1500. But if I use rho = 1550 or 1450, it converges. Also if I refine it 4 times, it converges for rho = 1500.

Do you know where the issue might come from ?

My code :

V = VectorFunctionSpace(mesh, "Lagrange", 1)

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

# Define Dirichlet boundary (x = 0 or x = 1)
c = Constant((0, 0, 0))
bcl = DirichletBC(V, c, left)

# Material & Elasticity parameters

g = 9.81

E = 1e7
nu = 0.3
rho = 1500.
mu = Constant(E/(2*(1 + nu)))
lmbda = Constant(E*nu/((1 + nu)*(1 - 2*nu)))

# 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, -rho * g))  # 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
F = I + grad(u)             # Deformation gradient
C = F.T*F                   # Right Cauchy-Green tensor

# Invariants of deformation tensors
Ic = tr(C)
J  = det(F)

# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2

# Total potential energy
Pi = psi*dx - dot(B, u)*dx - dot(T, u)*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, bcl, J=J,
      form_compiler_parameters=ffc_options)
2

Hey Clement, did you solve your problem? I have a possibly alike problem.
My simulation diverges, when displacement rises over meshing width.