Variational formulation for hyperelastic material based on spectral decomposition

variational formulation
1

Hi, all. I’m trying to solve a problem of a neo-hookean hyperelastic model based on the spectral decompostion, but 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
Can you tell me where I’m wrong and how I can change it? Thanks a lot in advance!
The code is attached below
from future import division
from dolfin import *
from mshr import *
from scipy import optimize

# ----------------------------------------------------------------------------
# Parameters for DOLFIN and SOLVER
# ----------------------------------------------------------------------------
set_log_level(LogLevel.WARNING)  # 20, // information of general interest
# set some dolfin specific parameters
# ----------------------------------------------------------------------------
parameters["form_compiler"]["representation"]="uflacs"
parameters["form_compiler"]["optimize"]=True
parameters["form_compiler"]["cpp_optimize"]=True
parameters["form_compiler"]["quadrature_degree"]=4
info(parameters,True)

ut      = 1.0 # reference value for the loading (imposed displacement)

# Numerical parameters of the alternate minimization
maxiteration = 2000
AM_tolerance = 1e-4

# Mesh generation
mesh = UnitSquareMesh.create(1,1, CellType.Type.quadrilateral)

#define function space
W=VectorFunctionSpace(mesh,'CG',1)#displacement

u,v=TrialFunction(W),TestFunction(W)
u_new=Function(W)
u_old=Function(W)


E=10
nu=0.3
mu=Constant(E/(2*(1 + nu)))
K0=Constant(E/(3*(1-2*nu)))

#define boundary condition
support = CompiledSubDomain("on_boundary && near(x[1], 0, tol)", tol=1e-14)
force_face = CompiledSubDomain("on_boundary && near(x[1], 1, tol)", tol=1e-14)
applied_disp = Expression("disp", disp=0.0, degree=2)
# --------------------------------------------------------------------------
bc_bot = DirichletBC(W, Constant((0.0, 0.0)), support)
bc_top = DirichletBC(W.sub(1), applied_disp, force_face)
bcu = [bc_bot, bc_top]

boundaries=MeshFunction("size_t",mesh,mesh.topology().dim()-1)
boundaries.set_all(0)
force_face.mark(boundaries,1)
ds=Measure("ds")(subdomain_data=boundaries)

#define some continuum formulation
d=len(u_new)
I=Identity(d)

# Kinematics
d = len(u_new)
I = Identity(d)
# Identity tensor
def DeformationGradient(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)

def delta_E(v,u):
    F = DeformationGradient(u)
    return (dot(F.T,grad(v))+dot(grad(v).T,F))/2

def eigv_eigw_epsilon(e): #calculate the eigenvalue and eigenvector of a tensor
    tr_e  = tr(e)
    det_e = det(e)
    v1    = 0.5*tr_e + sqrt(tr_e*tr_e/4.0 - det_e)
    v2    = 0.5*tr_e - sqrt(tr_e*tr_e/4.0 - det_e)
    w1    = conditional(gt(abs(e[0,1]),1E-10), as_vector([v1-e[1,1], e[0,1]]), as_vector([1.0, 0.0]))
    w2    = conditional(gt(abs(e[0,1]),1E-10), as_vector([v2-e[1,1], e[0,1]]), as_vector([0.0, 1.0]))
    norm1 = sqrt(dot(w1,w1))
    norm2 = sqrt(dot(w2,w2))
    w1    = w1/norm1
    w2    = w2/norm2
    return v1, v2, w1, w2

#energy density of compressible neo-hookean solid
#W=mu*(I1-3-2ln(J))/2+0.5*K0*(ln(J))^2
# where I1=tr(F.T*F),J=det(F)

def pk2(u,K0,mu):  #calculate the 2nd piola-kirchhoff stress tensor
    C = RightCauchyGreen(u)
    J = Jacobian(u)
    C1=inv(C)
    eig12,eig22,w1,w2=eigv_eigw_epsilon(C)
    piola2=mu*(1-1/eig12)*outer(w1,w1)+mu*(1-1/eig22)*outer(w2,w2)+K0*ln(J)*C1
    return piola2

F_u=inner(pk2(u_new,K0,mu),delta_E(v,u_new))*dx
J_u=derivative(F_u,u_new,u)

# Variational problem for the displacement
problem_u = NonlinearVariationalProblem(F_u, u_new, bcu, J=J_u)
# Set up the solvers
solver_u  = NonlinearVariationalSolver(problem_u)

t=0
ur=1
deltaT=0.1
tol=1e-4
conc_f=File("./hyperelastic1/u.pvd")
fname=open('F_D_curve1.txt','w')
outfile=XDMFFile("output1.xdmf")

#staggered scheme
while t<=1.0:
    t+=deltaT
    if t>=0.4:
        deltaT=0.1
    applied_disp.disp=t*ur
    iter=0
    err=1
    while err>AM_tolerance and iter<maxiteration:
        iter +=1
        solver_u.solve()
        u_err = u_new.vector() - u_old.vector()
        err = u_err.norm('linf')
        u_old.assign(u_new)
        if err<tol:
            print('Iterations:',iter,',Total time',t)
            conc_f << u_new
            outfile.write(u_new, t)


print('simulation completed')
print(u_new.vector()[:])
1 Like