##calculate the Values of P and Q and small displacement component
## initialization of variables
import math
## part (b)
K1=2. ##N/mm (K1=E1A1/L1)
K2=3. ##N/mm (K2=E2A2/L2)
b1=400. ## mm (b1=h)
h=400. ## mm
b2=300. ##mm
u=30. ##mm
v=40. ##mm
## calculations
## Units convertion
K1=K1*10**3
K2=K2*10**3
b1=b1*10**-3
b2=b2*10**-3
h=h*10**-3
u=u*10**-3
v=v*10**-3
L1=math.sqrt(b1**2+h**2)
L2=math.sqrt(b2**2+h**2)
N1=math.sqrt((b1+u)**2+(h+v)**2)-L1
N2=math.sqrt((b1+u)**2+(h+v)**2)
N3=math.sqrt((b2-u)**2+(h+v)**2)-L2
N4=math.sqrt((b2-u)**2+(h+v)**2)
P=K1*(b1+u)*N1/N2-K2*(b2-u)*N3/N4
Q=K1*(h+v)*N1/N2+K2*(h+v)*N3/N4
## results
print('part (b)')
print"%s %.2f %s"%('\n P = ',P,' N')
print"%s %.2f %s"%('\n Q =',Q,' N')
##calculate the small displacement components
## initialization of variables
import math
## Loads and stresses and dimensions
P=10. ##kN
Q=30. ##kN
S0=70. ##MPa
eps0=0.001
b1=400. ##mm
h=400. ##mm
b2=300. ##mm
A1=300. ##mm^2
A2=300. ##mm^2
## calculations
## Units convertion
P=P*10**3
Q=Q*10**3
S0=S0*10**6
b1=b1*10**-3
b2=b2*10**-3
h=h*10**-3
A1=A1*10**-6
A2=A2*10**-6
L1=math.sqrt(b1**2+h**2)
L2=math.sqrt(b2**2+h**2)
a=L1*(Q*b2+P*h)/(A1*S0*h*(b1+b2))
b=L2*(Q*b1-P*h)/(A2*S0*h*(b1+b2))
c=L1**2*eps0/(b1+b2)
d=L2**2*eps0/(b1+b2)
u=c*math.sinh(a)-d*math.sinh(b)
v=b2/h*c*math.sinh(a)+b1/h*d*math.sinh(b)
## results
print"%s %.2f %s"%('u = ',u*10**3,' mm')
print"%s %.2f %s"%('\n v =',v*10**3,' mm')
##calculate the deflection of the free end of the beam
## initialization of variables
import math
## Material constants
E=200. ##GPa
G=77.5 ## GPa
Lh=5. ## Lh = L/h
## part (b)
rhs1=1.8*Lh*E/G
rhs2=7.*12.*Lh**3./16.
LHS=1.8*Lh*E/G+7.*12.*Lh**3./16.
e=rhs1/LHS*100.
print"%s %.2f %s"%('The error in neglecting small terms is ',e,' per cent')
##calculate the rotation of section B
## initialization of variables
import math
## Specifications
T=2. ##kN.m
E=72. ## G Pa
G=27. ## GPa
b=30. ##mm
h=40. ##mm
d=60. ##mm
l1=400. ##mm
l2=800. ##mm
## calculations
E=E*10**9
G=G*10**9
b=b*10**-3
h=h*10**-3
d=d*10**-3
l1=l1*10**-3
l2=l2*10**-3
T=T*10**3 ##N.m
Ix=b*h**3./12.
J=math.pi*d**4/32.
thB= 2.*l1**3./3.*0.001**2*T/(E*Ix)+T*l2/(G*J)
print"%s %.2f %s"%('The rotation of shaft B is th = ',thB,' rad')
## Wrong answer to an extent in the textbook
##calculate the deflection of the free end of curved beam in he direction and error deflection
## initialization of variables
import math
## specification
R=65. ##mm
E=200. ##GPa
G=77.5 ##GPa
v=0.29
P=6.##kN
##calculations
R=R*10**-3
E=E*10**9
G=G*10**9
P=P*10**3
A=30**2*10**-6
I=30**4./12.*10**-12
q_p1=3*math.pi*P*R/(4.*E*A)+1.2*3.*math.pi*P*R/(4.*G*A)+(9.*math.pi/4.+2.)*P*R**3./(E*I)
print('part (a)')
print"%s %.2f %s"%('\n qp = ',q_p1*10**3,' mm')
##part (b)
## if Un and Us are neglected
q_p2=(9.*math.pi/4.+2.)*P*R**3./(E*I)
e=(q_p1-q_p2)/q_p1*100.
print('\n part (b)')
print"%s %.2f %s"%('\n error = ',e,' per cent')
##calculate the effect of neglecting the strain energy Us due to shear and vertical deflection
## initialization of variables
import math
## part (b)
## Specifications
P=150. ##N
R=200. ##mm
d=20. ##mm
E=200. ##GPa
G=77.5 ##GPa
##calculations
R=R*10**-3.
d=d*10**-3.
E=E*10**9.
G=G*10**9.
r1=R+d/2.
r2=R-d/2.
A=314.*10**-6
I=7850.*10**-12 ##m**4
Ax=3.*math.pi/4.*P*R/(E*A)
Sh=3.*math.pi/4.*1.33*P*R/(G*A)
M=(7.*math.pi/4.+1.)*P*R**3./(E*I)
##qc=3*%pi/4*P*R/(E*A)+3*%pi/4*1.33*P*R/(G*A)+(7*%pi/4+1)*P*R**3/(E*I)
qc=Ax+Sh+M
print"%s %.2f %s %.2f %s %.2f %s %.2f %s"%('qc = ',qc*10**3,' mm 'and 'among which due to Axial is ',Ax*10**3,' mm'and '',Sh*10**3, 'mm'and ' is due to shear, and ',M*10**3,' mm is due to moment')
print('\n which means The concentrations of axial loads and shear are negligible')
##calculate the rotation of member BE caused by the loads P and Q
## initialization of variables
import math
## Material properties and dimensions
E=72. ##G Pa
P=10. ##kN
Q=5. ##kN
Aab=150. ##mm**2
Abc=900. ##mm**2
Acd=900. ##mm**2
Ade=900. ##mm**2
Abd=150. ##mm**2
Abe=150. ##mm**2
Lab=2. ##m
Lbc=2.5 ##m
Lbd=1.5 ##m
Lbe=2.5 ##m
Lcd=2. ##m
Lde=2. ##m
##calculations
E=E*10**9.
P=P*10**3.
Q=Q*10**3.
Aab=150.
Aab=Aab*10**-6.
Abc=Abc*10**-6.
Acd=Acd*10**-6.
Ade=Ade*10**-6.
Abd=Abd*10**-6.
Abe=Abe*10**-6.
M=0.
Nab=4./3.*(Q+2.*P)-5.*M/(3.*Lbe)
dNab=-5./(3.*Lbe)
Nbc=-5./3.*(Q+P)
dNbc=0.
Nbd=Q
dNbd=0.
Nbe=5.*P/3.-4./3.*M/Lbe
dNbe=-4./(3.*Lbe)
Ncd=-4.*P/3.+5./3.*M/Lbe
dNcd=5./(3.*Lbe)
Nde=Ncd
thBE=Nab*Lab*dNab/(E*Aab)+Nbc*Lbc*dNbc/(E*Abc)+Nbd*Lbd*dNbd/(E*Abd)+Nbe*Lbe*dNbe/(E*Abe)+2.*Ncd*Lcd*dNcd/(E*Lcd)
print"%s %.2f %s"%('The rotation of member BE is ',thBE,' rad')
## Wrong answer in the text