Newtonsolver does not converge - Hyperelasticity

1

Hello everyone,

I want to simulate a simple strain Test with the help of a hyperelastic Material.

For the simulation of the load I defined a boundary condition for the displacement which is increased for each time step. At first everything is fine, but after the 34 th timestep the newton solver doesn’t converge anymore and I get the error Message:

*** Error:   Unable to solve nonlinear system with NewtonSolver.
*** Reason:  Newton solver did not converge because maximum number of iterations reached.
*** Where:   This error was encountered inside NewtonSolver.cpp.
*** Process: 0

this has happened when I use a 2D mesh as well for a 3D one. Can you tell me what I did wrong and how I can change that?
Thanks a lot in advance!

Here is my code

from fenics import *
import time
import matplotlib.pyplot as plt

# Optimization options for the form compiler
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True

# Elasticity parameters
E = 220*1000
nu = 0.45
beta = Constant(nu/(1-2*nu))
mu= Constant(E/(2*(1 + nu)))

# Zeitvariablen
T = 10          # final time
num_steps = 100    # number of time steps
dt = T / num_steps # time step size


##############################################################################
#                                   Create mesh                               #
##############################################################################

# Create mesh and define function space 
nx0 = 0
ny0 = 0
nz0 = 0
nx1 = 0.11  # in m
ny1 = 0.015/2  # in m 
nz1 = 0.003/2 # in m
### number of elements in x.y und z-direction
nx = 20
ny = 6
nz = 3


mesh = BoxMesh(Point(nx0,ny0,nz0),Point(nx1,ny1,nz1), nx,ny, nz)


##############################################################################
#                                   FunctionSpace                      #
##############################################################################


P2 = VectorElement('Lagrange', mesh.ufl_cell(), 2)      
V = FunctionSpace(mesh, P2)



##############################################################################
#                                Boundary Stuff                              #
##############################################################################

############################### x-direction ####################################
class rechts(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[0],nx1,tol)
class links(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[0],nx0,tol)
# x    
right = rechts()
left = links()



boundaries = MeshFunction('size_t',mesh,mesh.topology().dim()-1)



strainrate = 1*10**(-3) # in m

u_RB_lin = Expression(("strainrate*t"), degree=2, strainrate=strainrate, t=0)
      
bcr = DirichletBC(V.sub(0),  u_RB_lin, right)
bcl = DirichletBC(V, Constant((0, 0, 0)), left)

bcs = [bcl, bcr]


##############################################################################
#                                Functions                                   #
##############################################################################

# Define functions
du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement from previous iteration

# Kinematics
d = u.geometric_dimension()
I = Identity(d)             # Identity tensor
F = I + grad(u)             # Deformation gradient
invF = inv(F)
C = F.T*F                   # Right Cauchy-Green tensor

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

# stress tensor 
Piola = mu*(-J**(-2*beta)*invF.T)+mu*F

# Momentum balance (no external forces)
Imp = inner(Piola,nabla_grad(v))*dx

eqn = Imp

##############################################################################
#                                Solving                                   #
##############################################################################


t = 0   # Definieren des Zeitpunkts t0
for n in range(num_steps):
    start = time.time()

    # Update current time
    t += dt
    u_RB_lin.t=t
 
    # Compute Jacobian of F
    J = derivative(eqn, u, du)
    
    # Solve variational problem
    solve(eqn == 0, u, bcs, J=J,solver_parameters={"newton_solver":
                                            {"relative_tolerance": 1e-6}})
2

Hi, the error is quite explicit. Try increasing the number of time steps to see if it helps.

3

Hello,
thank you for your answer!

I changed the number of time steps (num_steps) from 100 to 5000 and for the 2D Modell it helped, but for the 3D one it still didn’t converge after some time-steps.
Is there a solver-option or something similar, that helps with that kind of problems?

4

Dear Martin

I am not sure if your formulation is the main issue, try to use another material model, like neo-Hooke. Technically, length scale in meter means that the stress is in Pa and Young’s modulus in your code is 220 kPa. So it is a soft matter, I would expect that Poisson’s ratio is near to 0.5 (isochoric material). In this case the parameter beta goes to infinity and your system will be impossible to solve.

Best,
Emek
http://bilenemek.abali.org

5

Hello Emek,

I have used a Neo Hook Material. I used the free Energy in the form that was presented by Holzapfel:
Nonlinear Solid Mechanics: A Continuum Approach for Enineering: A Continuum Approach for Engineering – 23. März 2000 p.247

Psi = c_1/beta (J^{-2beta}-1) + c_1(tr(C ) -3 )
c_1 = mu/2
beta = nu/(1-2nu))
C = F*F^T
mu - 2. Lamee Constant
nu - Poissonratio
F - Deformationgradient

and derivated the 1st Piola from this energyfunction and then implemented it into fenics.

Or could the error be, that I have to assign my displacement from the previous time step to the new time step? Or is that not necessary?

6

It seems the problem was the number of elements of the mesh.
I increased the number of elements and then it works. :slight_smile: