 # Hyperelastic Stress calculation

#1

Dear all,

I am trying to calculate Second Piola-Kirchhoff stress for a hyperelastic neo-Hookean material. However, my results are not correct. It is an incremental displacement problem for a 3D block where bottom surface is fixed and incremental displacement from the top. Any recommendation or suggestion would be greatly appreciated. Here is my complete code:

from dolfin import *
from fenics import *
import matplotlib.pyplot as plt
import numpy as np

#Initialize Mesh and Geometry
lx = 10
ly = 10
lz = 5
mesh = BoxMesh(Point(0,0,0), Point(lx, ly, lz),10, 10, 10)
V = VectorFunctionSpace(mesh, “Lagrange”, 1)
W = TensorFunctionSpace(mesh, ‘DG’, 0)
Z = FunctionSpace(mesh, ‘DG’, 0)
File(“my_mesh.xml.gz”) << mesh
plt.figure(1)
plot(mesh, title = “Original mesh”)
#plt.show()
#print("num elem: ", mesh.num_cells())

def main():
# Elasticity parameters
lmbda = 830
mu=50
du = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
u = Function(V) # Displacement from previous iteration

``````top   =   CompiledSubDomain("near(x, side) && on_boundary", side = lz)
bottom   =   CompiledSubDomain("near(x, side) && on_boundary", side = 0.0)

#Define surface boundary
boundaries = MeshFunction("size_t", mesh, 2)
boundaries.set_all(0)
top.mark(boundaries, 3)
bottom.mark(boundaries, 4)

# Surface force
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

disp = Expression(("-1.0"), degree=1)
bcs10 = DirichletBC(V.sub(2), disp , boundaries, 3)
bcs1 = DirichletBC(V, Constant([0.0, 0.0, 0.0]), boundaries, 4)
bcs = [bcs1, bcs10]

def Energyfunction(C):
I = Identity(3)             # Identity tensor dimension 3
F = I + grad(u)             # Deformation gradient
C = variable(F.T*F)                   # Right Cauchy-Green tensor
J  = det(F)
I1 = tr(C)
I2 = 0.5*(tr(C)**2-tr(C*C))
I1bar = J**(-2.0/3.0)*I1
psiV = (lmbda/2)*(J - 1)**2
psiD = (mu/2)*(I1bar - 3)
psi = psiV+psiD
return (psi, C)

psi, C=Energyfunction(u)
# Total potential energy
Pi = psi*dx
# Compute first variation of Pi (directional derivative about u in the direction of v)
FF = derivative(Pi, u, v)
# Compute Jacobian of F
J = derivative(FF, u, du)

solver_parameters = {"nonlinear_solver": "snes",
"snes_solver":
{"linear_solver": "lu",
"line_search": "basic",
"maximum_iterations": 30,
"relative_tolerance": 1e-6,
"error_on_nonconvergence": False}}

fid = XDMFFile("Solution_3d_incD/u.xdmf")
#fidS = XDMFFile("Solution_3d_incD/S.xdmf")
solve(FF == 0, u, bcs, J=J, bcs=bcs, solver_parameters=solver_parameters)
fid.write(u)
top = (lx/2.0,ly/2.0,lz)
top_disp = u(top)
print('Top Mid Displacement of the beam: ', top_disp)

fidS = XDMFFile("Solution_3d_incD/S.xdmf")
S=2*diff(psi,C)
S_project=project(S,W)
print('Stress: ', S_project.vector().get_local().min())
fidS.write(S_project)
``````

if name == “main”:
main()