#59 Optimization of beam cross-section

import os
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import minimize_scalar
os.chdir(r'D:\projects\wordpress\ex59')
os.getcwd()
# The trapezoid is formed by removing the top from a traingle of base
# B=48 mm and height H=60 mm. The task is to find the height y of the trapezoid
# that maximizes the section modulus S=Ibar/c, where Ibar is
# second moment of area about xbar-axis that passes through the centroid C
# of the cross-section. By optimizing the section modulus, we
# minimize the maximum bending stress sigma_max = M/S in the beam, where M is
# the bending moment. 

def f(y):
    B = 48.0
    H = 60.0
    a = B*(H-y)/H #base of rectangle
    b = (B-a)/2.0 #base of triangle
    A = (B + a)*y/2.0 #area
    Qx = (a*y**2)/2.0 + (b*y**2)/3.0 #first moment of area about x-axis
    d = Qx/A #location of centroid
    c = y-d #distance involved in S
    I = (a*y**3)/3.0 + (b*y**3)/6.0 #second moment of area about x-axis
    Ibar = I - A*d**2 #parallel axis theorem
    S = Ibar/c #section modulus
    S = - S # we want maximize S or minimize -S
    return S

y = np.arange(1, 60, 1)
plt.plot(y,-f(y))
bounds= [1.0,60.0]
res = minimize_scalar(f,method='bounded',bounds=bounds)
ystar=res.x
S_trapezoid=-f(res.x)
plt.plot(ystar,S_trapezoid,'r*',markersize=15)
print("Optimal y={}, optimal S={}".format(ystar,S_trapezoid))
S_triangle=-f(60.0)
print("S of triangle = {}".format(S_triangle))
improvement = (S_trapezoid-S_triangle)/S_triangle*100
print("S of triangle = {} %".format(improvement))
plt.savefig("example59.png", dpi=100)
plt.show()
plt.close()