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

174

views

0

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:

This one works, either with quadrilaterals:

Martin

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)
F = I + dolfin.grad(U)
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

### 1 Answer

3

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:

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)
####### Added #######
F -= dolfin.inner(T,dV)*ds(2)
#####################
J = dolfin.derivative(F, U, dU)
dolfin.solve(F==0, U, bcs, J=J)
file = dolfin.File("displacement.pvd")
file << U
```

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