Modal analysis of a solid cantilever beam¶

Open In Colab Binder

This example computes the natural frequencies and mode shapes of a 3-D cantilever beam using Solid.solve_modes().

In [1]:
import torch

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

torch.set_default_dtype(torch.float64)

Model setup¶

We build a cantilever beam with length $L = 10$ and square cross-section $H = W = 0.5$, clamped at one end.

In [2]:
# Geometric parameters
L = 10.0
H = 0.5
W = 0.5

# Material parameters
E = 1000.0
nu = 0.3
rho = 1.0

# Build model
material = IsotropicElasticity3D(E=E, nu=nu, rho=rho)
nodes, elements = cube_hexa(81, 5, 5, Lx=L, Ly=H, Lz=W)

model = Solid(nodes, elements, material)
left = nodes[:, 0] == 0.0
model.constraints[left, :] = True

model.plot()

Modal analysis¶

Solve for the first four natural frequencies and mode shapes. The ABAQUS reference uses C3D8 elements.

Note on the expected discrepancy. Abaqus C3D8 unconditionally applies selective reduced integration with $\bar{B}$: the volumetric part of $\mathbf{K}$ is computed with a single centroid Gauss point, while the deviatoric part uses the full 2×2×2 rule. torch-fem's Hexa1 uses full 2×2×2 Gauss for both parts, producing a slightly stiffer volumetric response and therefore higher frequencies on coarse meshes (~12% in $\omega^2$ for this 30×2×2 mesh). The two formulations converge to the same analytical answer under mesh refinement.

In [3]:
N = 4
omega_sq, modes = model.solve_modes(n_modes=N)
In [4]:
# ABAQUS reference values for the first 4 modes
omega_sq_ref = [0.02595, 0.02595, 0.99589, 0.99589]

header = f"{'Mode':>4}  {'ω² (torch-fem)':>14}  {'ω² (ABAQUS)':>12}"
print(header)
print("-" * len(header))
for k in range(N):
    print(f"{k + 1:>4}  {omega_sq[k]:>14.6f}  {omega_sq_ref[k]:>12.6f}")

Mode  ω² (torch-fem)   ω² (ABAQUS)
----------------------------------
   1        0.026699      0.025950
   2        0.026699      0.025950
   3        1.025721      0.995890
   4        1.025721      0.995890

Visualise mode shapes¶

The plotting call opens an interactive 3-D viewer for each mode.

In [5]:
for k in range(N):
    mode_k = modes[k]
    max_mag = torch.norm(mode_k, dim=1).max()
    u_plot = mode_k / max_mag
    model.plot(
        u=u_plot,
        node_property={"|u|": torch.norm(u_plot, dim=1)},
        show_undeformed=True,
        cmap="viridis",
    )