Problem with simple 1D nonlinear elasticity problem

Question

Hi,

I am testing out fenics to see if it would be ok for a non-mathematical PhD student of mine. But Ive run into a few problems with a simple 1D example.

The problem is to find the displacement of a beam $u$ of length $L$ help at one end (x=0), lying on a nonlinear foundation and placed under pressure. The potential energy is the strain energy elastic strain energy minus the work done by the load. i.e

$$W = \frac{1}{3} \int_0^L u^3\;dx - q \int_0^Lu\;dx$$

So find $u$ which minimises $W$ such that $u(0) = 0$ which minimises $W$. Adapted from the hyper elasticity demo here is my code.

from dolfin import *

import math
import ufl

import numpy as np

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

# User inputs
nel = 10

q = 0.0001;


# Create mesh and define function space
mesh = IntervalMesh(nel,0.0,1.0)

V = FunctionSpace(mesh, "Lagrange", 1) # Linear Lagrange Fintie Elements on 1D Mesh

# Mark Boundaries
left =  CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)

# Define Dirichlet boundary (u = 0 at x = 0)

bcl = DirichletBC(V, Constant(0.0), left)

# Define functions
du = TrialFunction(V)            # Incremental Displacement
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement from previous iteration

# Loading Looping


# Total Potential Energy
W =  (1.0/3.0) * (u ** 3) * dx(mesh) - q * u * dx(mesh)

# Compute first variation of Pi (directional derivative about u in the direction of v)
F = derivative(W, u, v)

# Compute Jacobian of F
J = derivative(F, u, du)

# Solve variational problem
solve(F == 0, u, bcl, J = J, form_compiler_parameters=ffc_options)

# Save solution in VTK format
file = File("u.pvd");
file << u;

Newton does not converge, with the following error

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at

---

*** fenics-support@googlegroups.com

---

*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.

---

*** -------------------------------------------------------------------------
*** 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: 2016.1.0
*** Git changeset: f4e22e59622b42182d16059f0212ddb1c6aa8712
*** -------------------------------------------------------------------------

If I make the problem linear i.e.

$$W = \frac{1}{3} \int_0^L u^2\;dx - q \int_0^Lu\;dx$$

Works fine (still thinks it is nonlinear)

If I put the 'power' of the foundation as only slightly nonlinear i.e 2.000001, it does not converge.

What have I missed!?

Thanks for you help in advance

Tim Dodwell (University of Exeter, UK)

Answer 1

Hi Tim, That is strange. Perhaps the problem is that u=0 is also a solution with those boundary conditions? With q very small the difference between u=0 and a nonzero solution may be too small, especially with such a coarse grid. Maybe you can scale the equation so that the nonzero solution is O(1)?

Answer 2

The functional has no minimum, but it has highly mesh-dependent local minima. Maybe you are missing something in the problem formulation?

The error here is misreported by dolfin - the Newton solver fails on the first iteration because the Hessian is singular at the initial guess. If you start with a positive initial guess it may converge to a local minimum.

Comments

Just to clear this up, this is correct.

Solved by adding tension to the model with is a term in the variational form of

\frac{1}{2}EA\int_0^Lu'^2\;

and works fine.

Thanks for you answer.