Cantilever plate

A plate of dimensions 100mm x 50mm x 0.5mm clamped at one end. It is subjected to concentrated forces of 5N at each corner at the tip.

A reference solution is computed with ABAQUS for the following shell elements:

Type	u_max	S_max
S3	30.36	423.6
S4	30.54	438.5
S4R	30.63	440.2

All reference values use the same mesh, and the plate is solved with both shell element types of torch-fem on those same nodes, so the triangles compare to S3 and the quadrilaterals to S4 and S4R.

import torch

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

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

# Material parameters
material = IsotropicElasticityPlaneStress(E=200000.0, nu=0.3)

Triangular shell elements

# Define nodes and elements of the plate
nodes, elements = rect_tri(10, 5, Lx=100.0, Ly=50.0, variant="zigzag")
nodes = torch.hstack([nodes, torch.zeros((nodes.size(0), 1))])

# Define Shell model
tri = Shell(nodes, elements, material, thickness=0.5)

# Define masks for boundaries
left = nodes[:, 0] < 0.01
right = nodes[:, 0] > 99.9
front = nodes[:, 1] < 0.01
back = nodes[:, 1] > 49.99

# Apply load and displacement boundary conditions
tri.forces[right & front, 2] = -5
tri.forces[right & back, 2] = -5
tri.constraints[left] = True

# Solve
u_tri, f, σ_tri, _, _ = tri.solve(aggregate_integration_points=False)

# Displacement on deformed configuration
tri.plot(
    u_tri[:, 0:3],
    node_property={"u": torch.linalg.norm(u_tri[:, 0:3], dim=1)},
    bcs=True,
)

v. Mises stress

# v. Mises stress at the top of the shell. Integration points run over the in-plane
# points and the thickness stations, of which the last one sits at +t/2.
s = σ_tri.reshape(-1, tri.n_z, tri.n_elem, 2, 2)[:, -1]
mises_tri = (
    torch.sqrt(
        s[..., 0, 0] ** 2
        - s[..., 0, 0] * s[..., 1, 1]
        + s[..., 1, 1] ** 2
        + 3 * s[..., 0, 1] ** 2
    )
    .max(dim=0)
    .values
)

tri.plot(
    u_tri[:, 0:3], element_property={"v. Mises stress (top)": mises_tri}, cmap="inferno"
)

Quadrilateral shell elements

# Define nodes and elements of the plate
nodes, elements = rect_quad(10, 5, Lx=100.0, Ly=50.0)
nodes = torch.hstack([nodes, torch.zeros((nodes.size(0), 1))])

# Define Shell model
quad = Shell(nodes, elements, material, thickness=0.5)

# Define masks for boundaries
left = nodes[:, 0] < 0.01
right = nodes[:, 0] > 99.9
front = nodes[:, 1] < 0.01
back = nodes[:, 1] > 49.99

# Apply load and displacement boundary conditions
quad.forces[right & front, 2] = -5
quad.forces[right & back, 2] = -5
quad.constraints[left] = True

# Solve
u_quad, f, σ_quad, _, _ = quad.solve(aggregate_integration_points=False)

# Displacement on deformed configuration
quad.plot(
    u_quad[:, 0:3],
    node_property={"u": torch.linalg.norm(u_quad[:, 0:3], dim=1)},
    bcs=True,
)

v. Mises stress

# v. Mises stress at the top of the shell. Integration points run over the in-plane
# points and the thickness stations, of which the last one sits at +t/2.
s = σ_quad.reshape(-1, quad.n_z, quad.n_elem, 2, 2)[:, -1]
mises_quad = (
    torch.sqrt(
        s[..., 0, 0] ** 2
        - s[..., 0, 0] * s[..., 1, 1]
        + s[..., 1, 1] ** 2
        + 3 * s[..., 0, 1] ** 2
    )
    .max(dim=0)
    .values
)

quad.plot(
    u_quad[:, 0:3],
    element_property={"v. Mises stress (top)": mises_quad},
    cmap="inferno",
)

Comparison

Both meshes carry the same nodes, so the quadrilaterals use half as many elements.

print(f"{'Type':<7} {'Elements':>9} {'u_max':>8} {'S_max':>8}")
print(
    f"{'Tria1':<7} {tri.n_elem:>9} "
    f"{torch.min(u_tri[:, 0:3]):8.2f} {torch.max(mises_tri):8.1f}"
)
print(
    f"{'Quad1':<7} {quad.n_elem:>9} "
    f"{torch.min(u_quad[:, 0:3]):8.2f} {torch.max(mises_quad):8.1f}"
)

Type     Elements    u_max    S_max
Tria1          72   -30.72    430.6
Quad1          36   -30.53    437.9