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[2], side) && on_boundary", side = lz)
bottom   =   CompiledSubDomain("near(x[2], 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)[2]
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()