Error while solving hyper-elastic problem using Neo-Hookean model based on principal stretches

Hello,

I am trying to solve a hyperelastic problem using the Neo-Hookean model based on principal stretches. The standard (V1) neo Hookean model is given by

\psi = \frac{\mu}{2} (I_{c} - 2 - 2\mu \ln(J) )+ \frac{\lambda}{2}\ln(J)^{2}

where

I_{c} = {\rm trace}(C)

A secod version (V2) of the model, which is based on the eigen values (\lambda_1, \lambda_2) of right Cauchy-Green tensor (C) is given by

\psi = \frac{\mu}{2} (\Lambda_{c} - 2 - 2\mu \ln(J_m) )+ \frac{\lambda}{2}\ln(J)^{2}

where

\Lambda_{c} = \lambda_1 + \lambda_2
J_{m} = \sqrt{\lambda_1 * \lambda_2}

I have attempted to formulate this formulation and the V1 Neo-Hookean model works fine. But, when I try to run the V2 model the code gives me the following error.

Solving nonlinear variational problem.
    Newton iteration 0: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
    Newton iteration 1: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
    Newton iteration 2: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
    Newton iteration 3: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
    Newton iteration 4: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
*** -------------------------------------------------------------------------
*** 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
*** 
*** DOLFIN version: 2019.1.0
*** Git changeset:  74d7efe1e84d65e9433fd96c50f1d278fa3e3f3f
*** -------------------------------------------------------------------------

Here is my minimum working code. You can run both models by changing the variable “model”.

from dolfin import *
from ufl import grad as ufl_grad

mesh = UnitSquareMesh(32, 32)
# --------------------------------------------------------------------------
U = VectorFunctionSpace(mesh, "CG", 1)
u, v = TrialFunction(U), TestFunction(U)
u_new = Function(U)
# --------------------------------------------------------------------------
support = CompiledSubDomain("on_boundary && near(x[0], 0, tol)", tol=1e-14)
force_face = CompiledSubDomain("on_boundary && near(x[0], 1, tol)", tol=1e-14)
applied_disp = Expression("disp", disp=0.0, degree=2)
# --------------------------------------------------------------------------
bc_sp = DirichletBC(U, Constant((0.0,0.0)), support)
bc_force_x = DirichletBC(U.sub(0), applied_disp, force_face)
bcs = [bc_sp, bc_force_x]
# --------------------------------------------------------------------------
dim = 2
mu = 1
lmbda = 1
# --------------------------------------------------------------------------
def Grad(v):
    return ufl_grad(v)

def InfinitesimalStrain(u):
    return variable(0.5 * (Grad(u) + Grad(u).T))

def SecondOrderIdentity(u):
    return variable(Identity(u.geometric_dimension()))

def DeformationGradient(u):
    I = SecondOrderIdentity(u)
    return variable(I + Grad(u))

def Jacobian(u):
    F = DeformationGradient(u)
    return variable(det(F))

def RightCauchyGreen(u):
    F = DeformationGradient(u)
    return variable(F.T * F)
# --------------------------------------------------------------------------

model = "nstandard"
def neo_hookean(u):

    J = Jacobian(u)
    C = RightCauchyGreen(u)

    eig1 = 0.5 * (tr(C) + sqrt(tr(C) ** 2 - 4 * det(C)))
    eig2 = 0.5 * (tr(C) - sqrt(tr(C) ** 2 - 4 * det(C)))
    
    if model == "standard":
        # Standard neo hookean parameters
        Ic = tr(C)
        Jm = J
    else:
        # Neo hookean parameters based on principal streches
        Ic = eig1 + eig2
        Jm = sqrt(eig1 * eig2)
    
    si = (mu / 2) * (Ic - dim - 2 * ln(Jm)) + (lmbda / 2) * (ln(J)) ** 2
    return si 
# --------------------------------------------------------------------------
total_energy = neo_hookean(u_new) * dx 
weak_form = derivative(total_energy, u_new, v)
jacobian = derivative(weak_form, u_new, u)
# --------------------------------------------------------------------------
problem_disp = NonlinearVariationalProblem(weak_form, u_new, bcs, jacobian)
solver_disp = NonlinearVariationalSolver(problem_disp)
# --------------------------------------------------------------------------
applied_disp.disp = 0.001
(iter_num, converged) = solver_disp.solve()
print(iter_num)

I would really appreciate any help to make the second version work.
Thanks

You need to be careful when using UFL’s symbolic differentiation of sqrt when the argument goes to zero. The following modified code using a regularized square root function converges:

def safeSqrt(x):
    return sqrt(x+DOLFIN_EPS)

model = "nstandard"
def neo_hookean(u):

    J = Jacobian(u)
    C = RightCauchyGreen(u)

    eig1 = 0.5 * (tr(C) + safeSqrt(tr(C) ** 2 - 4 * det(C)))
    eig2 = 0.5 * (tr(C) - safeSqrt(tr(C) ** 2 - 4 * det(C)))
    
    if model == "standard":
        # Standard neo hookean parameters
        Ic = tr(C)
        Jm = J
    else:
        # Neo hookean parameters based on principal streches
        Ic = eig1 + eig2
        Jm = safeSqrt(eig1 * eig2)
    
    si = (mu / 2) * (Ic - dim - 2 * ln(Jm)) + (lmbda / 2) * (ln(J)) ** 2
    return si 
1 Like

Thank you so much. I have been struggling with it for the past two days. This solved the issue.