Hyperelasticity problem - mesh convergence issue

[Amirr]

I am following the tutorial for hyperelasticity problem. I am trying to simulate plain strain bending problem of a block which is constrained from the left side and a body force is applied. I am able to solve the problem and newton solver converges. However, as I refine the mesh, the value of the components of displacement gradient at the nodes close to boundary keeps increasing while tip displacement of the block remains almost constant.

Does anyone know how to get converged results?
Thank you

Here is my code:

from dolfin import *

# 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}

# Create mesh and define function space
mesh = RectangleMesh(Point(0,0),Point(6,1),900, 150)
V = VectorFunctionSpace(mesh, "Lagrange", 2)

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

bcl = DirichletBC(V,  Constant((0.,0.)) , left)
bcs = [bcl]

# Define functions
du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement from previous iteration
B  = Constant((0.0, -0.01))  # Body force per unit volume
T  = Constant((0.0,  0.0))  # Traction force on the boundary

# Kinematics
d = u.geometric_dimension()
I = Identity(d)             # Identity tensor
F = I + grad(u)             # Deformation gradient
C = F.T*F                   # Right Cauchy-Green tensor

# Invariants of deformation tensors
Ic = tr(C)
J  = det(F)

# Elasticity parameters
E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2

# Total potential energy
Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds

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

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

# Solve variational problem
solve(Func == 0, u, bcs, J=J,
      form_compiler_parameters=ffc_options)



W=FunctionSpace(mesh,"CG",2)
project_of_grad_u = project(grad(u)[1,1],W)
file = File('gradu11.pvd')
file << project_of_grad_u

[dokken]

Please format your code using ``` so that it is runnable by others

[Amirr]

Thanks for your comment, Now it is runnablel.

[bleyerj]

You should not project strain components (which are inherently living in a DG-1 space when coming from a CG-2 displacement) over a CG 2 space.

[dokken]

As @bleyer said, you need to use a DG1 space for the strain.
To properly visualize DG or higher order CG spaces in paraview, use XDMFFile and write_checkpoint as

W=FunctionSpace(mesh,"DG",1)
project_of_grad_u = project(grad(u)[1,1],W)

XDMFFile("gradu11.xdmf").write_checkpoint(project_of_grad_u,"gradu11", 0.0)

[Amirr]

Thanks for your reply, Now I am using DG1 but still I see values keep increasing as mesh size decreases. Does it mean I have to use finer mesh because I think my mesh size is quite small already and tip displacement is converged.
Here is the results with seed numbers 900*150:

And results with seed numbers 600*100:

Thanks

[dokken]

You should have a look at the development on even finer meshes, as well as having a look at the deformation field as well (The deformation is quite large).