##calculate the depth of the Beams cross section
## initialization of variables
import math
E=200. ##G Pa
Y=250. ##M Pa
SF=1.9
w=1. ##kN/m
L=3. ##m
S_max=Y
## Calculations
E=E*10**9
Y=Y*10**6
w=w*10**3
Mx=-SF*w*L**2/2.
S_max=S_max*10**6
k=2 ## c_max=h/k
##Formula to be used
## S_max=abs(Mx)*c_max/Ix
## Note that c_max=h/2 and Ix=h^4/24
h=(abs(Mx)*24./(k*S_max))**(1./3.)
print"%s %.3f %s"%('h = ',h,' m')
##calculate the maxmum tensile and maximum normal stress
## initialization of variables
import math
P1=1.5 ##kN
P2=4.5 ##kN
## part (a)
A=1000. ##mm^2
A1=500. ##mm^2
A2=500. ##mm^2
##calculation
A=A*10**-6
A1=A1*10**-6
A2=A2*10**-6
y1=25*10**-3
y2=55*10**-3
c1=(A1*y1+A2*y2)/A
c2=60.*10**-3-c1 ## c1+c2=60 mm
y_1=c1-25.*10**-3
y_2=c2-5*10**-3
b1=50.*10**-3
h1=10.*10**-3
h2=50.*10**-3
b2=10.*10**-3
Ix=1/12.*b1*h1**3 + A1*y_1**2 + 1/12.*b2*h2**3 + A2*y_2**2
print('part (a)')
R1=2550. ##N
Vy=750. ##N
Mx=975. ##N.m
S_zzT=Mx*c1/Ix
S_zzC=Mx*(-c2)/Ix
print"%s %.2f %s"%('\n Maximum Tensile stress = ',S_zzT/10**6,'MPa')
print"%s %.2f %s"%('\n Maximum Compressive stress =',S_zzC/10**6,' MPa')
##calculate the maximum tensile and compressive stress and stress magnitudes
## initialization of variables
import math
P=12. ##kN
Phi=math.pi/3.
## calculations
L=3. ##m
P=12. ##kN
A=10000. ##mm^2
Ix=39.69*10**6 ##mm^4
yo=82. ##mm
Iy=30.73*10**6 ##mm^4
Ixy=0.
P=P*10**3
Ix=Ix*10**-12
Iy=Iy*10**-12
alpha=math.atan(-Ix/(Iy*math.tan(Phi)))
M=-L*P
Mx=M*math.sin(Phi)
yA=-118.*10**-3 ##m
xA=-70.*10**-3 ##m
xB=-xA
yB=82*10**-3 ##m
S_A=Mx*(yA-xA*math.tan(alpha))/(Ix-Ixy*math.tan(alpha))
S_B=Mx*(yB-xB*math.tan(alpha))/(Ix-Ixy*math.tan(alpha))
print"%s %.2f %s"%(' Sigma A = ',S_A/10**6,' M Pa \n')
print"%s %.2f %s"%(' Sigma B = ',S_B/10**6,' M Pa')
##calculate the load stress relations derived for nonsymmetrical bending
## initialization of variables
import math
P=4. ##kN
L=1.2 ##m
A=1900. ##mm**2
Ix=2.783*10**6 ##mm**4
Iy=1.003*10**6 ##mm**4
Ixy=-0.973*10**6 ##mm**4
P=P*10**3
Ix=Ix*10**-12
Iy=Iy*10**-12
Ixy=Ixy*10**-12
A=1900. ##mm**2
xo=19.74 ##mm
yo=39.74 ##mm
Phi=2.*math.pi/3.
Nr=Ixy-Ix/math.tan(Phi)
Dr=Iy-Ixy/math.tan(Phi)
alpha=math.atan(Nr/Dr)
M=L*P
Mx=M*math.sin(Phi)
yA=39.74*10**-3 ##m
xA=-60.26*10**-3 ##m
xB=19.74*10**-3
yB=-80.26*10**-3 ##m
S_A=Mx*(yA-xA*math.tan(alpha))/(Ix-Ixy*math.tan(alpha))
S_B=Mx*(yB-xB*math.tan(alpha))/(Ix-Ixy*math.tan(alpha))
print('part (a)')
print"%s %.2f %s"%('\n Sigma A = ',S_A/10**6,' M Pa \n')
print"%s %.2f %s"%(' Sigma B = ',S_B/10**6,' M Pa')
## part (b)
th=1/2.*math.atan(-2.*Ixy/(Ix-Iy))
th1=0.415 ##rad
th2=-1.156 ##rad
IX=Ix*(math.cos(th1))**2+Iy*(math.sin(th1))**2-2*Ixy*math.sin(th1)*math.cos(th1)
IY=Ix+Iy-IX
Phi=2*math.pi/3-th1
alphA=-IX/(IY*math.tan(Phi))
alpha=alphA+th1
XA=xA*math.cos(th1)+yA*math.sin(th1)
YA=yA*math.cos(th1)-xA*math.sin(th1)
XB=xB*math.cos(th1)+yB*math.sin(th1)
YB=yB*math.cos(th1)-xB*math.sin(th1)
MX=M*math.sin(Phi)
MY=-M*math.cos(Phi)
S_A=MX*YA/IX-MY*XA/IY
S_B=MX*YB/IX-MY*XB/IY
print('\n part (b)')
print"%s %.2f %s"%('\n Sigma A = ',S_A/10**6,' M Pa \n')
print"%s %.2f %s"%(' Sigma B = ',S_B/10**6,' M Pa')
##calculate the maximum load
## initialization of variables
import math
A=3085.9 ##mm**2
Ix=29.94e-6 ##m**4
Iy=4.167e-6 ##m**4
Ixy=0.
ybar=207.64 ##mm
tau_max=165e6 ##Pa
##calculations
A=A*1e-6
ybar=ybar*1e-3
Mxk=-6.1*math.cos(math.pi/6.) ## Mx=Mxk*P
Myk=-6.1*math.sin(math.pi/6.) ##My=Myk*P
## Equation to be followed
## S_zz=Mx*y/Ix-My*x/Iy
## At A x=100 mm y=-92.36 mm
x=100e-3
y=-92.36e-3
S_zzA=Mxk*y/Ix-Myk*x/Iy ##Sigma_zz=S_zz*P
## At B x=-100 mm y=-92.36 mm
x=-100e-3
y=-92.36e-3
S_zzB=Mxk*y/Ix-Myk*x/Iy ##Sigma_zz=S_zz*P
## At C x=-3.125 mm y=207.64 mm
x=-3.125e-3
y=207.64e-3
S_zzC=Mxk*y/Ix-Myk*x/Iy ##Sigma_zz=S_zz*P
## To find P
P=2*tau_max/max(S_zzA,S_zzB,S_zzC)
print"%s %.2f %s"%('P = ',P/10**3,' kN')
##calculate the maximum deflection of the beam
## initialization of variables
import math
P=35. ##kN
Phi=5.*math.pi/9.
E=72e9 ##Pa
L=3. ##m
Ix=39.69*10**6 ##mm**4
Iy=30.73*10**6 ##mm**4
Ixy=0.
##calculations
P=P*1e3
Ix=Ix*10**-12
Iy=Iy*10**-12
alpha=math.atan(-Ix/(Iy*math.tan(Phi)))
M=P*L/4.
Mx=M*math.sin(Phi)
yA=-118.*10**-3 ##m
xA=70.*10**-3 ##m
xB=-xA
yB=82.*10**-3 ##m
S_comp=Mx*(yA-xA*math.tan(alpha))/(Ix-Ixy*math.tan(alpha))
S_tens=Mx*(yB-xB*math.tan(alpha))/(Ix-Ixy*math.tan(alpha))
print"%s %.2f %s"%(' Tensile strength = ',S_tens/10**6,' M Pa \n')
print"%s %.2f %s"%(' Compressive Strength = ',S_comp/10**6,' M Pa')
v=P*L**3*math.sin(Phi)/(48.*E*Ix)
u=-v*math.tan(alpha)
delta=math.sqrt(u**2+v**2)
print"%s %.2f %s"%('\n The total deflection is ',delta*10**3,' mm')
##calculate the magnitude of P necessary ti intiate yeilding in the beam
## initialization of variables
import math
L=3. ##m
Ix=56.43e6 ##mm**4
Iy=18.11e6 ##mm**4
Ixy=22.72e6 ##mm**4
Phi=math.pi/3.
E=200e9 ##Pa
Y=300e6 ##Pa
##calculations
Ix=Ix*10**-12
Iy=Iy*10**-12
Ixy=Ixy*10**-12
yA=-120*10**-3 ##m
xA=-91.*10**-3 ##m
Nr=Ixy-Ix/math.tan(Phi)
Dr=Iy-Ixy/math.tan(Phi)
alpha=math.atan(Nr/Dr)
## M=-L*P To know P we do the following
Mxk=-L*math.sin(Phi)##Mx=Mxk*P
P=Y*(Ix-Ixy*math.tan(alpha))/(Mxk*(yA-xA*math.tan(alpha)))
print"%s %.2f %s"%('P =',P/10**3,' kN \n')
v=P*L**3*math.sin(Phi)/(3*E*(Ix-Ixy*math.tan(alpha)))
u=-v*math.tan(alpha)
delta=math.sqrt(u**2+v**2)
print"%s %.2f %s"%(' deflection = ',delta*10**3,' mm')
## Wrong calculation starting from v in Textbook
##calculate the neutral axis orientation and ratio of the maximum tensile stress
## initialization of variables
import math
Ix=937e+06 ##mm^4
Iy=18.7e+6 ##mm^4
Ixy=0.
yA=305 ##mm
xA=90.5 ##mm
Phi=1.5533 ##rad
##calculations
Ix=Ix*10**-12
Iy=Iy*10**-12
Ixy=Ixy*10**-12
yA=yA*10**-3 ##m
xA=xA*10**-3 ##m
alpha=math.atan(-Ix/(Iy*math.tan(Phi)))
Mxk=math.sin(Phi) ## Mx=Mxk*M
Sigma_Ak1=Mxk*(yA-xA*math.tan(alpha))/(Ix-Ixy*math.tan(alpha))
##Sigma_A=Aigma_Ak*M
## When the plane of the loads coincide with the y axes
Sigma_Ak2=yA/Ix
ratio=Sigma_Ak1/Sigma_Ak2
percent=(ratio-1.)*100.
print"%s %.2f %s"%('alpha = ',alpha,' rad')
print"%s %.2f %s"%('\n The maximum stress in the beam is increased ',percent,' percent when the plane of the loads is merely 1 degre from the symmetrical vertical plane')
## Wrong alpha given in the textbook
##calculate the orientation of the neutral axis and plane of the loads
## initialization of variables
import math
Y=280. ##MPa
AB=40. ##mm
BC=60. ##mm
##calculations
Y=Y*10**6
alpha=math.atan(BC/AB)
C11=20./3. ##mm
C12=-10. ##mm
C21=-20/3. ##mm
C22=10. ##mm
Beta=math.atan((C11-C21)/(C22-C12))
Phi=math.pi/2.+Beta
d=math.sqrt((AB/2.-C11)**2+(BC/2.-C22)**2)
d=d*10**-3 ##m
At=1/2.*AB*BC/2.*10**-6
Mp=At*Y*d
print"%s %.2f %s"%('Mp = ',Mp/10**3,' kN.m')