Skip to content

Instantly share code, notes, and snippets.

@terjehaukaas

terjehaukaas/ThinWalledCrossSectionAnalysis.py

Last active February 4, 2020 05:43
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save terjehaukaas/de05752d02b5518e9ae22192925648da to your computer and use it in GitHub Desktop.
Save terjehaukaas/de05752d02b5518e9ae22192925648da to your computer and use it in GitHub Desktop.
Download ZIP
Thin-walled Cross-section Analysis
Raw
ThinWalledCrossSectionAnalysis.py
# ------------------------------------------------------------------------
# The following Python code is implemented by Professor Terje Haukaas at
# the University of British Columbia in Vancouver, Canada. It is made
# freely available online at terje.civil.ubc.ca together with notes,
# examples, and additional Python code. Please be cautious when using
# this code; it may contain bugs and comes without warranty of any form.
# ------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
# ------------------------------------------------------------------------
# INPUT (The last one is used by the code)
# ------------------------------------------------------------------------
# ------------------------------------------------------------------------
# S-SHAPE
t = 10.0
b = 200.0
h = 200.0
vertices = np.array([(-b, 0),
(0, 0),
(0, h),
(b, h)])
edges = np.array([(1, 2, t),
(2, 3, t),
(3, 4, t)])
# ------------------------------------------------------------------------
# ONE CELL
b = 500.0
h = 200.0
ttop = 10.0
tbottom = 10.0
tleft = 15.0
tright = 20.0
vertices = np.array([(0.0, 0.0),
(b, 0.0),
(0.0, h),
(b, h)])
edges = np.array([(1, 2, tbottom),
(3, 4, ttop),
(1, 3, tleft),
(2, 4, tright)])
# ------------------------------------------------------------------------
# LAMBDA-SHAPE
t = 10.0
b = 100.0
h = 200.0
vertices = np.array([(0, 0),
(0, h),
(b, h)])
edges = np.array([(1, 2, t),
(2, 3, t)])
# ------------------------------------------------------------------------
# CELL WITH FLANGES
b = 200.0
h = 150.0
a = 80.0
t = 10.0
vertices = np.array([(-a, h),
(0, h),
(b, h),
(b+a, h),
(0, 0),
(b, 0)])
edges = np.array([(1, 2, t),
(2, 3, t),
(3, 4, t),
(5, 2, t),
(6, 3, t),
(5, 6, t)])
# ------------------------------------------------------------------------
# I-BEAM
ttop = 15.0
tbottom = 15.0
tweb = 15.0
btop = 100.0
bbottom = 150.0
h = 250.0
vertices = np.array([(-btop/2, h),
(0.0, h),
(btop/2, h),
(-bbottom/2, 0.0),
(0.0, 0.0),
(bbottom/2, 0.0)])
edges = np.array([(1, 2, ttop),
(2, 3, ttop),
(2, 5, tweb),
(4, 5, tbottom),
(5, 6, tbottom)])
# ------------------------------------------------------------------------
# L-SHAPE
t = 10.0
b = 100.0
h = 200.0
vertices = np.array([(0, 0),
(0, h),
(b, 0)])
edges = np.array([(1, 2, t),
(1, 3, t)])
# ------------------------------------------------------------------------
# CHANNEL PROFILE
t = 20.0
b = 80.0
h = 200.0
vertices = np.array([(b, 0),
(0, 0),
(0, h),
(b, h)])
edges = np.array([(1, 2, t),
(2, 3, t),
(3, 4, t)])
# ------------------------------------------------------------------------
# DEFINE FUNCTION THAT SEARCHES THE OMEGA DIAGRAM
# ------------------------------------------------------------------------
def getOmegaValue(diagram, node):
for i in range(len(diagram)):
index = 0
value = 0.0
if int((diagram[i])[0]) == int(node):
return int(node), (diagram[i])[1]
return index, value
# ------------------------------------------------------------------------
# IDENTIFY CELLS AND FLANGES
# ------------------------------------------------------------------------
# Initialize lists of paths and cells in the cross-section
pathList = [[int(edges[0, 0]), int(edges[0, 1])]]
cellList = []
# Keep walking around in the cross-section, until all paths are covered
stillPathsToWalk = True
while stillPathsToWalk:
foundMatch = False
# Loop over all paths that have been walked
for i in range(len(pathList)):
# Extract path
currentPath = pathList[i].copy()
# Leave cells alone
isCell = False
for k in range(1, len(currentPath)+1):
num = currentPath.count(k)
if num > 1:
isCell = True
if not isCell:
# Initialize before looping over edges
leftCandidate = []
rightCandidate = []
# Loop over all edges and try to append the current path
for j in range(edges.shape[0]):
# Length of the current path
lengthOfCurrentPath = len(currentPath)
# Get the left-most items in the path
leftEnd = currentPath[0]
nextToLeftEnd = currentPath[1]
# Get the right-most items in the path
rightEnd = currentPath[lengthOfCurrentPath - 1]
nextToRightEnd = currentPath[lengthOfCurrentPath - 2]
# Detect link on the left-hand side, from first column of edges
if int(edges[j, 0]) == leftEnd:
# Make sure it is not the same edge
if int(edges[j, 1]) != nextToLeftEnd:
leftCandidate.append(int(edges[j, 1]))
foundMatch = True
# Is it a cell?
if int(edges[j, 1]) in currentPath:
cellCandidate = currentPath.copy()
cellCandidate.insert(0, int(edges[j, 1]))
for k in range(len(cellCandidate)):
if cellCandidate[0] != int(edges[j, 1]):
cellCandidate.pop(0)
cellList.append(cellCandidate)
# Detect link on the left-hand side, from second column of edges
if int(edges[j, 1]) == leftEnd:
# Make sure it is not the same edge
if int(edges[j, 0]) != nextToLeftEnd:
leftCandidate.append(int(edges[j, 0]))
foundMatch = True
# Is it a cell?
if int(edges[j, 0]) in currentPath:
cellCandidate = currentPath.copy()
cellCandidate.insert(0, int(edges[j, 0]))
for k in range(len(cellCandidate)):
if cellCandidate[0] != int(edges[j, 0]):
cellCandidate.pop(0)
cellList.append(cellCandidate)
# Detect link on the right-hand side, from first column of edges
if int(edges[j, 0]) == rightEnd:
# Make sure it is not the same edge
if int(edges[j, 1]) != nextToRightEnd:
rightCandidate.append(int(edges[j, 1]))
foundMatch = True
# Is it a cell?
if int(edges[j, 1]) in currentPath:
cellCandidate = currentPath.copy()
cellCandidate.append(int(edges[j, 1]))
for k in range(len(cellCandidate)):
if cellCandidate[0] != int(edges[j, 1]):
cellCandidate.pop(0)
cellList.append(cellCandidate)
# Detect link on the right-hand side, from second column of edges
if int(edges[j, 1]) == rightEnd:
# Make sure it is not the same edge
if int(edges[j, 0]) != nextToRightEnd:
rightCandidate.append(int(edges[j, 0]))
foundMatch = True
# Is it a cell?
if int(edges[j, 0]) in currentPath:
cellCandidate = currentPath.copy()
cellCandidate.append(int(edges[j, 0]))
for k in range(len(cellCandidate)):
if cellCandidate[0] != int(edges[j, 0]):
cellCandidate.pop(0)
cellList.append(cellCandidate)
# See if we have any candidates for appending/branching the current path
numLeftCandidates = len(leftCandidate)
numRightCandidates = len(rightCandidate)
if numLeftCandidates == 1:
currentPath.insert(0, leftCandidate[0])
elif numLeftCandidates > 1:
for k in range(1, numLeftCandidates):
newPath = currentPath.copy()
newPath.insert(0, leftCandidate[0])
currentPath.insert(0, leftCandidate[0])
# Insert the new path into the path list
pathList.append(newPath)
if numRightCandidates == 1:
currentPath.append(rightCandidate[0])
elif numRightCandidates > 1:
for k in range(1, numRightCandidates):
newPath = currentPath.copy()
newPath.append(rightCandidate[k])
currentPath.append(rightCandidate[0])
# Insert the new path into the path list
pathList.append(newPath)
# Put the current path back into the path list
pathList[i] = currentPath
# Stop if there were no links anywhere
if not foundMatch:
stillPathsToWalk = False
# Identify all flanges, i.e., segments not part of cells
else:
flangeList = []
for i in range(len(pathList)):
currentList = pathList[i]
for j in range(len(pathList[i])-1):
if len(cellList) == 0:
flangeList.append([currentList[j], currentList[j + 1]])
else:
segmentAppearsInCell = False
for k in range(len(cellList)):
if currentList[j] in cellList[k] and currentList[j+1] in cellList[k]:
segmentAppearsInCell = True
if not segmentAppearsInCell:
flangeList.append([currentList[j], currentList[j+1]])
# Create list of duplicate flanges
flangesThatShouldBeRemoved = []
for i in range(len(flangeList)):
list1 = sorted(flangeList[i])
for j in range(i+1, len(flangeList)):
list2 = sorted(flangeList[j])
if list1 == list2:
flangesThatShouldBeRemoved.append(j)
# Remove the duplicate flanges
flangesThatShouldBeRemoved = np.unique(flangesThatShouldBeRemoved)
uniqueFlanges = []
for i in range(len(flangeList)):
if i not in flangesThatShouldBeRemoved:
uniqueFlanges.append(flangeList[i])
# Create list of duplicate cells
cellsThatShouldBeRemoved = []
for i in range(len(cellList)):
list1 = np.unique(cellList[i])
for j in range(i+1, len(cellList)):
list2 = np.unique(cellList[j])
if np.array_equal(list1, list2):
cellsThatShouldBeRemoved.append(j)
# Remove the duplicate cells
cellsThatShouldBeRemoved = np.unique(cellsThatShouldBeRemoved)
uniqueCells = []
for i in range(len(cellList)):
if i not in cellsThatShouldBeRemoved:
uniqueCells.append(cellList[i])
# Determine max-coordinates of the cells
yMaxPerCell = []
zMaxPerCell = []
yMinPerCell = []
zMinPerCell = []
for i in range(len(uniqueCells)):
currentCell = uniqueCells[i]
yValues = []
zValues = []
for j in range(len(currentCell)):
yValues.append(vertices[int(currentCell[j])-1, 0])
zValues.append(vertices[int(currentCell[j])-1, 1])
yMaxPerCell.append(np.max(yValues))
zMaxPerCell.append(np.max(zValues))
yMinPerCell.append(np.min(yValues))
zMinPerCell.append(np.min(zValues))
# Create list of cells that encompass other cells
cellsThatShouldBeRemoved = []
for i in range(len(uniqueCells)):
for j in range(i+1, len(uniqueCells)):
if yMaxPerCell[j] >= yMaxPerCell[i] and yMinPerCell[j] <= yMinPerCell[i] and \
zMaxPerCell[j] >= zMaxPerCell[i] and zMinPerCell[j] <= zMinPerCell[i]:
cellsThatShouldBeRemoved.append(j)
# Remove the duplicate cells
cellsThatShouldBeRemoved = np.unique(cellsThatShouldBeRemoved)
finalCells = []
for i in range(len(uniqueCells)):
if i not in cellsThatShouldBeRemoved:
finalCells.append(uniqueCells[i])
# Output
print('\n'"Paths:", pathList)
print('\n'"Cells:", finalCells)
print('\n'"Flanges:", uniqueFlanges)
# ------------------------------------------------------------------------
# AREA AND CENTROID
# ------------------------------------------------------------------------
A = 0
yA = 0
zA = 0
areas = []
centroids = []
thetas = []
lengths = []
thicknesses = []
for i in range(len(edges)):
# Coordinates of end 1
end1 = edges[i, 0]
y1 = vertices[int(end1)-1, 0]
z1 = vertices[int(end1)-1, 1]
# Coordinates of end 2
end2 = edges[i, 1]
y2 = vertices[int(end2)-1, 0]
z2 = vertices[int(end2)-1, 1]
# Midpoint coordinates
yMidpoint = 0.5 * (y1 + y2)
zMidpoint = 0.5 * (z1 + z2)
centroids.append([yMidpoint, zMidpoint])
# Length of segment
delta_y = y2-y1
delta_z = z2-z1
L = np.sqrt(delta_y**2 + delta_z**2)
lengths.append(L)
# Angle relative to y-axis
if delta_y != 0:
theta = np.arctan(delta_z / delta_y)
else:
theta = 0.5 * np.pi
thetas.append(theta)
# Area of segment
thickness = edges[i, 2]
area = L * thickness
areas.append(area)
thicknesses.append(thickness)
# Accumulate total area
A += area
# Accumulate A * y and A * z
yA += yMidpoint * area
zA += zMidpoint * area
# Centroid location for the cross-section
yCentroid = yA/A
zCentroid = zA/A
# ------------------------------------------------------------------------
# MOMENTS OF INERTIA and PRINCIPAL AXES ORIENTATION
# ------------------------------------------------------------------------
Iy = 0
Iz = 0
Iyz = 0
for i in range(len(edges)):
# Local moments of inertia
Iyo = lengths[i]**3 * thicknesses[i] * np.sin(thetas[i])**2 / 12.0
Izo = lengths[i]**3 * thicknesses[i] * np.cos(thetas[i])**2 / 12.0
Iyzo = lengths[i]**3 * thicknesses[i] * np.cos(thetas[i]) * np.sin(thetas[i]) / 12.0
# Add contributions from the parallel axis theorem (Steiner's sats)
Iy += Iyo + ((centroids[i])[1] - zCentroid)** 2 * areas[i]
Iz += Izo + ((centroids[i])[0] - yCentroid)** 2 * areas[i]
Iyz += Iyzo + ((centroids[i])[0] - yCentroid) * ((centroids[i])[1] - zCentroid) * areas[i]
# Apply caution with the principal axes rotation: what if Iy=Iz
thetaPrincipal = 0.0
if np.abs(Iyz) < 1e-10:
thetaPrincipal = 0.0
elif np.abs(Iy-Iz) < 1e-10:
if np.abs(Iyz) < 1e-10:
thetaPrincipal = 0.0
else:
thetaPrincipal = np.sign(Iyz) * 0.25 * np.pi
else:
thetaPrincipal = 0.5 * np.arctan(2.0*Iyz/(Iy-Iz))
# Let counter-clockwise rotation of the system be positive
thetaPrincipal *= -1
# Transform the moments of inertia
IyRotated = 0.5*(Iy+Iz) + 0.5*np.sqrt((Iz-Iy)**2 + 4.0* Iyz**2)
IzRotated = 0.5*(Iy+Iz) - 0.5*np.sqrt((Iz-Iy)**2 + 4.0* Iyz**2)
# ------------------------------------------------------------------------
# SAINT VENANT CONSTANT J
# ------------------------------------------------------------------------
J = 0
cellArea = 0
if len(finalCells) > 1:
print('\n'"Multi-cell cross-sections not considered in St. Venant torsion calculations")
# Open cross-section
elif len(finalCells) == 0:
for i in range(len(edges)):
# Pick up contributions to J
J += lengths[i] * thicknesses[i]**3 / 3.0
# Single-cell cross-sections
elif len(finalCells) == 1:
# Loop over edges to calculate the area
circleIntegral = 0
for i in range(len(finalCells[0])-1):
# Get the two points that are connected to this edge
vertex1 = (finalCells[0])[i]
vertex2 = (finalCells[0])[i+1]
# Calculate the length of this edge
y1 = vertices[vertex1-1, 0]
z1 = vertices[vertex1-1, 1]
y2 = vertices[vertex2-1, 0]
z2 = vertices[vertex2-1, 1]
delta_y = y2 - y1
delta_z = z2 - z1
length = np.sqrt( delta_y**2 + delta_z**2)
# Get the thickness of this edge
for i in range(edges.shape[0]):
node1 = int(edges[i, 0])
node2 = int(edges[i, 1])
if node1 == vertex1 and node2 == vertex2:
thickness = edges[i, 2]
elif node1 == vertex2 and node2 == vertex1:
thickness = edges[i, 2]
# Add to the denominator integral
circleIntegral += length / thickness
# Add to the cell area, using the cross-product
yCoordVector1 = y1 - yCentroid
zCoordVector1 = z1 - zCentroid
yCoordVector2 = y2 - yCentroid
zCoordVector2 = z2 - zCentroid
cellArea += 0.5 * np.abs(yCoordVector1*zCoordVector2-yCoordVector2*zCoordVector1)
J = 4.0 * cellArea**2 / circleIntegral
Jcell = J
# Check if there are free flanges too
for i in range(len(uniqueFlanges)):
# Get the two points that are connected to this edge
vertex1 = (uniqueFlanges[i])[0]
vertex2 = (uniqueFlanges[i])[1]
# Calculate the length of this edge
y1 = vertices[vertex1-1, 0]
z1 = vertices[vertex1-1, 1]
y2 = vertices[vertex2-1, 0]
z2 = vertices[vertex2-1, 1]
delta_y = y2 - y1
delta_z = z2 - z1
length = np.sqrt(delta_y**2 + delta_z**2)
# Get the thickness of this edge
for i in range(edges.shape[0]):
node1 = int(edges[i, 0])
node2 = int(edges[i, 1])
if node1 == vertex1 and node2 == vertex2:
thickness = edges[i, 2]
elif node1 == vertex2 and node2 == vertex1:
thickness = edges[i, 2]
J += length * thickness**3 / 3.0
# Contributions to J expressed as fractions
cellFractionOfJ = Jcell/J
# ------------------------------------------------------------------------
# ROTATE AXES IF Iyz HAS VALUE
# ------------------------------------------------------------------------
# Shift and rotate axis system to align with centroid and principal axes
if np.abs(thetaPrincipal) > 1e-10:
for i in range(len(vertices)):
y = vertices[i-1, 0] - yCentroid
z = vertices[i-1, 1] - zCentroid
vertices[i - 1, 0] = y * np.cos(thetaPrincipal) + z * np.sin(thetaPrincipal)
vertices[i - 1, 1] = -y * np.sin(thetaPrincipal) + z * np.cos(thetaPrincipal)
storeYcentroid = yCentroid
storeZcentroid = zCentroid
yCentroid = 0.0
zCentroid = 0.0
# ------------------------------------------------------------------------
# WARPING DIAGRAM PLUS ysc, zsc, Cw
# ------------------------------------------------------------------------
Cw = 0
if len(finalCells) > 1:
print('\n'"Multi-cell cross-sections not considered in warping torsion calculations")
# Open cross-section
elif len(finalCells) == 0:
# Walk through the first path
omegaTrial = []
veryFirstEnd = (pathList[0])[0]
omegaTrial.append([veryFirstEnd, 0.0])
for k in range(len(pathList[0]) - 1):
# Coordinates of end 1
end1 = (pathList[0])[k]
y1 = vertices[int(end1) - 1, 0]
z1 = vertices[int(end1) - 1, 1]
# Coordinates of end 2
end2 = (pathList[0])[k + 1]
y2 = vertices[int(end2) - 1, 0]
z2 = vertices[int(end2) - 1, 1]
# Cross-product
cross = y2 * z1 - y1 * z2
# Set omega values
omegaTrial.append([end2, (omegaTrial[k])[1] + cross])
# Walk through the rest of the paths
for i in range(1, len(pathList)):
# Find the first node in this path that already has an omega value
foundOmegaValue = False
for k in range(len(pathList[i])):
# Check if there is an omega value at this node
node = (pathList[i])[k]
omegaIndex, value = getOmegaValue(omegaTrial, node)
if omegaIndex != 0:
foundOmegaValue = True
locationInPath = k
# Walk forward from first-found-node-with-omega-value
for kk in range(locationInPath, len(pathList[i])-1):
# Nothing to do if next node has omega value too
end2 = (pathList[i])[kk+1]
omegaIndex, value = getOmegaValue(omegaTrial, end2)
if omegaIndex == 0:
# Coordinates of end 1
end1 = (pathList[i])[kk]
y1 = vertices[int(end1)-1, 0]
z1 = vertices[int(end1)-1, 1]
# Omega value at end 1
omegaIndex, value = getOmegaValue(omegaTrial, end1)
# Coordinates of end 2
y2 = vertices[int(end2)-1, 0]
z2 = vertices[int(end2)-1, 1]
# Cross-product
cross = y2 * z1 - y1 * z2
# Amend omega diagram
omegaTrial.append([end2, value + cross])
# Walk backward from first-found-node-with-omega-value
for kk in range(len(pathList[i])-1, locationInPath-1, -1):
# Nothing to do if next node has omega value too
end2 = (pathList[i])[kk-1]
omegaIndex, value = getOmegaValue(omegaTrial, end2)
if omegaIndex == 0:
# Coordinates of end 1
end1 = (pathList[i])[kk]
y1 = vertices[int(end1)-1, 0]
z1 = vertices[int(end1)-1, 1]
# Omega value at end 1
omegaIndex, value = getOmegaValue(omegaTrial, end1)
# Coordinates of end 2
y2 = vertices[int(end2)-1, 0]
z2 = vertices[int(end2)-1, 1]
# Cross-product
cross = y2 * z1 - y1 * z2
# Amend omega diagram
omegaTrial.append([end2, value + cross])
if foundOmegaValue == False:
print("Error: Found a path without any omega values")
# Single-cell cross-sections
elif len(finalCells) == 1:
# Figure out if the path around the cell is clockwise
clockwise = False
ySum = 0
zSum = 0
for i in range(len(finalCells[0])-1):
node = (finalCells[0])[i]
ySum += vertices[node-1, 0]
zSum += vertices[node-1, 1]
yCentre = ySum/float(len(finalCells[0])-1)
zCentre = zSum/float(len(finalCells[0])-1)
y1 = vertices[(finalCells[0])[0]-1, 0]
y2 = vertices[(finalCells[0])[1]-1, 0]
z1 = vertices[(finalCells[0])[0]-1, 1]
z2 = vertices[(finalCells[0])[1]-1, 1]
cross = (y1-yCentre) * (z2-zCentre) - (y2-yCentre) * (z1-zCentre)
if cross < 0:
clockwise = True
# Walk around the cell
omegaTrial = []
veryFirstEnd = (finalCells[0])[0]
omegaTrial.append([veryFirstEnd, 0.0])
for i in range(len(finalCells[0])-2):
# Length of this part
end1 = (finalCells[0])[i]
end2 = (finalCells[0])[i+1]
y1 = vertices[end1-1, 0]
z1 = vertices[end1-1, 1]
y2 = vertices[end2-1, 0]
z2 = vertices[end2-1, 1]
delta_y = y2 - y1
delta_z = z2 - z1
length = np.sqrt( delta_y**2 + delta_z**2)
# Thickness of this part
for k in range(edges.shape[0]):
node1 = int(edges[k, 0])
node2 = int(edges[k, 1])
if node1 == end1 and node2 == end2:
thickness = edges[k, 2]
elif node1 == end2 and node2 == end1:
thickness = edges[k, 2]
# Cross-product
cross = y1 * z2 - y2 * z1
# Correction for cell walls
theSign = -1
if clockwise:
theSign = 1
correction = theSign * cellFractionOfJ * J * length / (2.0 * thickness * cellArea)
# Set omega values
omegaTrial.append([end2, (omegaTrial[i])[1] + cross + correction])
# Now address free flanges obeying kinematic compatibility with the cell they are connected to
for i in range(len(uniqueFlanges)):
# Get the two points that are connected to this flange
end1 = (uniqueFlanges[i])[0]
end2 = (uniqueFlanges[i])[1]
y1 = vertices[end1-1, 0]
z1 = vertices[end1-1, 1]
y2 = vertices[end2-1, 0]
z2 = vertices[end2-1, 1]
# Find the cell node that is connected to this flange
for j in range(len(finalCells[0])-1):
node = (finalCells[0])[j]
yNode = vertices[node - 1, 0]
zNode = vertices[node - 1, 1]
if end1 == node:
# Must now find the cell-node next to "node" that is in the same plane as End 1
forwardNode = 0
backNode = 0
if j == 0:
backNode = (finalCells[0])[len(finalCells[0])-2]
forwardNode = (finalCells[0])[j + 1]
else:
forwardNode = (finalCells[0])[j+1]
backNode = (finalCells[0])[j-1]
yBack = vertices[backNode - 1, 0]
zBack = vertices[backNode - 1, 1]
yForward = vertices[forwardNode - 1, 0]
zForward = vertices[forwardNode - 1, 1]
if zBack == z1 and z1 == z2:
omegaIndex, nodeValue = getOmegaValue(omegaTrial, node)
omegaIndex, backValue = getOmegaValue(omegaTrial, backNode)
change = (nodeValue-backValue)/(yNode-yBack)*(y2-y1)
omegaTrial.append([end2, nodeValue + change])
break
elif yBack == y1 and y1 == y2:
omegaIndex, nodeValue = getOmegaValue(omegaTrial, node)
omegaIndex, backValue = getOmegaValue(omegaTrial, backNode)
change = (nodeValue-backValue)/(zNode-zBack)*(z2-z1)
omegaTrial.append([end2, nodeValue + change])
break
elif zForward == z1 and z1 == z2:
omegaIndex, nodeValue = getOmegaValue(omegaTrial, node)
omegaIndex, forwardValue = getOmegaValue(omegaTrial, forwardNode)
change = (nodeValue-forwardValue)/(yNode-yForward)*(y2-y1)
omegaTrial.append([end2, nodeValue + change])
break
elif yForward == y1 and y1 == y2:
omegaIndex, nodeValue = getOmegaValue(omegaTrial, node)
omegaIndex, forwardValue = getOmegaValue(omegaTrial, forwardNode)
change = (nodeValue-forwardValue)/(zNode-zForward)*(z2-z1)
omegaTrial.append([end2, nodeValue + change])
break
elif end2 == node:
# Now find the cell-node next to "node" that is in the same plane as End 2
forwardNode = 0
backNode = 0
if j == 0:
backNode = (finalCells[0])[len(finalCells[0]) - 2]
forwardNode = (finalCells[0])[j + 1]
else:
forwardNode = (finalCells[0])[j + 1]
backNode = (finalCells[0])[j - 1]
yBack = vertices[backNode - 1, 0]
zBack = vertices[backNode - 1, 1]
yForward = vertices[forwardNode - 1, 0]
zForward = vertices[forwardNode - 1, 1]
if zBack == z1 and z1 == z2:
omegaIndex, nodeValue = getOmegaValue(omegaTrial, node)
omegaIndex, backValue = getOmegaValue(omegaTrial, backNode)
change = (nodeValue - backValue) / (yNode - yBack) * (y1 - y2)
omegaTrial.append([end1, nodeValue + change])
break
elif yBack == y1 and y1 == y2:
omegaIndex, nodeValue = getOmegaValue(omegaTrial, node)
omegaIndex, backValue = getOmegaValue(omegaTrial, backNode)
change = (nodeValue - backValue) / (zNode - zBack) * (z1 - z2)
omegaTrial.append([end1, nodeValue + change])
break
elif zForward == z1 and z1 == z2:
omegaIndex, nodeValue = getOmegaValue(omegaTrial, node)
omegaIndex, forwardValue = getOmegaValue(omegaTrial, forwardNode)
change = (nodeValue - forwardValue) / (yNode - yForward) * (y1 - y2)
omegaTrial.append([end1, nodeValue + change])
break
elif yForward == y1 and y1 == y2:
omegaIndex, nodeValue = getOmegaValue(omegaTrial, node)
omegaIndex, forwardValue = getOmegaValue(omegaTrial, forwardNode)
change = (nodeValue - forwardValue) / (zNode - zForward) * (z1 - z2)
omegaTrial.append([end1, nodeValue + change])
break
# Flip the sign of the omega diagram if the walk was counterclockwise
for i in range(len(omegaTrial)):
(omegaTrial[i])[1] = theSign * (omegaTrial[i])[1]
# Calculate integrals needed to find final omega diagram
intOmega = 0
intyOmega = 0
intzOmega = 0
for i in range(len(edges)):
# End 1 of this edge
end1 = edges[i, 0]
y1 = vertices[int(end1) - 1, 0] - yCentroid
z1 = vertices[int(end1) - 1, 1] - zCentroid
# End 2 of this edge
end2 = edges[i, 1]
y2 = vertices[int(end2) - 1, 0] - yCentroid
z2 = vertices[int(end2) - 1, 1] - zCentroid
# Omega values for this edge
index, omegaTrialEnd1 = getOmegaValue(omegaTrial, end1)
index, omegaTrialEnd2 = getOmegaValue(omegaTrial, end2)
# Pick up stored length and thickness
L = lengths[i]
t = thicknesses[i]
# Integrals
intOmega += 0.5*(omegaTrialEnd1+omegaTrialEnd2)*L*t
intyOmega += L*t/6.0*(y1*(2.0*omegaTrialEnd1+omegaTrialEnd2)+y2*(2.0*omegaTrialEnd2+omegaTrialEnd1))
intzOmega += L*t/6.0*(z1*(2.0*omegaTrialEnd1+omegaTrialEnd2)+z2*(2.0*omegaTrialEnd2+omegaTrialEnd1))
# Shear centre coordinates and normalizing constant
C = -intOmega/A
if np.abs(thetaPrincipal) > 1e-10:
ysc = -intzOmega/IyRotated
zsc = intyOmega/IzRotated
else:
ysc = -intzOmega / Iy
zsc = intyOmega / Iz
# Final omega diagram
omegaFinal = []
minWarp = 0
maxWarp = 0
for i in range(len(omegaTrial)):
node = (omegaTrial[i])[0]
yHere = vertices[node-1, 0] - yCentroid
zHere = vertices[node-1, 1] - zCentroid
trialValue = (omegaTrial[i])[1]
finalValue = trialValue + C + ysc * zHere - zsc * yHere
omegaFinal.append([node, finalValue])
if finalValue < minWarp:
minWarp = finalValue
if finalValue > maxWarp:
maxWarp = finalValue
# Cross-section constant for warping torsion
for i in range(len(edges)):
# Nodes of this edge
end1 = edges[i, 0]
end2 = edges[i, 1]
# Omega values for this edge
index, omegaFinalEnd1 = getOmegaValue(omegaFinal, end1)
index, omegaFinalEnd2 = getOmegaValue(omegaFinal, end2)
# Length and thickness of this edge
L = lengths[i]
t = thicknesses[i]
# Accumulate the integral
if np.sign(omegaFinalEnd1) == np.sign(omegaFinalEnd2):
Cw += L*t/3.0*(omegaFinalEnd1**2 + omegaFinalEnd2**2 + omegaFinalEnd1*omegaFinalEnd2)
else:
Cw += L*t/3.0*(omegaFinalEnd1**2 + omegaFinalEnd2**2 - np.abs(omegaFinalEnd1*omegaFinalEnd2))
# ------------------------------------------------------------------------
# ROTATE AXES BACK
# ------------------------------------------------------------------------
if np.abs(thetaPrincipal) > 1e-10:
yCentroid = storeYcentroid
zCentroid = storeZcentroid
y = ysc
z = zsc
ysc = y * np.cos(thetaPrincipal) - z * np.sin(thetaPrincipal) + yCentroid
zsc = y * np.sin(thetaPrincipal) + z * np.cos(thetaPrincipal) + zCentroid
for i in range(len(vertices)):
y = vertices[i-1, 0]
z = vertices[i-1, 1]
vertices[i - 1, 0] = y * np.cos(thetaPrincipal) - z * np.sin(thetaPrincipal) + yCentroid
vertices[i - 1, 1] = y * np.sin(thetaPrincipal) + z * np.cos(thetaPrincipal) + zCentroid
# ------------------------------------------------------------------------
# PRINT RESULTS
# ------------------------------------------------------------------------
print('\n'"Area:", A)
print('\n'"Centroid y-coordinate: %.2f" % yCentroid)
print('\n'"Centroid z-coordinate: %.2f" % zCentroid)
print('\n'"Moment of inertia Iy: %.2f" % Iy)
print('\n'"Moment of inertia Iz: %.2f" % Iz)
if np.abs(thetaPrincipal) > 1e-10:
print('\n'"Principal axis rotation (degrees): %.2f" % (thetaPrincipal/np.pi*180.0))
print('\n'"Product of inertia Iyz: %.2f" % Iyz)
print('\n'"New principal moment of inertia Iy: %.2f" % IyRotated)
print('\n'"New principal moment of inertia Iz: %.2f" % IzRotated)
print('\n'"Shear centre y-coordinate in original coordinate system: %.2f" % ysc)
print('\n'"Shear centre z-coordinate in original coordinate system: %.2f" % zsc)
else:
print('\n'"Shear centre y-coordinate: %.2f" % ysc)
print('\n'"Shear centre z-coordinate: %.2f" % zsc)
print('\n'"St. Venant torsion constant J: %.2f" % J)
print('\n'"Warping torsion constant Cw: %.2f" % Cw)
# ------------------------------------------------------------------------
# PLOT
# ------------------------------------------------------------------------
for i in range(len(edges)):
# End 1 of this edge
end1 = edges[i, 0]
y1 = vertices[int(end1)-1, 0]
z1 = vertices[int(end1)-1, 1]
# End 2 of this edge
end2 = edges[i, 1]
y2 = vertices[int(end2)-1, 0]
z2 = vertices[int(end2)-1, 1]
# Omega values for this edge
if Cw != 0:
index, omegaFinalEnd1 = getOmegaValue(omegaFinal, end1)
index, omegaFinalEnd2 = getOmegaValue(omegaFinal, end2)
# Discretize and plot the line
segments = 10
for j in range(segments):
y1Here = y1 + j/float(segments)*(y2-y1)
z1Here = z1 + j/float(segments)*(z2-z1)
y2Here = y1 + (j+1)/float(segments)*(y2-y1)
z2Here = z1 + (j+1)/float(segments)*(z2-z1)
if Cw > 1e-5:
warpValueHere = omegaFinalEnd1 + j/float(segments) * (omegaFinalEnd2-omegaFinalEnd1)
colorCode = (warpValueHere - minWarp) / (maxWarp - minWarp)
plt.plot(np.array([y1Here, y2Here]), np.array([z1Here, z2Here]), color=(colorCode, 0, 1 - colorCode), linewidth=6.0, zorder=1)
else:
colorCode = 0.0
plt.plot(np.array([y1Here, y2Here]), np.array([z1Here, z2Here]), color=(0.52, 0.52, 0.52), linewidth=6.0, zorder=1)
# Create the plot
max = np.amax(np.amax((vertices)))
min = np.amin(np.amin(vertices))
yMin = min - (max-min)*0.4
border = (max - min) * 0.1
plt.plot(np.array([yCentroid, max+border]), np.array([zCentroid, zCentroid+(np.abs(max)+border-yCentroid)*np.tan(thetaPrincipal)]), 'k-',linewidth=1.0)
plt.plot(np.array([yCentroid, yMin-border]), np.array([zCentroid, zCentroid-(np.abs(yMin)+border+yCentroid)*np.tan(thetaPrincipal)]), 'k-',linewidth=1.0)
plt.plot(np.array([yCentroid, yCentroid-(np.abs(max)+border-zCentroid)*np.tan(thetaPrincipal)]), np.array([zCentroid, max+border]), 'k-',linewidth=1.0)
plt.plot(np.array([yCentroid, yCentroid+(np.abs(min)+border+zCentroid)*np.tan(thetaPrincipal)]), np.array([zCentroid, min-border]), 'k-',linewidth=1.0)
plt.scatter([yCentroid], [zCentroid], s=100, c='k', marker='x', zorder=2)
plt.scatter([ysc], [zsc], s=100, c='k', marker='x', zorder=2)
plt.axis([yMin-border, max+border, min-border, max+border])
plt.title("Centroid + Shear Centre + Warping")
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment