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

variational formulation
1

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

2

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
3 Likes
3

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