### How to correctly pull back applied external load in hyperelasticity?

174
views
0
3 months ago by
Dear all,

I am trying to impose some surface load (Neumann boundary condition) to a hyperelastic problem. Compared to examples I found online, I am given the applied surface load in the deformed configuration and need to pull it back to the reference configuration. Here is a slightly modified hyperelasticity demo:

import dolfin

dolfin.parameters["form_compiler"]["representation"] = "uflacs"
dolfin.parameters["form_compiler"]["cpp_optimize"] = True

mesh = dolfin.UnitSquareMesh.create(10, 10, dolfin.CellType.Type_triangle)
# mesh = dolfin.UnitSquareMesh.create(10, 10, dolfin.CellType.Type_quadrilateral)

mesh_normals = dolfin.FacetNormal(mesh)

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

left   = dolfin.CompiledSubDomain("near(x[0], x0) && on_boundary", x0=0.)
right  = dolfin.CompiledSubDomain("near(x[0], x1) && on_boundary", x1=1.)
bottom = dolfin.CompiledSubDomain("near(x[1], y0) && on_boundary", y0=0.)
top    = dolfin.CompiledSubDomain("near(x[1], y1) && on_boundary", y1=1.)

boundaries = dolfin.MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)
bottom.mark(boundaries, 3)
top.mark(boundaries, 4)

dx = dolfin.Measure("dx", mesh)
ds = dolfin.Measure("ds", mesh, subdomain_data=boundaries)

bcl = dolfin.DirichletBC(V.sub(0), 0., boundaries, 1)
bcb = dolfin.DirichletBC(V.sub(1), 0., boundaries, 3)
bcs = [bcl, bcb]

U  = dolfin.Function(V)
dU = dolfin.TrialFunction(V)
dV = dolfin.TestFunction(V)

dim = U.geometric_dimension()
I = dolfin.Identity(dim)
J = dolfin.det(F)
J0 = dolfin.det(I)
C = dolfin.transpose(F) * F
IC = dolfin.tr(C)
IC0 = dolfin.tr(I)

E, nu = 1.0, 0.3
mu = dolfin.Constant(E/(2*(1 + nu)))
lmbda = dolfin.Constant(E*nu/((1 + nu)*(1 - (dim-1)*nu)))

Psi = (mu/2) * (IC - IC0 - 2*dolfin.ln(J)) + (lmbda/2) * (J**2 - J0**2 - 2*dolfin.ln(J))
Pi = Psi * dx

T = dolfin.Constant((0.1, 0.0))

# FmTN = dolfin.dot(dolfin.transpose(dolfin.inv(F)), mesh_normals)
# T = J * dolfin.sqrt(dolfin.inner(FmTN, FmTN)) * T

Pi -= dolfin.inner(T, U) * ds(2)

F = dolfin.derivative(Pi, U, dV)
J = dolfin.derivative(F, U, dU)
dolfin.solve(F==0, U, bcs, J=J)

file = dolfin.File("displacement.pvd")
file << U
​

This one works, either with quadrilaterals:

or triangles:

Now, if you uncomment the two lines defining FmTN and T toward the bottom, the solution seems perturbed at the upper right corner, both with quads:

and triangles:

Are there some precautions to take in computing subdomains or mesh normals toward edges? Any help would be super appreciated.

Martin
Community: FEniCS Project

3
3 months ago by
First, with quads in 2017.2, many non-trivial boundary integrals produce incorrect results, so I would avoid quads/hexes in problems like this for now.  (Based on comments from other treads, I believe the developers are already aware of this, and working on solutions.)

In the implementation with triangles, it looks like the pulled-back traction is correctly derived from Nanson's formula, but you want to include it directly in the momentum balance residual rather than the energy functional.  The modifications below produce a displacement solution more like what you are probably expecting:

####### (Everything the same above here) #######

T = dolfin.Constant((0.1, 0.0))

FmTN = dolfin.dot(dolfin.transpose(dolfin.inv(F)), mesh_normals)
T = J * dolfin.sqrt(dolfin.inner(FmTN, FmTN)) * T

####### Removed #######
#Pi -= dolfin.inner(T,U) * ds(2)
#######################

F = dolfin.derivative(Pi, U, dV)

​