Cantilever with a composite laminate

A 100 mm x 50 mm x 1 mm cantilever plate, clamped at one end and loaded by transverse forces at the two tip corners. The plate is made of a linear-elastic cross-ply laminate.

The Shell measures ply angles from a global orientation (projected onto each element), so we can sweep that reference direction to rotate the laminate's reference direction and study how the tip deflection depends on the fiber angle.

import math

import matplotlib.pyplot as plt
import torch

from torchfem import Laminate, Shell
from torchfem.materials import OrthotropicElasticityPlaneStress
from torchfem.mesh import rect_quad, rect_tri

torch.set_default_dtype(torch.float64)

Material and laminate

# Unidirectional GFRP
gfrp = OrthotropicElasticityPlaneStress(
    E_1=54000.0,
    E_2=9400.0,
    nu_12=0.33,
    G_12=5500.0,
    G_13=5500.0,
    G_23=3000.0,
    rho=1.8e-9,
)

# Cross-ply laminate
layup = Laminate(
    materials=[gfrp] * 4,
    thicknesses=[0.25] * 4,
    angles=[0.0, torch.pi / 2, torch.pi / 2, 0.0],
)

layup.plot()

# Reference orientations to sweep
angles = torch.linspace(0.0, 180.0, 19)

Triangular shell elements

Rotate the global reference direction (the 0° ply direction) from 0° to 180° and record the maximum tip deflection for each angle.

# 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 masks for boundaries
left = nodes[:, 0] < 0.01
corner_1 = (nodes[:, 0] > 99.9) & (nodes[:, 1] < 0.01)
corner_2 = (nodes[:, 0] > 99.9) & (nodes[:, 1] > 49.99)

# Sweep the reference orientation and record the maximum tip deflection
deflections_tri = []
for angle in angles:
    # Build orientation vector from angle in degrees
    theta = math.radians(angle.item())
    orientation = torch.tensor([math.cos(theta), math.sin(theta), 0.0])

    # Build cantilever problem and solve
    plate = Shell(nodes, elements, layup, orientation=orientation)
    plate.constraints[left] = True
    plate.forces[corner_1, 2] = -2.0
    plate.forces[corner_2, 2] = -2.0
    u, _, _, _, _ = plate.solve()
    deflections_tri.append(u[:, 2].abs().max().item())

deflections_tri = torch.tensor(deflections_tri)

Quadrilateral shell elements

The same sweep on a quadrilateral mesh carrying the same nodes.

# 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 masks for boundaries
left = nodes[:, 0] < 0.01
corner_1 = (nodes[:, 0] > 99.9) & (nodes[:, 1] < 0.01)
corner_2 = (nodes[:, 0] > 99.9) & (nodes[:, 1] > 49.99)

# Sweep the reference orientation and record the maximum tip deflection
deflections_quad = []
for angle in angles:
    # Build orientation vector from angle in degrees
    theta = math.radians(angle.item())
    orientation = torch.tensor([math.cos(theta), math.sin(theta), 0.0])

    # Build cantilever problem and solve
    plate = Shell(nodes, elements, layup, orientation=orientation)
    plate.constraints[left] = True
    plate.forces[corner_1, 2] = -2.0
    plate.forces[corner_2, 2] = -2.0
    u, _, _, _, _ = plate.solve()
    deflections_quad.append(u[:, 2].abs().max().item())

deflections_quad = torch.tensor(deflections_quad)

Comparison

The deflection is minimal if the 0° outer plies are oriented with the plate. The maximum deflection does not occur at 90°, because a slight off-angle introduces twisting which increases the maximum deflection.

# ABAQUS reference deflections, with triangular and quadrilateral elements
abq_s3_half = [6.60, 8.09, 10.66, 13.85, 17.06, 19.68, 21.16, 21.60, 21.50, 21.13]
abq_s4_half = [6.62, 8.20, 10.92, 14.28, 17.59, 20.12, 21.41, 21.73, 21.57, 21.15]
abq_s3_deflection = abq_s3_half + abq_s3_half[-2::-1]
abq_s4_deflection = abq_s4_half + abq_s4_half[-2::-1]

plt.plot(angles, abq_s3_deflection, ".-", color="gray", label="ABAQUS (S3)")
plt.plot(angles, abq_s4_deflection, ".-", color="black", label="ABAQUS (S4)")
plt.plot(angles, deflections_tri, ".--", color="deeppink", label="torch-fem (Tria1)")
plt.plot(angles, deflections_quad, ".--", color="teal", label="torch-fem (Quad1)")
plt.xlabel("Reference orientation")
plt.ylabel("Maximum tip deflection")
plt.grid(True)
plt.legend()
plt.xticks(range(0, 181, 30), [f"{t}°" for t in range(0, 181, 30)])
plt.show()