Adding a zeros-average condition on hyperelasticity example

variational formulation
1

Hi,
I would like to add a zero-average condition on the fluctuation field on the hyperelasticity demo, by considering an additional vectorial Lagrange multiplier λ. (as in this example neumann-poisson_demo )
But I doesn’t see how to change the derivative of potential energy in order to add the new average condition. Here is my code :

from dolfin import *

# The behavior of the form compiler FFC can be adjusted by prescribing
# various parameters. Here, we want to use the UFLACS backend of FFC::

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"

# Create mesh and define function space
mesh = UnitCubeMesh(24, 16, 16)

Re = VectorElement("R", mesh.ufl_cell(), 0)
Ve = VectorElement("Lagrange", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement([Ve,Re]))

w=Function(W)

(u,lamb) = split(w)

# Define functions
v,Tlamb  = TestFunctions(W)             # Test function
du, dlamb = TrialFunctions(W)

B  = Constant((0.0, 0.0, 0.0))  # Body force per unit volume
T  = Constant((0.0,  0.0, 0.0))  # Traction force on the boundary

# Kinematics
d = len(u)
I = Identity(d)             # Identity tensor
F = I + grad(u) + Constant([[0,0.2,0],[0.2,0,0],[0,0,0]]) # Deformation gradient + constant macro strain
C = F.T*F                   # Right Cauchy-Green tensor

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

# Elasticity parameters
E, nu = 1e6, 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)
F = derivative(Pi, u, v) + dot(Tlamb,du)*dx + dot(dlamb,v)*dx

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

# Solve variational problem
solve(F == 0, w, [], J=J)

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

The most concise approach would be to add a constraint term directly to Pi, then take the derivative for F with respect to the mixed solution w. To do this, you would define

#Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds
Pi = psi*dx - dot(B, u)*dx - dot(T, u)*ds + dot(lamb,u)*dx

then take derivatives

#F = derivative(Pi, u, v) + dot(Tlamb,du)*dx + dot(dlamb,v)*dx
F = derivative(Pi,w)

and

#J = derivative(F, u, du)
J = derivative(F, w)

If you want to specify a TestFunction Tw and TrialFunction dw as third arguments to derivative, as done with displacement in your MWE, you would need to make the modifications

#v,Tlamb  = TestFunctions(W)
Tw = TestFunction(W)
v,Tlamb = split(Tw)
#du, dlamb = TrialFunctions(W)
dw = TrialFunction(W)
du,dlamb = split(dw)

earlier in the script, but this is not technically necessary. An equivalent approach adding the constraint to F instead of Pi would be to add

+ dot(lamb,v)*dx + dot(Tlamb,u)*dx

i.e., using the Function w for the solution, not the TrialFunction dw, as done in the MWE. This is consistent with the rest of F, where the displacement is u, not du. There is no TrialFunction until differentiating again, for the Jacobian J.

For the output, you’ll also need to add u,lamb = w.split() before writing u to a file.

2 Likes