Coupling a lattice, a solid and a shell

A cantilever beam of bar elements in an open lattice at the clamped root, solid elements in the middle third, and shell elements to the tip. An Assembly ties the three parts together, and a line load along the free edge of the shell bends the beam.

import torch

from torchfem import Assembly, ReferencePoint, Shell, Solid, Truss
from torchfem.elements import linear_to_quadratic
from torchfem.materials import (
    IsotropicElasticity1D,
    IsotropicElasticity3D,
    IsotropicElasticityPlaneStress,
)
from torchfem.mesh import cube_hexa, mesh_to_lattice, rect_tri

torch.set_default_dtype(torch.float64)

# Overall dimensions
length = 10.0
width = 1.0
thickness = 0.4

# Interface positions
truss_end = length / 3
solid_end = 2 * length / 3

# Material properties
E = 1000.0
nu = 0.3

truss_material = IsotropicElasticity1D(E)
solid_material = IsotropicElasticity3D(E, nu)
shell_material = IsotropicElasticityPlaneStress(E, nu)

Lattice

# Four bays of one hexahedral cell each, braced with a diagonal cross per face
lattice, cells = cube_hexa(5, 2, 2, truss_end, width, thickness)
bars, members = mesh_to_lattice(lattice, cells, "cross")

# Offset nodes
bars += torch.tensor([0.0, 0.0, -thickness / 2])

truss = Truss(bars, members, truss_material)
truss.areas[:] = width * thickness / 12
truss.constraints[bars[:, 0] == 0.0] = True

truss.plot()

Solid

# Quadratic elements to avoid shear locking
nodes, elements = cube_hexa(7, 3, 3, solid_end - truss_end, width, thickness)
nodes, elements = linear_to_quadratic(nodes, elements)

# Offset nodes
nodes += torch.tensor([truss_end, 0.0, -thickness / 2])

# Solid model
solid = Solid(nodes, elements, solid_material)
solid.plot()

Shell

plate, triangles = rect_tri(11, 3, length - solid_end, width / 2)
# Offset nodes
plate += torch.tensor([solid_end, width / 4])
# Add z coordinate
plate = torch.hstack([plate, torch.zeros((len(plate), 1))])

# Shell model
shell = Shell(plate, triangles, shell_material, thickness=thickness)

# Line load along the free edge
tip = plate[:, 0] == length
shell.forces[tip, 2] = -0.05 / tip.sum()
shell.plot(bcs=True)

Assembly

# Truss and solid are both connected via a reference point to ensure that rotations are
# properly transferred. The reference point is located at the center of the flange.
flange = ReferencePoint([truss_end, width / 2, 0.0])

# Create the assembly
assembly = Assembly([truss, solid, shell, flange])

# Couplings
assembly.coupling(truss, bars[:, 0] == truss_end, flange)
assembly.coupling(solid, nodes[:, 0] == truss_end, flange)
assembly.coupling(solid, nodes[:, 0] == solid_end, shell, plate[:, 0] == solid_end)

# Plot assembled structure with boundary conditions
assembly.plot(bcs=True)

Solution

# Solve without aggregation to keep the shell integration points
u, _, sigma, _, _ = assembly.solve(aggregate_integration_points=False)

# [1D in truss /  3x3 in solid /  [n_int, 3x3] in shell / None in reference point]
sigma = [sigma[0], sigma[1].mean(dim=0)[:, 0, 0], sigma[2][-1][:, 0, 0], None]
assembly.plot(u=u, element_property={"Stress": sigma}, bcs=True)