Preventing Newton-Solver Divergence for Nearly Incompressible Hyperelastic Materials

27 days ago by
I'm working with a nearly incompressible Neo-Hookean model with this potential energy function:

$W=\frac{1}{2}\left(mu\right)\left(tr\left(C_{hat}\right)-3\right)+\frac{1}{2}K\left(J-1\right)^2$W=12 (mu)(tr(Chat)3)+12 K(J1)2

When I impose any more than an extremely small deformation, the Newton solver diverges. I believe it is because, when it deforms the mesh even slightly, the incompressibility term makes the potential energy increase by a very large amount. I am currently using a K/mu ratio of 1000. How do I prevent this problem? I have experimented around with other nonlinear solvers, but I can't seem to get them to work as well as the Newton solver, and I need them to work better. I could try starting out with a better initial guess, but that would be pretty difficult to code, and I'd rather it be a last resort. Here's a snippet of my code:

body = Constant((0.0, 0.0, 0.0))
traction = Constant((0.0, 0.0, 0.0))

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

du = TrialFunction(V)
v = TestFunction(V)
u = Function(V)

d = u.geometric_dimension()
I = Identity(d)
F = variable(I + grad(u))
J = variable(det(F))
C = variable((F.T*F)*J**(-2.0/3.0))
I1 = variable(tr(C))

W = (mu/2)*(I1 - 3) + (1/2)*(K*((J - 1)**2))

Pi = W*dx - dot(body, u)*dx - dot(traction, u)*ds

F = derivative(Pi, u, v)
Jac = derivative(F, u, du)

problem = NonlinearVariationalProblem(F, u, dirichlet, J=Jac)
solver = NonlinearVariationalSolver(problem)
solver.parameters['nonlinear_solver'] = 'newton'
solver.parameters['newton_solver']['maximum_iterations'] = 20

Thank you so much for the help.
Community: FEniCS Project

1 Answer

27 days ago by
Hi, technically, what you do is completely correct. However, there are other material models giving a better numerical convergence. Try to change the model (and also use a stepwise increment of the load such that the last result is used as the initial guess). You may like to try the neo Hookean energy density:

\( W = \lambda\frac12 \ln^2(J) + \mu \frac12 ( F_{ij} F_{ij} - \delta_{ii} - \ln(J) ) \)

I know that you intend to use an incompressible model. The parameter lambda is basically the one like K controlling the incompressibility. It is difficult to see it there; but if you find out the stress:

\( P_{ji} = \frac{\partial W}{\partial F_{ij}} = ( \lambda \ln(J) - \mu ) F^{-1}_{ji} + \mu F_{ij} \)

then the role of lambda becomes clear. Details and implementation can be found in Section 1.3 of Abali, B. E. (2016). Computational Reality. Springer-Nature Singapore.
Thank you so much! I re-wrote my code so that it uses the stepwise increment that you suggested, and it now runs flawlessly with the original energy equation, even for larger deformations. Much appreciated.
written 27 days ago by jhsansom  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »