import numpy as np
import openmdao.api as om
from openaerostruct.structures.utils import norm, unit, norm_d, unit_d, cross_d
[docs]
class VonMisesTube(om.ExplicitComponent):
"""Compute the von Mises stress in each element.
parameters
----------
nodes[ny, 3] : numpy array
Flattened array with coordinates for each FEM node.
radius[ny-1] : numpy array
Radii for each FEM element.
disp[ny, 6] : numpy array
Displacements of each FEM node.
Returns
-------
vonmises[ny-1, 2] : numpy array
von Mises stress magnitudes for each FEM element.
"""
def initialize(self):
self.options.declare("surface", types=dict)
def setup(self):
self.surface = surface = self.options["surface"]
ny = self.ny = surface["mesh"].shape[1]
self.add_input("nodes", val=np.zeros((ny, 3)), units="m")
self.add_input("radius", val=np.zeros((ny - 1)), units="m")
self.add_input("disp", val=np.zeros((ny, 6)), units="m")
self.add_output("vonmises", val=np.zeros((ny - 1, 2)), units="N/m**2")
self.E = surface["E"]
self.G = surface["G"]
row = np.concatenate([np.zeros(6), np.ones(6)])
rows = np.tile(row, ny - 1) + np.repeat(2 * np.arange(ny - 1), 12)
col = np.tile(np.arange(6), 2)
cols = np.tile(col, ny - 1) + np.repeat(3 * np.arange(ny - 1), 12)
self.declare_partials("*", "nodes", rows=rows, cols=cols)
rows = np.arange(2 * (ny - 1))
cols = np.repeat(np.arange(ny - 1), 2)
self.declare_partials("*", "radius", rows=rows, cols=cols)
row = np.concatenate([np.zeros(12), np.ones(12)])
rows = np.tile(row, ny - 1) + np.repeat(2 * np.arange(ny - 1), 24)
col = np.tile(np.arange(12), 2)
cols = np.tile(col, ny - 1) + np.repeat(6 * np.arange(ny - 1), 24)
self.declare_partials("*", "disp", rows=rows, cols=cols)
def compute(self, inputs, outputs):
dtype = float
if self.under_complex_step:
dtype = complex
self.T = np.zeros((3, 3), dtype=dtype)
self.x_gl = np.array([1, 0, 0], dtype=dtype)
radius = inputs["radius"]
disp = inputs["disp"]
nodes = inputs["nodes"]
T = self.T
E = self.E
G = self.G
x_gl = self.x_gl
num_elems = self.ny - 1
for ielem in range(num_elems):
P0 = nodes[ielem, :]
P1 = nodes[ielem + 1, :]
L = norm(P1 - P0)
x_loc = unit(P1 - P0)
y_loc = unit(np.cross(x_loc, x_gl))
z_loc = unit(np.cross(x_loc, y_loc))
T[0, :] = x_loc
T[1, :] = y_loc
T[2, :] = z_loc
u0x, u0y, u0z = T.dot(disp[ielem, :3])
r0x, r0y, r0z = T.dot(disp[ielem, 3:])
u1x, u1y, u1z = T.dot(disp[ielem + 1, :3])
r1x, r1y, r1z = T.dot(disp[ielem + 1, 3:])
tmp = np.sqrt((r1y - r0y) ** 2 + (r1z - r0z) ** 2)
sxx0 = E * (u1x - u0x) / L + E * radius[ielem] / L * tmp
sxx1 = E * (u0x - u1x) / L + E * radius[ielem] / L * tmp
sxt = G * radius[ielem] * (r1x - r0x) / L
outputs["vonmises"][ielem, 0] = np.sqrt(sxx0**2 + 3 * sxt**2)
outputs["vonmises"][ielem, 1] = np.sqrt(sxx1**2 + 3 * sxt**2)
def compute_partials(self, inputs, partials):
radius = inputs["radius"]
disp = inputs["disp"]
nodes = inputs["nodes"]
T = self.T
E = self.E
G = self.G
x_gl = self.x_gl
num_elems = self.ny - 1
for ielem in range(num_elems):
# Compute the coordinate delta between the two element end points
P0 = nodes[ielem, :]
P1 = nodes[ielem + 1, :]
dP = P1 - P0
# Compute the derivative of element length
L = norm(dP)
dLddP = norm_d(dP)
# unit function converts a vector to a unit vector
# calculate the transormation to the local element frame.
# We use x_gl to provide a reference axis to reference
x_loc = unit(dP)
dxdP = unit_d(dP)
y_loc = unit(np.cross(x_loc, x_gl))
dtmpdx, _ = cross_d(x_loc, x_gl)
dydtmp = unit_d(np.cross(x_loc, x_gl))
dydP = dydtmp.dot(dtmpdx).dot(dxdP)
z_loc = unit(np.cross(x_loc, y_loc))
dtmpdx, dtmpdy = cross_d(x_loc, y_loc)
dzdtmp = unit_d(np.cross(x_loc, y_loc))
dzdP = dzdtmp.dot(dtmpdx).dot(dxdP) + dzdtmp.dot(dtmpdy).dot(dydP)
T[0, :] = x_loc
T[1, :] = y_loc
T[2, :] = z_loc
u0x = x_loc.dot(disp[ielem, :3])
r0x, r0y, r0z = T.dot(disp[ielem, 3:])
u1x = x_loc.dot(disp[ielem + 1, :3])
r1x, r1y, r1z = T.dot(disp[ielem + 1, 3:])
# #$$$$$$$$$$$$$$$$$$$$$$$$$$
# The derivatives of the above code wrt displacement all boil down to sections of the T matrix
dxddisp = T[0, :]
dyddisp = T[1, :]
dzddisp = T[2, :]
# The derivatives of the above code wrt T all boil down to sections of the #displacement vector
du0dloc = disp[ielem, :3]
dr0dloc = disp[ielem, 3:]
du1dloc = disp[ielem + 1, :3]
dr1dloc = disp[ielem + 1, 3:]
# $$$$$$$$$$$$$$$$$$$$$$$$$$
# Original code
# $$$$$$$$$$$$
tmp = np.sqrt((r1y - r0y) ** 2 + (r1z - r0z) ** 2) + 1e-50 # added eps to avoid 0 disp singularity
sxx0 = E * (u1x - u0x) / L + E * radius[ielem] / L * tmp
sxx1 = E * (u0x - u1x) / L + E * radius[ielem] / L * tmp
sxt = G * radius[ielem] * (r1x - r0x) / L
# $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
dtmpdr0y = 1 / tmp * (r1y - r0y) * -1
dtmpdr1y = 1 / tmp * (r1y - r0y)
dtmpdr0z = 1 / tmp * (r1z - r0z) * -1
dtmpdr1z = 1 / tmp * (r1z - r0z)
# Combine all of the derivtives for tmp
dtmpdDisp = np.zeros(12)
dr0xdDisp = np.zeros(12)
dr1xdDisp = np.zeros(12)
# r0 term
dr0xdDisp[3:6] = dxddisp
dr1xdDisp[9:12] = dxddisp
dtmpdDisp[3:6] = dtmpdr0y * dyddisp
dtmpdDisp[3:6] += dtmpdr0z * dzddisp
dtmpdDisp[9:12] = dtmpdr1y * dyddisp
dtmpdDisp[9:12] += dtmpdr1z * dzddisp
# x_loc, y_loc and z_loc terms
# (dttmpx_loc is zeros, so don't compute with it)
dtmpdy_loc = dtmpdr0y * dr0dloc + dtmpdr1y * dr1dloc
dtmpdz_loc = dtmpdr0z * dr0dloc + dtmpdr1z * dr1dloc
dtmpdP = dtmpdy_loc.dot(dydP) + dtmpdz_loc.dot(dzdP)
dsxx0dtmp = E * radius[ielem] / L
dsxx0du0x = -E / L
dsxx0du1x = E / L
dsxx0dL = -E * (u1x - u0x) / (L * L) - E * radius[ielem] / (L * L) * tmp
dsxx1dtmp = E * radius[ielem] / L
dsxx1du0x = E / L
dsxx1du1x = -E / L
dsxx1dL = -E * (u0x - u1x) / (L * L) - E * radius[ielem] / (L * L) * tmp
dsxx0dP = (
dsxx0dtmp * dtmpdP + dsxx0du0x * du0dloc.dot(dxdP) + dsxx0du1x * du1dloc.dot(dxdP) + dsxx0dL * dLddP
)
dsxx1dP = (
dsxx1dtmp * dtmpdP + dsxx1du0x * du0dloc.dot(dxdP) + dsxx1du1x * du1dloc.dot(dxdP) + dsxx1dL * dLddP
)
# Combine sxx0 and sxx1 terms
# Start with the tmp term
dsxx0dDisp = dsxx0dtmp * dtmpdDisp
dsxx1dDisp = dsxx1dtmp * dtmpdDisp
# Now add the direct u dep
dsxx0dDisp[0:3] = dsxx0du0x * dxddisp
dsxx0dDisp[6:9] = dsxx0du1x * dxddisp
dsxx1dDisp[0:3] = dsxx1du0x * dxddisp
dsxx1dDisp[6:9] = dsxx1du1x * dxddisp
# Combine sxt term
dsxtdr0x = -G * radius[ielem] / L
dsxtdr1x = G * radius[ielem] / L
dsxtdL = -G * radius[ielem] * (r1x - r0x) / (L * L)
dsxtdP = dsxtdr0x * (dr0dloc.dot(dxdP)) + dsxtdr1x * (dr1dloc.dot(dxdP)) + dsxtdL * dLddP
# disp
dsxtdDisp = dsxtdr0x * dr0xdDisp + dsxtdr1x * dr1xdDisp
# radius derivatives
dsxxdrad = E / L * tmp
dsxtdrad = G * (r1x - r0x) / L
fact = 1.0 / (np.sqrt(sxx0**2 + 3 * sxt**2))
dVm0dsxx0 = sxx0 * fact
dVm0dsxt = 3 * sxt * fact
fact = 1.0 / (np.sqrt(sxx1**2 + 3 * sxt**2))
dVm1dsxx1 = sxx1 * fact
dVm1dsxt = 3 * sxt * fact
ii = 2 * ielem
partials["vonmises", "radius"][ii] = dVm0dsxx0 * dsxxdrad + dVm0dsxt * dsxtdrad
partials["vonmises", "radius"][ii + 1] = dVm1dsxx1 * dsxxdrad + dVm1dsxt * dsxtdrad
ii = 24 * ielem
partials["vonmises", "disp"][ii : ii + 12] = dVm0dsxx0 * dsxx0dDisp + dVm0dsxt * dsxtdDisp
partials["vonmises", "disp"][ii + 12 : ii + 24] = dVm1dsxx1 * dsxx1dDisp + dVm1dsxt * dsxtdDisp
# Compute terms for the nodes
ii = 12 * ielem
dVm0_dnode = dVm0dsxx0 * dsxx0dP + dVm0dsxt * dsxtdP
partials["vonmises", "nodes"][ii : ii + 3] = -dVm0_dnode
partials["vonmises", "nodes"][ii + 3 : ii + 6] = dVm0_dnode
dVM1_dnode = dVm1dsxx1 * dsxx1dP + dVm1dsxt * dsxtdP
partials["vonmises", "nodes"][ii + 6 : ii + 9] = -dVM1_dnode
partials["vonmises", "nodes"][ii + 9 : ii + 12] = dVM1_dnode