For this session, you will need the Python package Truss that can be download here:
https://github.com/lcharleux/truss/archive/master.zip
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.
In [1]:
%load_ext autoreload
%autoreload 2
import numpy as np
import matplotlib.pyplot as plt
%matplotlib nbagg
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:
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.
In [3]:
E = 210.e9 # Young Modulus [Pa]
rho = 7800. # Density [kg/m**3]
A = 5.e-2 # Cross section [m**2]
sigmay = 400.e6 # Yield Stress [Pa]
# Model definition
model = truss.core.Model() # Model definition
# NODES
nA = model.add_node((0.,0.), label = "A")
nC = model.add_node((3.,0.), label = "C")
nD = model.add_node((3.,3.), label = "D")
nE = model.add_node((6.,0.), label = "E")
nF = model.add_node((6.,3.), label = "F")
nG = model.add_node((9.,0.), label = "G")
nH = model.add_node((9.,3.), label = "H")
# BOUNDARY CONDITIONS
nA.block[1] = True
nG.block[0] = True
nH.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)
# STRUCTURAL LOADING
nG.force = np.array([0., -1.e6])
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.e-6, forces = True)
plt.xlim(xlim)
plt.ylim(ylim)
plt.grid()
plt.xlabel("Axe $x$")
plt.ylabel("Axe $y$")
Out[3]:
In [5]:
model.data(at = "nodes")
Out[5]:
In [6]:
model.data(at = "bars")
Out[6]:
In [7]:
m0 = model.mass()
m0 * 1.e-3 # Mass in tons !
Out[7]:
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.
In [10]:
# Example:
model.data(at = "bars").state.failure.values
#...
Out[10]:
Question 2: Modify all the cross sections at the same time in order to minimize weight while keeping acceptable stress level.
In [20]:
new_section = 9e-3
for bar in model.bars:
bar.section = new_section
model.solve()
model.data(at = "bars").state.failure.values
Out[20]:
In [21]:
m1 = model.mass()
m1 *1.e-3 # mass in tons.
Out[21]:
In [22]:
print("Reduced mass by {0} kg".format(m0 - m1))
Question 2: 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.
In [23]:
def alpha(X):
"""
Cost function
"""
nD.coords[1] = X[0]
nF.coords[1] = X[1]
nH.coords[1] = X[2]
model.solve()
return abs(nG.displacement[1]) * model.mass()
X0 = [nD.coords[1], nF.coords[1], nH.coords[1]]
optimize.minimize(alpha, X0, method = "nelder-mead")
Out[23]:
In [24]:
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.e-6, forces = True)
plt.xlim(xlim)
plt.ylim(ylim)
plt.grid()
plt.xlabel("Axe $x$")
plt.ylabel("Axe $y$")
Out[24]:
In [25]:
model.mass()
Out[25]:
Question 3: Same question with displacements also along $\vec x$ of $C$, $D$, $E$ and $F$. Is it better ?
In [26]:
def alpha(X):
"""
Cost function
"""
nD.coords[1] = X[0]
nF.coords[1] = X[1]
nH.coords[1] = X[2]
nC.coords[0] = X[3]
nD.coords[0] = X[4]
nE.coords[0] = X[5]
nF.coords[0] = X[6]
model.solve()
return abs(nG.displacement[1]) * model.mass()
X0 = [nD.coords[1], nF.coords[1], nH.coords[1], nC.coords[0], nD.coords[0], nE.coords[0], nF.coords[0],]
sol = optimize.minimize(alpha, X0, method = "nelder-mead", options = {"maxfev": 1000})
In [27]:
xlim, ylim = model.bbox(deformed = False)
fig = plt.figure()
plt.clf()
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.e-6, forces = True)
plt.xlim(xlim)
plt.ylim(ylim)
plt.grid()
plt.xlabel("Axe $x$")
plt.ylabel("Axe $y$")
Out[27]:
Question 4: 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.
In [28]:
model2 = truss.core.Model() # Model definition
# NODES
nA = model2.add_node((0.,0.), label = "A")
#nC = model2.add_node((3.,0.), label = "C")
#nD = model2.add_node((3.,3.), label = "D")
#nE = model2.add_node((6.,0.), label = "E")
nF = model2.add_node((6.,3.), label = "F")
nG = model2.add_node((9.,0.), label = "G")
nH = model2.add_node((9.,3.), label = "H")
# BOUNDARY CONDITIONS
nA.block[1] = True
nG.block[0] = True
nH.block[0] = True
#G.force = np.array([0., -1.])*1.e4
#C.force[0] = .2
#C.force[1] = -0.
# BARS
AG = model2.add_bar(nA, nG, modulus = E, density = rho, section = A, yield_stress = sigmay)
FG = model2.add_bar(nF, nG, modulus = E, density = rho, section = A, yield_stress = sigmay)
HG = model2.add_bar(nH, nG, modulus = E, density = rho, section = A, yield_stress = sigmay)
AF = model2.add_bar(nA, nF, modulus = E, density = rho, section = A, yield_stress = sigmay)
FH = model2.add_bar(nF, nH, modulus = E, density = rho, section = A, yield_stress = sigmay)
# STRUCTURAL LOADING
nG.force = np.array([0., -1.e6])
def alpha(X):
"""
Cost function
"""
nF.coords[1] = X[0]
nH.coords[1] = X[1]
nF.coords[0] = X[2]
model2.solve()
return abs(nG.displacement[1]) * model2.mass()
X0 = [nF.coords[1], nH.coords[1], nF.coords[0]]
optimize.minimize(alpha, X0, method = "nelder-mead")
xlim, ylim = model2.bbox(deformed = False)
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.set_aspect("equal")
#ax.axis("off")
model2.draw(ax, deformed = False, field = "stress", label = True, force_scale = 1.e-6, forces = True)
plt.xlim(xlim)
plt.ylim(ylim)
plt.grid()
plt.xlabel("Axe $x$")
plt.ylabel("Axe $y$")
model2.data("bars")
Out[28]:
Question 5: 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.
In [31]:
def beta(X):
"""
Cost function
"""
for i in range(len(X)):
model2.bars[i].section = X[i]
model2.solve()
return ((abs(model2.data("bars").state.stress.values) - sigmay)**2).sum()
X0 = [bar.section for bar in model2.bars]
sol = optimize.minimize(beta, X0, method = "nelder-mead", options = {"maxfev": 10000, "maxiter":10000})
print(sol)
In [32]:
xlim, ylim = model2.bbox(deformed = False)
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.set_aspect("equal")
#ax.axis("off")
model2.draw(ax, deformed = False, field = "stress", label = True, force_scale = 1.e-6, forces = True)
plt.xlim(xlim)
plt.ylim(ylim)
plt.grid()
plt.xlabel("Axe $x$")
plt.ylabel("Axe $y$")
model2.data("bars")
Out[32]:
In [33]:
model2.mass()
Out[33]:
As a conclusion, it is observed that the mass reduction is huge in this case.
In [ ]: