Hyperelasticity Stress Calculation

Question

Hi, I am doing some test runs with Fenics demo program on hyperelasticity. The displacement is easily obtained but I need to calculate the stress. According to the hyperelasticity theory, stress is obtained by differentiation of an energy function with strain. Given the displacement u, my question is about calculation of stress (say PK2). Instead of a lengthy expression, I'd like to have the energy function differentiated with strain tensor C directly. See PK2 function below. Any help would be appreciated.

from dolfin import *

# Create mesh and define function space
mesh = UnitSquareMesh(16, 16)
Sc = FunctionSpace(mesh,"Lagrange",1)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
TT = TensorFunctionSpace(mesh,'DG',0)

# Mark boundary subdomians
left =  CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)

# Define Dirichlet boundary (x = 0 or x = 1)
c = Expression(("0.0", "0.0"))
r = Expression(("1.0", "0.0"))

bcl = DirichletBC(V, c, left)
bcr = DirichletBC(V, r, right)
bcs = [bcl, bcr]

# Define functions
du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement from previous iteration
B  = Constant((0.0, 0.0))  # Body force per unit volume
T  = Constant((0.0, 0.0))  # Traction force on the boundary

# Kinematics
d = u.geometric_dimension()
I = Identity(d)             # Identity tensor
F = I + grad(u)             # Deformation gradient
C = F.T*F                   # Right Cauchy-Green tensor

E = (C - I)/2               # Euler-Lagrange strain tensor
E = variable(E)

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

# Elasticity parameters
E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(Ic - 2) - mu*ln(J) + (lmbda/2)*(ln(J))**2

C = variable(C)
S = 2.0*diff(psi,C)

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

# Compute first variation of Pi (directional derivative about u in the direction of v)
F = derivative(Pi, u, v)

# Compute Jacobian of F
J = derivative(F, u, du)

# Solve variational problem
solve(F == 0, u, bcs, J=J,
      form_compiler_parameters=ffc_options)

# Define Green-Lagrange Finite Strain Tensor
def GLFS(u):
    F = I +grad(u)
    strain = 0.5*(F.T*F - I)    
    return strain

def PK2(u): #Piola-Kirchhoff Stress
    F = I +grad(u)
    C = 0.5*(F.T*F - I) 
    Ic = tr(C)
    J  = det(F)
    psi = (mu/2)*(Ic - 2) - mu*ln(J) + (lmbda/2)*(ln(J))**2
    C = variable(C)
    Ss = 2.0*diff(psi,C) 
    return Ss

PK2Tensor = project(PK2(u),TensorFunctionSpace(mesh,'DG',0)) 

Comments

Demo_Hyperelasticity

Hi Dave,
I had previously run the hyperelasticity demo with success. When I ran now, it produced some errors. I looked into the Forum and found that you had raised some question regarding the same demo on a different issue. The two statements that produced the errors for me are given below.

left, right = compile_subdomains(["(std::abs(x[0]) < DOLFIN_EPS) && on_boundary",
"(std::abs(x[0] - 1.0) < DOLFIN_EPS) && on_boundary"])
...................................................................................

Kinematics

I = Identity(V.cell().d) # Identity tensor``

I replaced the above, from your code sample and it was successful as follows:

left = CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)

Kinematics

d = u.geometric_dimension()
I = Identity(d) # Identity tensor

My question is though I knew that certain statements get updated from time to time in FEniCS, I do not know, as a beginner, where to look for? Even at present the code on the web for demo_hyperelasticity is the previous version which I used.

Any guidance will be appreciated for such situations in general.

Thanks

---

Hi Dave,

Did you ever figure that out? Thanks!

Martin

Answer 1

Hi!
You define:

E = (C - I)/2               # Euler-Lagrange strain tensor

and

E = 10.0

with the same name.
Could that be the problem?