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
S3	-3.773
STRI3	-4.041
STRI65	-4.486

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

# Constants
L = 100.0
F = -16.367

# Material parameters (Cs is transverse shear stiffness)
mat = IsotropicElasticityPlaneStress(E=10000.0, nu=0.3)

# Define 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
plate = 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
plate.forces[0, 2] = F
plate.constraints[top] = True
plate.constraints[right] = True
plate.constraints[left, 0] = True
plate.constraints[left, 4] = True
plate.constraints[left, 5] = True
plate.constraints[bottom, 1] = True
plate.constraints[bottom, 3] = True
plate.constraints[bottom, 5] = True

Solve

u, f, σ, _, _ = plate.solve(aggregate_integration_points=False)

# Extract stresses along thickness integration points
σ_top = σ[0]
σ_middle = σ[1]
σ_bottom = σ[2]

# Print maximum displacement
print(f"u_max: {u.min():.2f}")

u_max: -3.98

Visualize displacement

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

factor = 10.0
plate.plot(
    factor * u[:, 0:3],
    node_property={"u": torch.linalg.norm(u[:, 0:3], dim=1)},
    mirror=(True, True, False),
    bcs=True,
)

Visualize v. Mises stress

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

# Compute v. Mises stress at top (z=0.25, integration point 0)
mises_top = torch.sqrt(
    σ_top[:, 0, 0] ** 2
    - σ_top[:, 0, 0] * σ_top[:, 1, 1]
    + σ_top[:, 1, 1] ** 2
    + 3 * σ_top[:, 0, 1] ** 2
)

plate.plot(
    u[:, 0:3],
    element_property={"v. mises Stress (top)": mises_top},
    cmap="inferno",
    thickness=True,
    mirror=(True, True, False),
)