Finite strain

We compare linear and non-linear geometric models for a cantilever beam subjected to a concentrated force at its tip.

import torch
import pyvista
import matplotlib.pyplot as plt

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

# Set default data type to double precision
torch.set_default_dtype(torch.float64)

Create the beam model

# Define linear elastic small strain material
material = IsotropicElasticity3D(E=1000.0, nu=0.3)

# Generate beam
L = 5.0
nodes, elements = cube_hexa(21, 5, 5, L, 1.0, 1.0)
beam = Solid(nodes, elements, material)

# Assign boundary conditions
P = -10.0
tip = (nodes[:, 0] == L) & (nodes[:, 1] == 0.5) & (nodes[:, 2] == 0.5)
beam.constraints[nodes[:, 0] == 0, :] = True
beam.forces[tip, 2] = P

# Increments
increments = torch.linspace(0.0, 1.0, 21)

Small strain isotropic elasticity without geometric non-linearity

# Solve
u_lg, _, σ_lg, F_lg, _ = beam.solve(increments=increments, return_intermediate=True)

# Show result
beam.plot(
    u=u_lg[-1],
    element_property={"Cauchy stress": σ_lg[-1, :, 0, 0]},
    show_undeformed=True,
)

Logarithmic strain isotropic elasticity with geometric non-linearity

# Change material model
beam.material = IsotropicHencky3D(E=1000.0, nu=0.3).vectorize(beam.n_elem)

# Solve
u_ng, _, σ_ng, F_ng, _ = beam.solve(
    increments=increments, return_intermediate=True, nlgeom=True
)

# Show result
beam.plot(
    u=u_ng[-1],
    element_property={"Cauchy stress": σ_ng[-1, :, 0, 0]},
    show_undeformed=True,
)

Reference solutions

# Reference time increments
increments_abq = torch.linspace(0.0, 1.0, 11)

# ABAQUS displacements at the tip with linear elastic material and NLGEOM=NO
u_abq_lg = torch.tensor(
    [
        [0.00e-08, 0, -0.0000],
        [0.63e-08, 0, -0.5097],
        [1.26e-08, 0, -1.0195],
        [1.89e-08, 0, -1.5293],
        [2.52e-08, 0, -2.0391],
        [3.16e-08, 0, -2.5489],
        [3.79e-08, 0, -3.0587],
        [4.42e-08, 0, -3.5685],
        [5.05e-08, 0, -4.0782],
        [5.68e-08, 0, -4.5880],
        [6.32e-08, 0, -5.0978],
    ]
)
# ABAQUS displacements at the tip with linear elastic material and NLGEOM=YES
u_abq_ng = torch.tensor(
    [
        [-0.000, 0, -0.00000],
        [-0.029, 0, -0.50353],
        [-0.111, 0, -0.98141],
        [-0.232, 0, -1.41558],
        [-0.377, 0, -1.79791],
        [-0.536, 0, -2.12832],
        [-0.697, 0, -2.41129],
        [-0.856, 0, -2.65328],
        [-1.008, 0, -2.86087],
        [-1.152, 0, -3.03998],
        [-1.288, 0, -3.19567],
    ]
)

fig, ax = plt.subplots(1, 2, figsize=(10, 4))
ax[0].plot(increments_abq, u_abq_lg[:, 0], "-", color="lightgray", label="ABAQUS")
ax[0].plot(increments_abq, u_abq_ng[:, 0], "-", color="gray", label="ABAQUS (NLGEOM)")
ax[0].plot(increments, u_lg[:, tip, 0], ".-", color="deeppink", label="torch-fem")
ax[0].plot(increments, u_ng[:, tip, 0], ".-", color="blue", label="torch-fem (nlgeom)")
ax[0].set_xlabel("Time")
ax[0].set_ylabel("Displacement X")
ax[0].legend()
ax[0].grid()
ax[1].plot(increments_abq, u_abq_lg[:, 2], "-", color="lightgray", label="ABAQUS")
ax[1].plot(increments_abq, u_abq_ng[:, 2], "-", color="gray", label="ABAQUS (NLGEOM)")
ax[1].plot(increments, u_lg[:, tip, 2], ".-", color="deeppink", label="torch-fem")
ax[1].plot(increments, u_ng[:, tip, 2], ".-", color="blue", label="torch-fem (nlgeom)")
ax[1].set_xlabel("Time")
ax[1].set_ylabel("Displacement Z")
ax[1].legend()
ax[1].grid()
plt.show()

Comparison

plotter = pyvista.Plotter()
beam.plot(color="gray", plotter=plotter, opacity=0.1)
beam.plot(u=u_lg[-1], color="cyan", opacity=0.8, plotter=plotter)
beam.plot(u=u_ng[-1], color="orange", opacity=0.8, plotter=plotter)
plotter.show(jupyter_backend="html")