[Legacy 2019.1.0] unexpected non-zero displacement in hyperplasticity problem

General
1

Hello,

Following the Hyperplasticity tutorial, with a change of potential energy function W on line 26 in the following example code. Even though the applied BCs are undeformed, the domain appears to be distorted (while the potential is zero), and everything seems to point to fenics.ln() , does any one know if there is any known issues/fixes regarding this ?

Figure_1
Figure_1640×480 18.7 KB

import fenics as fe
import matplotlib.pyplot as plt
q_degree = 5
fe.dx = fe.dx(metadata={'quadrature_degree': q_degree})
fe.set_log_level(fe.LogLevel.ERROR)

mesh = fe.UnitSquareMesh(10,10)

def Bottom(x, on_boundary):return (on_boundary and fe.near(x[1], 0.0))
def Top(x, on_boundary):return (on_boundary and fe.near(x[1], 1.0))

V = fe.VectorFunctionSpace(mesh, "CG", 2)
bcs = [ fe.DirichletBC(V, fe.Constant((0.0, 0.0)), Bottom),  fe.DirichletBC(V, fe.Constant((0.0, 0.0)), Top)]

del_u = fe.TrialFunction(V)
u_test = fe.TestFunction(V)
u = fe.Function(V)

I = fe.Identity(2)
F = I + fe.grad(u) 
F = fe.variable(F)
C = F.T*F
I_1 = fe.tr(C)+fe.Constant(1.)
J = fe.det(F)

W= 1. * (I_1 - 3 - 2 * fe.ln(J)) + 1. * fe.ln(J)
form = fe.derivative(W* fe.dx, u, u_test)
problem = fe.NonlinearVariationalProblem(form, u, bcs, fe.derivative(form, u, del_u))
solver = fe.NonlinearVariationalSolver(problem)
solver.solve()

plt.figure()
fig=fe.plot(u, mode="displacement",title='displacement')
plt.colorbar(fig)
plt.show()

Additionally by allowing log printing, the Newton solver appears to perform multiple iteration even though it is on undeformed configuration:

Solving nonlinear variational problem.
  Newton iteration 0: r (abs) = 3.300e-01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
  Newton iteration 1: r (abs) = 4.970e-02 (tol = 1.000e-10) r (rel) = 1.506e-01 (tol = 1.000e-09)
  Newton iteration 2: r (abs) = 1.176e-03 (tol = 1.000e-10) r (rel) = 3.563e-03 (tol = 1.000e-09)
  Newton iteration 3: r (abs) = 8.236e-07 (tol = 1.000e-10) r (rel) = 2.496e-06 (tol = 1.000e-09)
  Newton iteration 4: r (abs) = 3.872e-13 (tol = 1.000e-10) r (rel) = 1.173e-12 (tol = 1.000e-09)
  Newton solver finished in 4 iterations and 4 linear solver iterations.
2

That’s highly unlikely. It is much more likely that there is a bug in your code, or you are not solving the problem you expected.

Solid mechanics is not my field, but that is not the definition of the first principal invariant I remember. Why is there an addend 1? That will result in a forcing term, which will make the solution non-zero, and hence …

… should be totally expected.

3

Hi Francesco,

Thanks you for your comments! The added one in I_1 = fe.tr(C)+fe.Constant(1.) is due to the fact that the problem is 2D and the trace of C in undeformed configuration would be 2 instead of 3, this is reflected in the energy function with I_1 - 3. At undeformed configuration, I use print('W = ',fe.project(W).vector().get_local()) to print the values of W and they are all zero as expected.