Optimizing a bridge structure#
Installation of the truss package#
For this session, you will need the Python package Truss that can be download here:
Dowload it, extract the content of the archive and put it in your work directory. Once this is completed, execute the cell below to check the install is working.
%load_ext autoreload
%autoreload 2
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
%matplotlib widget
import sys, copy, os
from scipy import optimize
sys.path.append("truss-master")
try:
import truss
print("Truss is correctly installed")
except:
print("Truss is NOT correctly installed !")
A short truss tutorial is available here:
Building the bridge structure#
In this session, we will modelled a bridge structure using truss and optimize it using various criteria. The basic structure is introduced below. It is made of steel bars and loaded with one vertical force on \(G\). The bridge is symmetrical so only the left half is modelled.
E = 210.0e9 # Young Modulus [Pa]
rho = 7800.0 # Density [kg/m**3]
A = 0.9e-2 # Cross section [m**2]
sigmay = 400.0e6 # Yield Stress [Pa]
# Model definition
model = truss.core.Model() # Model definition
# NODES
# Deck
nA = model.add_node((0.0, 0.0), label="A")
nB = model.add_node((1.0, 0.0), label="B")
nC = model.add_node((2.0, 0.0), label="C")
nD = model.add_node((3.0, 0.0), label="D")
nE = model.add_node((4.0, 0.0), label="E")
nF = model.add_node((5.0, 0.0), label="F")
nG = model.add_node((6.0, 0.0), label="G")
# Upper structure
nH = model.add_node((1.0, 2.0), label="H")
nI = model.add_node((2.0, 2.0), label="I")
nJ = model.add_node((3.0, 2.0), label="J")
nK = model.add_node((4.0, 2.0), label="K")
nL = model.add_node((5.0, 2.0), label="L")
nM = model.add_node((6.0, 2.0), label="M")
# BOUNDARY CONDITIONS
nA.block[1] = True
nG.block[0] = True
nM.block[0] = True
# BARS
AC = model.add_bar(nA, nC, modulus=E, density=rho, section=A, yield_stress=sigmay)
CD = model.add_bar(nC, nD, modulus=E, density=rho, section=A, yield_stress=sigmay)
AD = model.add_bar(nA, nD, modulus=E, density=rho, section=A, yield_stress=sigmay)
CE = model.add_bar(nC, nE, modulus=E, density=rho, section=A, yield_stress=sigmay)
DF = model.add_bar(nD, nF, modulus=E, density=rho, section=A, yield_stress=sigmay)
DE = model.add_bar(nD, nE, modulus=E, density=rho, section=A, yield_stress=sigmay)
EF = model.add_bar(nE, nF, modulus=E, density=rho, section=A, yield_stress=sigmay)
EG = model.add_bar(nE, nG, modulus=E, density=rho, section=A, yield_stress=sigmay)
FH = model.add_bar(nF, nH, modulus=E, density=rho, section=A, yield_stress=sigmay)
FG = model.add_bar(nF, nG, modulus=E, density=rho, section=A, yield_stress=sigmay)
GH = model.add_bar(nG, nH, modulus=E, density=rho, section=A, yield_stress=sigmay)
AI = model.add_bar(nA, nI, modulus=E, density=rho, section=A, yield_stress=sigmay)
CI = model.add_bar(nC, nI, modulus=E, density=rho, section=A, yield_stress=sigmay)
IJ = model.add_bar(nI, nJ, modulus=E, density=rho, section=A, yield_stress=sigmay)
EJ = model.add_bar(nE, nJ, modulus=E, density=rho, section=A, yield_stress=sigmay)
CJ = model.add_bar(nC, nJ, modulus=E, density=rho, section=A, yield_stress=sigmay)
GK = model.add_bar(nG, nK, modulus=E, density=rho, section=A, yield_stress=sigmay)
JK = model.add_bar(nJ, nK, modulus=E, density=rho, section=A, yield_stress=sigmay)
EK = model.add_bar(nE, nK, modulus=E, density=rho, section=A, yield_stress=sigmay)
# STRUCTURAL LOADING
nG.force = np.array([0.0, -1.0e6]) / 2.0
nE.force = np.array([0.0, -1.0e6])
nC.force = np.array([0.0, -1.0e6])
# model.solve()
xlim, ylim = model.bbox(deformed=False)
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect("equal")
# ax.axis("off")
model.draw(
ax, deformed=False, field="stress", label=True, force_scale=1.0e-6, forces=True
)
plt.xlim(xlim)
plt.ylim(ylim)
plt.grid()
plt.xlabel("Axe $x$")
plt.ylabel("Axe $y$")
plt.show()
Detailed results at the nodes#
model.data(at="nodes")
Detailed results on the bars#
model.data(at="bars")
Dead (or structural) mass#
m0 = model.mass()
m0 * 1.0e-3 # Mass in tons !
Questions#
Question 1: Verify that the yield stress is not exceeded anywhere, do you think this structure has an optimimum weight ? You can use the state/failure data available on the whole model.
# Example:
model.data(at="bars").state.failure
# ...
Question 2: Modify all the cross sections at the same time in order to minimize weight while keeping acceptable stress level.
model.data(at="bars").state.stress / 1e6
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[1], line 1
----> 1 model.data(at="bars").state.stress / 1e6
NameError: name 'model' is not defined
Cs = np.abs(model.data(at="bars").state.stress) / sigmay
Cs = Cs.values.astype(np.float64)
Cs
np.exp((Cs - 1) * 100)
def cost(X):
nI.coords[1] = X[0]
nJ.coords[1] = X[1]
nK.coords[1] = X[2]
nD.coords[1] = X[3]
nF.coords[1] = X[4]
nH.coords[1] = X[5]
# nI.coords[1] = X[3]
# nJ.coords[1] = X[4]
# nH.coords[1] = X[5]
model.solve()
m = model.mass()
# Cs = np.abs(model.data(at='bars').state.stress) / sigmay
# Cs = Cs.values.astype(np.float64)
# faillure = np.exp((Cs-1)*1000)
# faillure = np.sum(faillure) + 1
# print(faillure)
return m * np.abs(nG.displacement[1])
import scipy
x0 = np.array([-1.0, -1.0, -1.0, 3.0, 3.0, 3.0])
sol = scipy.optimize.minimize(cost, x0)
sol
model.solve()
xlim, ylim = model.bbox(deformed=False)
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect("equal")
# ax.axis("off")
model.draw(
ax, deformed=False, field="stress", label=True, force_scale=1.0e-6, forces=True
)
plt.xlim(xlim)
plt.ylim(ylim)
plt.grid()
plt.xlabel("Axe $x$")
plt.ylabel("Axe $y$")
x0 = np.array([-1.0, -1.0, -1.0, 3.0, 3.0, 3.0])
bounds = [(-2.0, 0.0), (-2.0, 0.0), (-2.0, 0.0), (0.0, 10.0), (0.0, 10.0), (0.0, 10.0)]
sol = scipy.optimize.minimize(cost, x0, bounds=bounds, method="L-BFGS-B")
sol
model.solve()
xlim, ylim = model.bbox(deformed=False)
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect("equal")
# ax.axis("off")
model.draw(
ax, deformed=False, field="stress", label=True, force_scale=1.0e-6, forces=True
)
plt.xlim(xlim)
plt.ylim(ylim)
plt.grid()
plt.xlabel("Axe $x$")
plt.ylabel("Axe $y$")
Question 3: We want to modify the position along the \(\vec y\) axis of the points \(D\), \(F\) and \(H\) in order to minimize the vertical displacement of the node \(G\) times the mass of the structure \(\alpha\):
\[
\alpha = |u_y(G)| m
\]
Where \(u_y(G)\) is the displacement of the node \(G\) along the \(\vec y\) axis and \(m\) the mass of the whole structure.
Do not further modify the sections determined in question 4. Comment the solution.
Question 4: Same question with displacements also along \(\vec x\) of \(C\), \(D\), \(E\) and \(F\). Is it better ?
Question 5: You can now try to perform topological optimization by removing/merging well chosen beams and nodes. In order to make the structure even more efficient.
Question 6: You are now asked to optimize the cross section along with the position of \(C\), \(D\), \(E\) and \(F\) in order to reach the yield stress in each individual beam.
Auto#
fix_points = np.array([[0.0, 0], [9.0, 0.0], [9.0, 3.0], [9.0, 6.0]])
ex_points = np.atleast_2d(np.arange(1, 9.0, 3))
ex_points = np.append(ex_points, np.zeros_like(ex_points), axis=0).T
fix_points = np.append(fix_points, ex_points, axis=0)
fix_points
nb_fix = len(fix_points)
nb_fix
Nb_inerpoints = 5
points = np.append(fix_points, np.random.rand(Nb_inerpoints, 2) * 9, axis=0)
tri = matplotlib.tri.Triangulation(points[:, 0], points[:, 1])
plt.figure()
plt.plot(points[:, 0], points[:, 1], "o")
plt.triplot(tri)
plt.show()
E = 210.0e9 # Young Modulus [Pa]
rho = 7800.0 # Density [kg/m**3]
A = 0.7e-2 # Cross section [m**2]
sigmay = 400.0e6 # Yield Stress [Pa]
# Model definition
model = truss.core.Model() # Model definition
# NODES
nodes = []
# nA = model.add_node((0.0, 0.0), label="A")
# nodes.append(nA)
# nG = model.add_node((9.0, 0.0), label="G")
# nodes.append(nG)
for i, p in enumerate(points):
nodes.append(model.add_node((p[0], p[1]), label="n{:03d}".format(i)))
for i, e in enumerate(tri.edges):
model.add_bar(
nodes[e[0]], nodes[e[1]], modulus=E, density=rho, section=A, yield_stress=sigmay
)
# BOUNDARY CONDITIONS
nodes[0].block[1] = True
nodes[1].block[0] = True
nodes[2].block[0] = True
nodes[3].block[0] = True
# STRUCTURAL LOADING
for i in range(nb_fix):
nodes[i].force = np.array([0.0, -1.0e6]) / 2.0
# model.solve()
xlim, ylim = model.bbox(deformed=False)
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect("equal")
# ax.axis("off")
model.draw(
ax, deformed=False, field="stress", label=True, force_scale=1.0e-6, forces=True
)
plt.xlim(xlim)
plt.ylim(ylim)
plt.grid()
plt.xlabel("Axe $x$")
plt.ylabel("Axe $y$")
model.solve()
np.abs(model.data(at="bars").state.stress).std() / 1e6
def cost(X):
for i, id_nodes in enumerate(range(nb_fix, len(nodes))):
nodes[id_nodes].coords[0] = X[i]
nodes[id_nodes].coords[1] = X[Nb_inerpoints + i]
model.solve()
m = model.mass()
# err = (np.abs(model.data(at='bars').state.stress) - sigmay)/1e7
# faillure = np.exp(err.values.astype(np.float)).sum()
# return m * np.abs(nodes[1].data().disp.uy) * np.abs(model.data(at='bars').state.stress).std()/1e6
return np.abs(model.data(at="bars").state.stress).std() / 1e6
x0 = np.zeros(2 * Nb_inerpoints)
for i, id_nodes in enumerate(range(nb_fix, len(nodes))):
x0[i] = nodes[id_nodes].coords[0]
x0[Nb_inerpoints + i] = nodes[id_nodes].coords[1]
print(x0)
cost(x0)
import scipy
sol = scipy.optimize.minimize(cost, x0)
sol
xlim, ylim = model.bbox(deformed=False)
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.set_aspect("equal")
# ax.axis("off")
model.draw(
ax, deformed=False, field="stress", label=False, force_scale=1.0e-6, forces=True
)
plt.xlim(xlim)
plt.ylim(ylim)
plt.grid()
plt.xlabel("Axe $x$")
plt.ylabel("Axe $y$")