Modal analysis of a fully clamped shell

This example computes natural frequencies and mode shapes of a flat fully clamped shell with the Shell model and solve_modes().

import torch

from torchfem import Shell
from torchfem.materials import IsotropicElasticityPlaneStress
from torchfem.mesh import rect_tri

torch.set_default_dtype(torch.float64)

Model setup

# Geometry
L = 100.0
W = 100.0
t = 1.0

# Material
E = 1000.0
nu = 0.3
rho = 1.0

material = IsotropicElasticityPlaneStress(E=E, nu=nu, rho=rho)
nodes, elements = rect_tri(21, 21, Lx=L, Ly=W, variant="zigzag")
nodes = torch.hstack([nodes, torch.zeros((nodes.size(0), 1))])

model = Shell(nodes, elements, material, thickness=t)
left = nodes[:, 0] == 0.0
model.constraints[left, :] = True
right = nodes[:, 0] == L
model.constraints[right, :] = True
top = nodes[:, 1] == W
model.constraints[top, :] = True
bottom = nodes[:, 1] == 0.0
model.constraints[bottom, :] = True

model.plot(bcs=True)

Modal analysis

Solve for the first four modes and print squared angular frequencies $\omega^2$.

N = 4
omega_sq, modes = model.solve_modes(n_modes=N)

Comparison to ABAQUS reference

The ABAQUS reference solution was obtained with S3 elements.

# ABAQUS reference values for the first 4 modes
omega_ref = [0.001226, 0.005194, 0.005194, 0.01188]

header = f"{'Mode':>4}  {'ω² (torch-fem)':>10}  {'ω² (ABAQUS)':>10}"
print(header)
print("-" * len(header))
for k in range(N):
    fem_val = float(omega_sq[k])
    abq_val = omega_ref[k]
    print(f"{k + 1:>4}  {fem_val:>10.4f}  {abq_val:>10.4f}")

Mode  ω² (torch-fem)  ω² (ABAQUS)
---------------------------------
   1      0.0012      0.0012
   2      0.0051      0.0052
   3      0.0051      0.0052
   4      0.0113      0.0119

Visualize mode shapes

for k in range(N):
    mode_k = modes[k, :, 0:3]
    mag = torch.linalg.norm(mode_k, dim=1).max()
    u_plot = 20 * mode_k / mag
    model.plot(
        u=u_plot,
        node_property={f"Amplitude_mode_{k + 1}": torch.linalg.norm(mode_k, dim=1)},
    )