# Hyperelastic beam

Hi, I’m solving the case of a 2d Mooney-Rivlin beam subjected to a distributed shear on the free edge x=6. I report my code that works correctly.

 from dolfin import *
import fenics as fe
import time

start = time.time()

parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}

mesh = RectangleMesh(Point(0.0, 0.0), Point(6.0, 1.5), 250, 50, "crossed")

V = VectorFunctionSpace(mesh, "Lagrange", 2)
left = CompiledSubDomain("near(x[0], side) && on_boundary", side=0.0, tol=10e-15)  # Code C++

c = Expression(("0.0", "0.0"), degree=2)
t = Expression(("0.0", "-1.0"), degree=2)

du = TrialFunction(V)            # Incremental displacement
v = TestFunction(V)             # Test function
u = Function(V)                 # Displacement from previous iteration

D = mesh.topology().dim()
print(D)
neumann_domain = MeshFunction("size_t", mesh, D-1)
neumann_domain.set_all(0)
CompiledSubDomain("near(x[0], side) && on_boundary", side=6.0, tol=10e-15).mark(neumann_domain, 1)
ds = Measure("ds", subdomain_data=neumann_domain)

bcs = DirichletBC(V, c, left)

d = u.geometric_dimension()
I = Identity(d)             # Identity tensor
C = F.T*F                   # Right Cauchy-Green tensor

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

# Elasticity parameters
J = sqrt(det(C))
I1 = tr(C)
I2 = 0.5 * (I1**2 - tr(C*C))
c1 = Constant(630)
c2 = Constant(-1.2)
c = Constant(100)
d = 2 * (c1 + 2 * c2)

#MR
psi = c*(J-1)**2 - d*ln(J) + c1*(I1-3) + c2*(I2-3)

# E, nu = 10**3, 0.3
# mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

# Neo-Hookean
# psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2

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

F = derivative(Pi, u, v)

J = derivative(F, u, du)
problem = NonlinearVariationalProblem(F, u, bcs, J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 1E-6
prm['newton_solver']['relative_tolerance'] = 1E-6
# prm['newton_solver']['maximum_iterations'] = 25
# prm['newton_solver']['relaxation_parameter'] = 1.0

prm['newton_solver']['linear_solver'] = 'petsc'

solver.solve()

L2 = inner(u, u) * dx
energynorm = sqrt(assemble(psi*dx))
# H1 = inner(F, P) * dx
L2norm = sqrt(assemble(L2))
H10norm = sqrt(assemble(H10))

print("L2 norm = %.10f" % norm(u, norm_type="L2"))
print("H1 norm = %.10f" % norm(u, norm_type="H1"))
print("H10 norm = %.10f" % norm(u, norm_type="H10"))
print("Running time = %.3f" % float(time.time()-start))


I decide to plot the Cauchy stress tensor, using this lines of code:

import matplotlib.pyplot as plt
C = F.T * F
J = sqrt(det(C))
I1 = tr(C)
I2 = 0.5 * (I1**2 - tr(C*C))
S = (2*c1 + 2*c2*I1)*I - 2*c2*C.T + (2*c*(J-1)*J - d)*inv(C.T)
s=1/J*F*S*F.T
pp=plot(s[1,1],mode='color',cmap='viridis')
plt.colorbar(pp)


That’s the component \sigma_{yy}:

I’m not able to understand if the non-zero stresses that I have in correspondence of the free edge x=6 are correct or they are due to an error in the code or a problem in the mesh (if I refine they remain).
Could you help me?
Thank you

You cannot expect that \sigma_{yy}=0 on this edge, you should rather have \sigma_{nn}=0 where n is the direction of the normal in the deformed configuration

Thank you very much.

One more thing, referring to what you say, is it correct to say that the component xx of the first Piola-Kirchhoff stress tensor is zero on the right side (since its the normal stress in the reference configuration)?

Yes this is correct.

Thank you very much.