Large stretch and compression

We compare different material models for uni-axial strain with a large stretch to twenty times the original length and a large compression to a twentieth of it.

import matplotlib.pyplot as plt
import numpy as np
import torch
from scipy.optimize import fsolve

from torchfem import Solid
from torchfem.materials import Hyperelastic3D, IsotropicElasticity3D
from torchfem.mesh import cube_hexa

torch.set_default_dtype(torch.float64)

# Material
E = 1000.0
nu = 0.3
lbd = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
mu = E / (2.0 * (1.0 + nu))

# Generate cube
nodes, elements = cube_hexa(5, 5, 5)

# Boundary sets
left = nodes[:, 0] == 0.0
right = nodes[:, 0] == 1.0
center1 = nodes[:, 1] == 0.5
center2 = nodes[:, 2] == 0.5

# Boundary conditions
constraints = torch.zeros_like(nodes, dtype=torch.bool)
constraints[left, 0] = True
constraints[right, 0] = True
constraints[center1, 1] = True
constraints[center2, 2] = True

# Load cases: displacement of the right face, applied in N increments
cases = {"Stretch": 19.0, "Compression": -0.95}
N = 21


def uniaxial(material, U):
    """Mean Cauchy stress σ_xx in each increment of a uniaxial strain to 1 + U."""
    box = Solid(nodes, elements, material)
    box.constraints = constraints
    box.displacements[right, 0] = U
    increments = torch.linspace(0.0, 1.0, N)
    _, _, σ, _, _ = box.solve(increments=increments, return_intermediate=True)
    return σ[..., 0, 0].mean(dim=1).detach()


def plot(ax, λ, σ_ref, σ, title):
    ax.plot(λ, σ_ref / E, ".-", color="gray", label="Reference")
    ax.plot(λ, σ / E, ".--", color="deeppink", label="torch-fem")
    ax.set_title(title)
    ax.set_xlabel("Stretch λ")
    ax.set_ylabel("Dimensionless Cauchy stress σ/E")
    ax.grid()
    ax.legend()

Linear elasticity with small strains

material = IsotropicElasticity3D(E=E, nu=nu)

σl = {}
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
for ax, (name, U) in zip(axes, cases.items()):
    λ = torch.linspace(1.0, 1.0 + U, N)
    σl[name] = uniaxial(material, U)
    plot(ax, λ, E * (λ - 1.0), σl[name], name)
plt.tight_layout()
plt.show()

Neo-Hookean hyperelasticity

The strain energy function for this Neo-Hookean material is chosen as

$$ \psi(C) = \frac{\mu}{2} \left( \textrm{tr}(C) - 3\right) - \mu \log(J) + \frac{\lambda}{2} \log(J)^2 $$

with the right Cauchy-Green tensor $$ C = F^\top F $$

and the Jacobian $$ J = \det(F). $$

def psi(F, params):
    """Neo-Hookean strain energy density function."""
    # Extract material parameters
    MU = params[0]
    LBD = params[1]
    # Compute the right Cauchy-Green deformation tensor
    C = F.transpose(-1, -2) @ F
    # Stable computation of the logarithm of the determinant
    logJ = 0.5 * torch.logdet(C)
    return MU / 2 * (torch.trace(C) - 3.0) - MU * logJ + LBD / 2 * logJ**2


def reference(λ):
    """Analytical Cauchy stress of the Neo-Hookean material under uniaxial strain."""
    λ = λ.numpy()
    J = fsolve(
        lambda J: mu * (J / λ - 1) + lbd * np.log(J), x0=np.ones(λ.shape), xtol=1e-10
    )
    J = np.asarray(J)
    return 1 / J * (lbd * np.log(J) + mu * (λ**2 - 1))


material = Hyperelastic3D(psi, params=[mu, lbd])

σn = {}
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
for ax, (name, U) in zip(axes, cases.items()):
    λ = torch.linspace(1.0, 1.0 + U, N)
    σn[name] = uniaxial(material, U)
    plot(ax, λ, reference(λ), σn[name], name)
plt.tight_layout()
plt.show()