Orientation of orthotropic materials¶
All unit cubes are clamped at one end and subjected to a uniaxial strain.
In [1]:
import torch
from torchfem import Solid
from torchfem.materials import OrthotropicElasticity3D
from torchfem.mesh import cube_hexa
from torchfem.rotations import axis_rotation
# Set default data type to double precision
torch.set_default_dtype(torch.float64)
# Orthotropic elasticity (pine wood)
material = OrthotropicElasticity3D(
6919.0, 271.0, 450.0, 0.388, 0.375, 0.278, 262.0, 354.0, 34.0
)
In [2]:
# Generate cube
nodes, elements = cube_hexa(5, 5, 5)
# Create model
box = Solid(nodes, elements, material)
# Set constraints
s = 0.1
box.constraints[nodes[:, 0] == 0.0, 0] = True
box.constraints[nodes[:, 0] == 1.0, 0] = True
box.constraints[nodes[:, 1] == 0.5, 1] = True
box.constraints[nodes[:, 2] == 0.5, 2] = True
box.displacements[nodes[:, 0] == 1.0, 0] = s
In [3]:
# Solve
u, f, σ, F, α = box.solve()
print(
f"Mean stress in x is {σ[:, 0, 0].mean():.2f}. Should be {(material.E_1 * s):.2f}."
)
Mean stress in x is 691.90. Should be 691.90.
In [4]:
# Rotate material direction
R = axis_rotation(torch.tensor([0.0, 1.0, 0.0]), torch.pi / 2)
box.material.rotate(R)
# Print engineering constants after rotation
print("Engineering constants after 90° rotation about y axis:")
print(f"E_1: {box.material.E_1[0].item():.2f}")
print(f"E_2: {box.material.E_2[0].item():.2f}")
print(f"E_3: {box.material.E_3[0].item():.2f}")
print(f"ν_12: {box.material.nu_12[0].item():.2f}")
print(f"ν_13: {box.material.nu_13[0].item():.2f}")
print(f"ν_23: {box.material.nu_23[0].item():.2f}")
print(f"G_12: {box.material.G_12[0].item():.2f}")
print(f"G_13: {box.material.G_13[0].item():.2f}")
print(f"G_23: {box.material.G_23[0].item():.2f}")
# Solve
u, f, σ, F, α = box.solve()
print(
f"Mean stress in x is {σ[:, 0, 0].mean():.2f}. Should be {(material.E_3 * s):.2f}."
)
Engineering constants after 90° rotation about y axis: E_1: 450.00 E_2: 271.00 E_3: 6919.00 ν_12: 0.46 ν_13: 0.02 ν_23: 0.02 G_12: 34.00 G_13: 354.00 G_23: 262.00 Mean stress in x is 45.00. Should be 45.00.
In [5]:
# Rotate material direction
R = axis_rotation(torch.tensor([0.0, 0.0, 1.0]), torch.pi / 2)
box.material.rotate(R)
# Print engineering constants after rotation
print("Engineering constants after additional 90° rotation about z axis:")
print(f"E_1: {box.material.E_1[0].item():.2f}")
print(f"E_2: {box.material.E_2[0].item():.2f}")
print(f"E_3: {box.material.E_3[0].item():.2f}")
print(f"ν_12: {box.material.nu_12[0].item():.2f}")
print(f"ν_13: {box.material.nu_13[0].item():.2f}")
print(f"ν_23: {box.material.nu_23[0].item():.2f}")
print(f"G_12: {box.material.G_12[0].item():.2f}")
print(f"G_13: {box.material.G_13[0].item():.2f}")
print(f"G_23: {box.material.G_23[0].item():.2f}")
# Solve
u, f, σ, ε, α = box.solve()
print(
f"Mean stress in x is {σ[:, 0, 0].mean():.2f}. Should be {(material.E_2 * s):.2f}."
)
Engineering constants after additional 90° rotation about z axis: E_1: 271.00 E_2: 450.00 E_3: 6919.00 ν_12: 0.28 ν_13: 0.02 ν_23: 0.02 G_12: 34.00 G_13: 262.00 G_23: 354.00 Mean stress in x is 27.10. Should be 27.10.
Rotation invariance of entire model¶
In [6]:
# Generate cube
nodes, elements = cube_hexa(10, 5, 5, 2.0, 1.0, 1.0)
# Create model
box = Solid(nodes, elements, material)
# Set constraints
s = 0.2
box.constraints[nodes[:, 0] == 0.0, 0] = True
box.constraints[nodes[:, 0] == 2.0, 0] = True
box.constraints[nodes[:, 1] == 0.5, 1] = True
box.constraints[nodes[:, 2] == 0.5, 2] = True
box.displacements[nodes[:, 0] == 2.0, 0] = s
# Solve
u, f, σ, F, α = box.solve()
ε = 0.5 * (F.transpose(-1, -2) + F) - torch.eye(3)
# Plot
box.plot(u=u, node_property={"Disp": u[:, 0]})
print(f"Mean stress in x is {σ[:, 0, 0].mean():.2f}.")
print(f"Mean strain in x is {ε[:, 0, 0].mean():.2f}.")
Mean stress in x is 691.90. Mean strain in x is 0.10.
In [7]:
# Rotate the model to different frame of reference
R = axis_rotation(torch.tensor([0.0, 1.0, 0.0]), torch.pi / 2)
nodes = nodes @ R.T
material.rotate(R)
# Create model
box = Solid(nodes, elements, material)
# Set constraints in rotated frame
s = 0.2
box.constraints[nodes[:, 2] <= 0.01, 2] = True
box.constraints[nodes[:, 2] == 2.0, 2] = True
box.constraints[nodes[:, 1] == 0.5, 1] = True
box.constraints[nodes[:, 0] == -0.5, 0] = True
box.displacements[nodes[:, 2] == 2.0, 2] = s
# Solve
u, f, σ, F, α = box.solve()
ε = 0.5 * (F.transpose(-1, -2) + F) - torch.eye(3)
# Rotate back
box.material.rotate(R.T)
box.nodes = nodes @ R
u = u @ R
f = f @ R
σ = torch.einsum("...ij,...ik,...jl->...kl", σ, R, R)
ε = torch.einsum("...ij,...ik,...jl->...kl", ε, R, R)
# Plot
box.plot(u=u, node_property={"Disp": u[:, 0]})
print(f"Mean stress in x is {σ[:, 0, 0].mean():.2f}.")
print(f"Mean strain in x is {ε[:, 0, 0].mean():.2f}.")
Mean stress in x is 691.90. Mean strain in x is 0.10.