### Preventing Newton-Solver Divergence for Nearly Incompressible Hyperelastic Materials

78

views

1

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$

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:

Thank you so much for the help.

$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 (`m``u`)(`t``r`(`C`_{hat})−3)+12`K`(`J`−1)^{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 `

`solver.solve()`

Thank you so much for the help.

Community: FEniCS Project

### 1 Answer

3

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.

\( 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.

written
27 days ago by
jhsansom

Please login to add an answer/comment or follow this question.