Center-loaded plate

A square plate of 100x100 is fully clamped at all edges and subjected to a concentrated force F=16.367 at the center. A symmetric quarter model is computed.

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

Type	u_max	S_max
S3	-3.773	72.60
S4	-3.987	89.90
S4R	-4.015	86.96

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)

# Constants
L = 100.0
F = -16.367

# Material parameters
mat = IsotropicElasticityPlaneStress(E=10000.0, nu=0.3)

Triangular shell elements

The symmetric quarter model is mirrored in x-direction and y-direction.

# Define nodes and elements of the plate
nodes, elements = rect_tri(9, 9, 0.5 * L, 0.5 * L, variant="zigzag")
nodes = torch.hstack([nodes, torch.zeros((len(nodes), 1))])

# Create Shell model
tri = Shell(nodes, elements, mat)

# Boundaries
top = nodes[:, 1] > (L / 2.0 - 0.01)
bottom = nodes[:, 1] < 0.01
left = nodes[:, 0] < 0.01
right = nodes[:, 0] > (L / 2.0 - 0.01)

# Force at the center, clamped outer edges and symmetry on the cut edges
tri.forces[0, 2] = F
tri.constraints[top] = True
tri.constraints[right] = True
tri.constraints[left, 0] = True
tri.constraints[left, 4] = True
tri.constraints[left, 5] = True
tri.constraints[bottom, 1] = True
tri.constraints[bottom, 3] = True
tri.constraints[bottom, 5] = True

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

# Displacement on deformed configuration, scaled by a factor of 10
tri.plot(
    10.0 * u_tri[:, 0:3],
    node_property={"u": torch.linalg.norm(u_tri[:, 0:3], dim=1)},
    mirror=(True, True, False),
    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",
    mirror=(True, True, False),
)

Quadrilateral shell elements

The symmetric quarter model is mirrored in x-direction and y-direction.

# Define nodes and elements of the plate
nodes, elements = rect_quad(9, 9, 0.5 * L, 0.5 * L)
nodes = torch.hstack([nodes, torch.zeros((len(nodes), 1))])

# Create Shell model
quad = Shell(nodes, elements, mat)

# Boundaries
top = nodes[:, 1] > (L / 2.0 - 0.01)
bottom = nodes[:, 1] < 0.01
left = nodes[:, 0] < 0.01
right = nodes[:, 0] > (L / 2.0 - 0.01)

# Force at the center, clamped outer edges and symmetry on the cut edges
quad.forces[0, 2] = F
quad.constraints[top] = True
quad.constraints[right] = True
quad.constraints[left, 0] = True
quad.constraints[left, 4] = True
quad.constraints[left, 5] = True
quad.constraints[bottom, 1] = True
quad.constraints[bottom, 3] = True
quad.constraints[bottom, 5] = True

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

# Displacement on deformed configuration, scaled by a factor of 10
quad.plot(
    10.0 * u_quad[:, 0:3],
    node_property={"u": torch.linalg.norm(u_quad[:, 0:3], dim=1)},
    mirror=(True, True, False),
    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",
    mirror=(True, True, False),
)

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} {u_tri.min():8.3f} {torch.max(mises_tri):8.2f}")
print(f"{'Quad1':<7} {quad.n_elem:>9} {u_quad.min():8.3f} {torch.max(mises_quad):8.2f}")

Type     Elements    u_max    S_max
Tria1         128   -3.999    78.74
Quad1          64   -3.986    89.76