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()