12  3D Linear Elasticity with P1 Tetrahedra

This notebook extends the 2D elasticity concepts to full 3D. We derive the element stiffness matrix K_e for a linear tetrahedron under isotropic elasticity, demonstrating:

  1. 3D Voigt notation for stress and strain (6 components instead of 3)
  2. Constitutive matrix D — the 6×6 isotropic matrix relating stress to strain
  3. Strain-displacement matrix B — how nodal displacements map to 6D Voigt strains
  4. Element stiffness K_e — the 12×12 matrix (3 DOFs × 4 nodes) from the integration K_e = V B^T D B
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import Math, display

# Import the 3D elasticity workflow
from symbolic_fem_workbench import build_elasticity_tetra_p1_3d
from symbolic_fem_workbench.elasticity import isotropic_3d_D

12.1 3D Constitutive Matrix

In 3D, stress and strain are both 6-component vectors under Voigt notation:

\[\boldsymbol{\sigma} = [\sigma_{xx}, \sigma_{yy}, \sigma_{zz}, \tau_{yz}, \tau_{xz}, \tau_{xy}]^T\] \[\boldsymbol{\varepsilon} = [\varepsilon_{xx}, \varepsilon_{yy}, \varepsilon_{zz}, \gamma_{yz}, \gamma_{xz}, \gamma_{xy}]^T\]

The isotropic constitutive matrix for linear elasticity is:

E, nu = sp.symbols("E nu", positive=True)
D = isotropic_3d_D(E, nu)

print("6×6 isotropic constitutive matrix:")
display(Math(r"D = " + sp.latex(D)))
print(f"\nMatrix shape: {D.shape}")

12.1.1 Material Parameters

The matrix is written in terms of: - E — Young’s modulus - ν — Poisson’s ratio

The off-diagonal shear blocks use the Lamé parameters: \[\lambda = \frac{E\nu}{(1+\nu)(1-2\nu)}, \quad \mu = \frac{E}{2(1+\nu)}\]

12.2 3D B-Matrix

The strain-displacement matrix B is 6×12 for a tetrahedron with 4 nodes (3 DOFs each). It maps nodal displacements u = [u₁, v₁, w₁, u₂, v₂, w₂, u₃, v₃, w₃, u₄, v₄, w₄]^T to Voigt strains:

\[\boldsymbol{\varepsilon} = \mathbf{B} \mathbf{u}\]

For P1 tetrahedra with constant shape function gradients, B is constant over the element.

# Build the complete 3D elasticity problem for a generic tetrahedron
data = build_elasticity_tetra_p1_3d()

print(f"B matrix shape: {data['B'].shape}")
print(f"\nFor the reference tetrahedron with vertices:")
for i, node in enumerate(data['reference_tetra'].nodes, 1):
    print(f"  Node {i}: {node}")

print(f"\n3D Jacobian matrix:")
display(Math(r"J = " + sp.latex(data["jacobian"])))

12.2.1 Volume Calculation

For an affine map, the volume is determined by the Jacobian:

print("Volume of tetrahedron:")
print(f"V = |det(J)| / 6 = {data['volume']}")

12.3 Element Stiffness K_e = V B^T D B

The element stiffness matrix is computed as:

\[\mathbf{K}_e = V \mathbf{B}^T \mathbf{D} \mathbf{B}\]

where V is the element volume. This is a 12×12 symmetric positive-definite matrix.

# Element stiffness for the unit right tetrahedron
Ke = data["Ke_unit_right_tetra"]

print(f"K_e shape: {Ke.shape}")
print(f"Symmetric: {sp.simplify(Ke - Ke.T) == sp.zeros(12, 12)}")

# Evaluate for steel-like material: E = 210 GPa, ν = 0.3
Ke_steel = Ke.subs({E: 210e9, nu: sp.Rational(3, 10)})
Ke_steel_num = np.array(Ke_steel, dtype=float)

print(f"\nK_e (E=210 GPa, ν=0.3, unit tet):")
print(f"  trace(K_e) = {np.trace(Ke_steel_num):.4e}")
print(f"  max|K_e| = {np.max(np.abs(Ke_steel_num)):.4e}")
print(f"  min eigenvalue = {np.min(np.linalg.eigvals(Ke_steel_num)):.4e}")

12.4 Symbolic Form

For the unit right tetrahedron with volume 1/6, the symbolic stiffness is:

# Show a compact form by reducing to essential parameters
print("Symbolic K_e (showing structure):")
print(f"Shape: {Ke.shape}")
print(f"First 3×3 block (DOFs at node 1):")
K11 = Ke[:3, :3]
display(Math(sp.latex(K11)))

12.5 Exercise

Compare how Poisson’s ratio affects the stiffness matrix. Plot the eigenvalues of K_e for different values of ν (say 0.0, 0.25, 0.49) and observe: - How the bulk modulus K = E / (3(1-2ν)) behaves as ν → 0.5 (incompressibility limit) - Which eigenvalues grow largest

# Exercise: Eigenvalue sensitivity to Poisson's ratio
nu_values = [0.0, 0.25, 0.3, 0.49]
E_val = 1.0  # Normalize to E = 1 for comparison

fig, ax = plt.subplots(figsize=(10, 6))

for nu_val in nu_values:
    Ke_num = np.array(Ke.subs({E: E_val, nu: nu_val}), dtype=float)
    eigs = np.sort(np.linalg.eigvals(Ke_num))[::-1]  # largest first
    ax.semilogy(range(1, len(eigs)+1), np.abs(eigs), marker='o', label=f'ν={nu_val:.2f}')

ax.set_xlabel('Eigenvalue index')
ax.set_ylabel('Eigenvalue magnitude')
ax.set_title('Eigenvalue Spectrum of K_e (unit tet, E=1) vs Poisson Ratio')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Print eigenvalue counts (rank)
print("Eigenvalue summary (E=1, unit right tet):")
for nu_val in nu_values:
    Ke_num = np.array(Ke.subs({E: E_val, nu: nu_val}), dtype=float)
    eigs = np.linalg.eigvals(Ke_num)
    rank = np.sum(np.abs(eigs) > 1e-10)
    print(f"  ν={nu_val:.2f}: rank={rank}, λ_max/λ_min={np.max(np.abs(eigs))/np.max(np.abs(eigs[eigs != 0])):.2e}")