# Incompressible Humphrey hyperelastic material in Fencis

Hi,

I am trying to implement Humphrey hyperelastic material (strain energy =` c1*(e^(c2*(I1-3))-1)` ) in Fenics. I am following the code given here (https://fenics-handson.readthedocs.io/en/latest/elasticity/doc.html). In this link, the stress has been derived as follows:

# Define constitutive law

``````def stress(u, p):
"""Returns 1st Piola-Kirchhoff stress and (local) mass balance
for given u, p."""
mu = Constant(E/(2.0*(1.0 + nu)))
J = det(F)
B = F * F.T
T = -p*I + mu*(B-I) # Cauchy stress
S = J*T*inv(F).T # 1st Piola-Kirchhoff stress
if nu == 0.5:
# Incompressible
pp = J-1.0
else:
# Compressible
lmbd = Constant(E*nu/((1.0 + nu)*(1.0 - 2.0*nu)))
pp = 1.0/lmbd*p + (J*J-1.0)
return S, pp
``````

If I want to implement incompressible Humphrey hyperelastic material, what modification should I make. I know for Humphrey material
Cauchy stress, `T = -p*I + 2*c1*c2*e^(c2*(I1-3)) * B`
But I donot know how to derive the hydrostatic pressure, pp, for this material. I tried the same expression , pp = J-1.0, but then I got error that â€śnon-linear solver fail to convergeâ€ť.

Could any one please let me know how to derive the hydrostatic pressure term?

I found the solution. There were a few issues with the variation formulation. Here is the complete code for uniaxial testing of incompressible NeoHookean hyperelastic material.

``````>
> from dolfin import *
> import matplotlib.pyplot as plt
> import os
> import dolfin as dlf
> import numpy as np
> import math
>
> parameters["form_compiler"]["cpp_optimize"] = True
> parameters["form_compiler"]["representation"] = "uflacs"
>
>
> # traction
> class Traction(UserExpression):
>     def __init__(self):
>         super().__init__(self)
>         self.t = 0.0
>     def eval(self, values, x):
>         values[0] = 0*self.t
>         values[1] = 0.0
>         values[2] = 0.0
>     def value_shape(self):
>         return (3,)
>
> # Kinematics
> def pk1Stress(u,pressure,E,nu):
>     G = E/(2*(1+nu))
>     c1 = G/2.0
>
>     I = Identity(V.mesh().geometry().dim())  # Identity tensor
>     C = F.T*F                # Right Cauchy-Green tensor
>     Ic = tr(C)               # Invariants of deformation tensors
>     J = det(F)
>     pk2 = 2*c1*I-pressure*inv(C) # second PK stress
>     return pk2, (J-1)
>
>
> def geometry_3d():
>     mesh = UnitCubeMesh(6, 6, 6)
>
>     boundary_parts = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
>     x0 = AutoSubDomain(lambda x: near(x[0], 0))
>     x1 = AutoSubDomain(lambda x: near(x[0], 1))
>     y0 = AutoSubDomain(lambda x: near(x[1], 0))
>     z0 = AutoSubDomain(lambda x: near(x[2], 0))
>     x0.mark(boundary_parts, 1)
>     y0.mark(boundary_parts, 2)
>     z0.mark(boundary_parts, 3)
>     x1.mark(boundary_parts, 4)
>     return boundary_parts
>
>
>
> # Create mesh and define function space ============================================
> facet_function = geometry_3d()
> mesh = facet_function.mesh()
> gdim = mesh.geometry().dim()
> dx = Measure("dx")
> ds = Measure("ds", subdomain_data=facet_function, subdomain_id=4)
> print('Number of nodes: ',mesh.num_vertices())
> print('Number of cells: ',mesh.num_cells())
>
>
> dx = dx(degree=4)
> ds = ds(degree=4)
>
>
> # Create function space
> element_v = VectorElement("P", mesh.ufl_cell(), 1) #
> element_s = FiniteElement("P", mesh.ufl_cell(), 1) #
> mixed_element = MixedElement([element_v, element_s]) #
> V = FunctionSpace(mesh, mixed_element)
>
>
> # Define test and trial functions
> dup = TrialFunction(V)
> _u, _p = TestFunctions(V)
>
> _u_p = Function(V)
> u, p = split(_u_p)
>
>
> # Create tensor function spaces for stress and stretch output
> W_DFnStress = TensorFunctionSpace(mesh, "DG", degree=0)
> PK1_stress = Function(W_DFnStress, name='PK1')
>
>
> # Displacement from previous iteration
> b = Constant((0.0,0.0, 0.0)) # Body force per unit mass
> h = Traction() # Traction force on the boundary
>
>
> # Define Dirichlet boundary
> bc0 = DirichletBC(V.sub(0).sub(0), Constant(0.), facet_function, 1)
> bc1 = DirichletBC(V.sub(0).sub(1), Constant(0.), facet_function, 2)
> bc2 = DirichletBC(V.sub(0).sub(2), Constant(0.), facet_function, 3)
>
> tDirBC = Expression(('1.0*time_'),time_=0.0 , degree=0)
> bc3 = DirichletBC(V.sub(0).sub(0), tDirBC, facet_function, 4)
> bcs = [bc0,bc1,bc2,bc3]
>
> # material parameters
> E, nu = 1e3, 0.5
>
> # # Total potential energy
> pkstrs, hydpress =  pk1Stress(u,p,E,nu)
> I = Identity(V.mesh().geometry().dim())
> dgF = I + grad(u)
> F1 = inner(dot(dgF,pkstrs), grad(_u))*dx - dot(b, _u)*dx - dot(h, _u)*ds
> F2 = hydpress*_p*dx
> F = F1+F2
> J = derivative(F, _u_p,dup)
>
> # Create nonlinear variational problem and solve
> problem = NonlinearVariationalProblem(F, _u_p, bcs=bcs, J=J)
> solver = NonlinearVariationalSolver(problem)
> solver.parameters['newton_solver']['relative_tolerance'] = 1e-6
> solver.parameters['newton_solver']['linear_solver'] = 'mumps'
>
>
> # Time stepping parameters
> dt = 0.1
> t, T = 0.0, 10*dt
>
> # Save solution in VTK format
> file_results.parameters["flush_output"] = True
> file_results.parameters["functions_share_mesh"] = True
>
> stretch_vec = []
> stress_vec=[]
> stress_ana=[]
>
> while t <= T:
>     print('time: ', t)
>
>     # Increase traction
>     h.t = t
>     tDirBC.time_ = t
>
>     # solve and save disp
>     solver.solve()
>
>     # Extract solution components
>     u_plot, p_plot = _u_p.split()
>     u_plot.rename("u", "displacement")
>     p_plot.rename("p", "pressure")
>
>     # get stretch at a point for plotting
>     point = (0.5,0.5,0)
>     DF = I + grad(u_plot)
>
>     # get stress at a point for plotting
>     PK1_s,thydpress = pk1Stress(u_plot,p_plot,E,nu)
>     PK1_stress.assign(project(PK1_s, W_DFnStress))
>     stress_vec.append(PK1_stress(point)[0])
>
>     # save xdmf file
>     file_results.write(u_plot, t)
>     file_results.write(PK1_stress, t)
>
>     # time increment
>     t += float(dt)
>
>
> # get analytical solution
> stretch_vec = np.array(stretch_vec)
> stress_vec = np.array(stress_vec)
> G = E/(2*(1+nu))
> c1 = G/2.0
> for i in range(len(stretch_vec)):
>     pk1_ana = 2*c1*(stretch_vec[i] - 1/(stretch_vec[i]*stretch_vec[i])) #PK1
>     pk2_ana = (1/stretch_vec[i])*pk1_ana # PK2
>     stress_ana.append(pk2_ana)
> stress_ana =  np.array(stress_ana)
>
> # plot results
> f = plt.figure(figsize=(12,6))
> plt.plot(stretch_vec, stress_vec,'r-')
> plt.plot(stretch_vec, stress_ana,'k.')

``````

Result:

1 Like

if i want to apply" tdirbs" on 2 different facets, meaning that two different facts with 2 different boundary conditions, i guess ill need to update the â€śdsâ€ť to be on those two subdomains right? any idea how can i do it?

for exp
bc2 = DirichletBC(V.sub(0), tDirBC, boundary_parts, 2)
bc4 = DirichletBC(V.sub(1), tDirBC, boundary_parts, 4)

and also, if i want to use two different â€śtDirBCâ€ť with different values, how do i consider that in my time loop?

and also, if i will change the the [0] to [1],[2] will i get the stress and stretch at y,z accordantly ?