Chapter 4 : Pure Bending

Example 4.01, Page number 232

In [4]:
import math

#Variable declaration
c=1.25                                                               # Radius(in)
Sy=36                                                                # Stress(ksi)
b=0.8                                                                # Breadth(in)
h=2.5                                                                # Height(in)

#Calculation         
I=(1/12.0)*(b)*(h)**3                                                  # Centroidal moment of inertia(in**4)
M=(I/c)*(Sy)                                                         # Bending moment(kip.in)

# Result
print ('Bending moment = %lf kip.in' %M)
Bending moment = 30.000000 kip.in

Example 4.02, Page number 233

In [5]:
import math

#Variable declaration
r=12                                                                     # Radius(mm)
p=2.5                                                                    # Mean radius(m)  
E=70                                                                     # Modulus of rigidity(GPa)
n=-1

#Calculation         
Y=(4*r)/(3*(math.pi))                                                    # Ordinate(mm)
c=r-Y                                                                    # Distance from the neutral axis to the point of crossection(mm) 
Em=(c*(10**-3))/p                                                        # Maximum absolute value of the strain
Sm=((E*(pow(10,9)))*Em)/(pow(10,6)*(1.0))                                      # Maximum tensile stress(MPa)
Scomp=(n)*(Y/c)*(Sm)                                                     # Maximum compressive stress(MPa) 

# Result
print ('Maximum tensile stress = %lf MPa' %Sm)
print ('Maximum compressive stress = %lf MPa' %Scomp)
Maximum tensile stress = 193.397171 MPa
Maximum compressive stress = -142.602829 MPa

SAMPLE PROBLEM 4.1, Page number 235

In [2]:
import math

#Variable declaration
sY=40                                                                   # Stress(ksi)
sU=60                                                                   # Stress(ksi)
E=(10.6)*(pow(10,6))                                                    # Modulus of rigidity(psi)
FS=3                                                                    # Factor of safety

#Calculation
#Moment of Inertia
E=(10.6)*(pow(10,6))                                                    # Modulus of rigidity(psi)
I=round(((1/12.0)*3.25*pow(5,3))-((1/12)*(2.75)*pow(4.5,3)),2)            # Centroidal moment of inertia of a rectangle
#Allowable Stress
sALL=(sU/FS)                                                            # Allowable stress(ksi)
#Case(a) Bending Moment
c=(1/2.0)*(5)                                                             # Radius(in)
M=((12.97)*(20))/2.5                                                    # Bending moment(kip.in)
#Case(b) Radius of Curvature
p=round((10.6*pow(10,6)*12.97)/(103.8*pow(10,3)),1)                     # Radius of curvature(in)
p=round((p*0.08333),1)                                                  # Converting into feet(ft)  
#Alternative Solution.
Em=(sALL/(E*(pow(10,-3))*(1.0)))                                              # Maximum strain(in./in)
p=(c/Em)                                                                # Radius of curvature(in)
p=round((p*0.08333),1)                                                  # Converting into feet(ft) 

# Result
print ('Bending moment M for which factor of safety is 3 = %lf kip.in' %M)
print ('Radius of curvature of tube = %lf ft' %p)
Bending moment M for which factor of safety is 3 = 103.760000 kip.in
Radius of curvature of tube = 110.400000 ft

SAMPLE PROBLEM 4.2, Page number 236

In [5]:
import math

#Variable declaration
n=-1

#Calculation
#Centroid
sumA=3000                                                                       # Summing up the area(mm**2) 
M=3                                                                             # Couple(kN.m) 
cA=0.022                                                                        # Distance(m)                 
Y=(114*pow(10,6))/(3000.0)                                                        # Distance(mm)
#Centroidal Moment of Inertia
Ix=((1/12.0)*(90)*(pow(20,3)) + (90*20*pow(12,2)) + ((1/12.0)*(30)*(pow(40,3))) + (30*40*pow(18,2)))/(pow(10,12)*(1.0)) # Centroidal moment of inertia(m**4)  
#Case(a) Maximum Tensile Stress
sA=((M*cA)/(Ix)*(1.0))/(1000.0)                                                          # Maximum tensile stress(MPa) 
#Maximum Compressive Stress
sB=n*(3*0.038)/((868*pow(10,-9)*pow(10,3)))                                      # Maximum compressive stress(MPa) 
#Case(b) Radius of Curvature
p=((165*868*(pow(10,-9)))/(3))*(pow(10,6))                                       # Radius of curvature(m)


# Result
print ('Maximum tensile stress = %lf MPa' %sA)
print ('Maximum compressive stress = %lf MPa' %sB)
print ('Radius of curvature = %lf ft' %p)
Maximum tensile stress = 76.036866 MPa
Maximum compressive stress = -131.336406 MPa
Radius of curvature = 47.740000 ft

Example 4.03, Page number 244

In [16]:
import math

#Variable declaration
Es=29*(pow(10,6))                                                          # Modulus of rigidity(psi)
Eb=15*(pow(10,6))*(1.0)                                                          # Modulus of rigidity(psi) 
M=40                                                                       # Bending moment(kip.in)
h=3                                                                        # Height(3)
b=2.25                                                                     # Breadth(in) 
c=1.5                                                                      # Distance(in)

#Calculation         
n=Es/Eb                                                                    # Ratio
W=0.75*n                                                                   # width(in)
I=(1/12.0)*(b)*((h)**3)                                                      # Moment of inertia of the transformed section(in**4)
Sm=(M*c)/(I)                                                               # Maximum stress in the transformed section(ksi)
Sbrass=Sm                                                                  # Maximum stress in brass portion(ksi)
Ssteel=1.933*(Sbrass)                                                      # Maximum stress in steel portion(ksi)

# Result
print ('Maximum stress in brass portion = %lf ksi' %Sbrass)
print ('Maximum stress in steel portion = %lf ksi' %Ssteel)
Maximum stress in brass portion = 11.851852 ksi
Maximum stress in steel portion = 22.909630 ksi

Example 4.04, Page number 247

In [20]:
import math

#Variable declaration
depth=10                                                                 # Depth(mm)
width=60                                                                 # Width(mm)
thickness=9                                                              # Thickness(mm)
Smax=150                                                                 # Maximum stress(MPa) 
M=180                                                                    # Bending moment(N.m) 

#Calculation         
d=width-(2*depth)                                                        # Distance(mm) 
c=(1/2.0)*d                                                                # Distance(mm) 
b=9                                                                      # Distance(mm)
I=(1/12.0)*(b*(pow(10,-3)))*((d*(pow(10,3)))**3)                           # Moment of inertia of the critical cross section(m**4)
Ratio=((M)*(c)*(pow(10,3)))/(I)                                          # Stress(MPa)
k=150/75.0                                                                 # Factor 
Ratio2=width/(d*1.0)                                                           # Ratio
r=0.13*40                                                                # Radius(mm)
wid=2*r                                                                  # Width(mm)

# Result
print ('Smallest allowable width of the groves = %lf mm' %wid)
Smallest allowable width of the groves = 10.400000 mm

SAMPLE PROBLEM 4.3, Page number 248

In [21]:
import math

# Variable declaration
Es=200                                                                    # Moduluss of rigidity(GPa)
Ew=12.5                                                                   # Moduluss of rigidity(GPa) 

#Transformed Section.
n=(Es/Ew)                                                                 # Ratio
#Neutral Axis
Y=round(((0.160)*(3.2*0.020))/(3.2*0.020+0.470*0.300),2)                  # Distance(m) 
#Centroidal Moment of Inertia
I=round(((1/12)*0.470*(pow(0.3,3)))+(0.470*0.3*(pow(0.05,2)))+((1/12)*(3.2)*(pow(0.020,3)))+(3.2*0.020*(pow(0.160-0.050,2))),5) # Centroidal Moment of Inertia 
#Maximum Stress in Wood
sW=((50*(pow(10,3)))*(0.200))/(2.19*pow(10,-3))                           # Maximum stress in wood(MPa)
sW=round((sW/(pow(10,6))),2)                                              # Rounding
#Stress in Steel
sS=((16)*(50*(pow(10,3)))*(0.120))/(2.19*(pow(10,-3)))                    # Stress in steel(MPa) 
sS=round((sS/(pow(10,6))),1)                                              # Rounding

# Result
print ('Maximum stress in the wood = %lf MPa' %sW)
print ('Stress in steel = %lf MPa' %sS)
Maximum stress in the wood = 4.570000 MPa
Stress in steel = 43.800000 MPa

SAMPLE PROBLEM 4.4, Page number 249

In [23]:
import math
from sympy import symbols, solve

# Variable declaration
x=symbols('x')                                                            # Variable declartion                                                       
d=(5/8.0)                                                                   # Diameter(in)
Es=(29*(pow(10,6)))                                                       # Modulus of elasticity for concrete(psi)
Ec=(3.6*(pow(10,6)))                                                      # Modulus of elasticity for stell(psi)

# Calculation
#Transformed Section
As=round((2)*(math.pi)*(pow(d/2,2)),3)                                    # Cross sectional area(in**2)
n=round(Es/Ec,2)                                                          # Ratio
TA=round(n*As,2)                                                          # Transformed steel area(in**2)
#Neutral Axis
x=(solve((12*x)*(x/2.0)-(4.95)*(4-x)))
#Moment of inertia
I=(1/3.0)*(12)*((1.450)**3) + 4.95*(4-1.450)**2                              # Centroidal moment of inertia(in**4) 
#Maximum Stress in Concrete
sC=round(((40)*(1.450))/(44.4),3)                                         # Maximum stress in concrete(ksi) 
#Stress in Steel
sS=round(((8.06)*(40)*(2.55))/(44.4),2)                                   # Stress in steel(ksi)

# Result
print ('Maximum stress in concrete = %lf MPa' %sC)
print ('Stress in steel = %lf MPa' %sS)
Maximum stress in concrete = 1.306000 MPa
Stress in steel = 18.520000 MPa

Example 4.05, Page number 260

In [24]:
import math
from sympy import symbols, solve 

#Variable declaration
M=36.8
c=60*(pow(10,-3))                                                             # Half length of rod(mm)
b=50*(pow(10,-3))                                                             # Breadth of rod(mm)
Sy=240                                                                        # Stress(MPa) 
yY=symbols('yY')                                                              # Variable declaration
E=200                                                                         # Modulus of elasticity(GPa)
n=-1

#Calculation         
# Case(a)
Rt=(2/3.0)*(b)*(c)**2                                                           # Ratio(m**3)
My=Rt*Sy                                                                      # Maximum elastic moment(kN.m) 
yY=solve(36.8-(3/2.0)*(28.8)*(1-((1/3.0)*(((yY)**2)/(60**2)))),yY)                # 
# Case(b)
Ey=(Sy*(pow(10,6)))/(E*(pow(10,9)*(1.0)))
p=(yY[0]*(pow(10,-3)))/(Ey)

# Result
print ('Thickness of elastic core = %lf mm' %(n*(2*yY[0])))
print ('Radius of curvature = %lf m' %(n*p))
Thickness of elastic core = 80.000000 mm
Radius of curvature = 33.333333 m

Example 4.06, Page number 262

In [27]:
import math

#Variable declaration
M=36.8                                                                           # Bending moment(kN)                              
Sy=240                                                                           # Yield strength(MPa)
yY=40                                                                            # Thickness of elastic core(mm) 
n=-1
Sx=n*35.5*(pow(10,6))                                                            # Stress(Pa)    
E=200*(pow(10,9))

#Calculation         
# Case(a)
Sml=((36.8)/(120*(pow(10,-6))))/(1000)                                           # Residual stress(MPa) 
# Case(b)
Ex=Sx/E                                                                          # Residual strain 
p=(n*(40*(pow(10,-3))))/(Ex)                                                     # Radius of Curvature after Unloading(m)

# Result
print ('Residual stress = %lf MPa' %Sml)
print ('Radius of curvature after unloading = %.lf m' %p)
Residual stress = 306.666667 MPa
Radius of curvature after unloading = 225 m

SAMPLE PROBLEM 4.5, Page number 263

In [28]:
import math

# variable declaration
E=(29*pow(10,6))                                                               # Modulus of elastoplasticity(psi)
sY=50                                                                          # Stress(ksi) 

# Calculation
#Case(a) Onset Of Yield
I=round((1/12.0)*(12)*(pow(16,3))-(1/12.0)*(12-0.75)*(pow(14,3)),0)                # Centroidal moment of inertia(in**4) 
#Bending Moment
sMAX=sY                                                                        # Stress(ksi)
c=8.0                                                                            # Distance(in)
My=(sY*I)/c                                                                    # Bending moment(kip.in)
#Radius of Curvature 
Ey=sY/(E*(1.0))                                                                        # Strain
pY=(c/Ey)/(1000.0)                                                               # Radius of curvature(in)
#Case(b) Flanges Fully Plastic
R1=50*12*1                                                                     # Compressive forces on top(kips) 
R4=R1                                                                          # Compressive forces on top(kips) 
R2=round((1/2.0)*(50)*(7)*(0.75)+0.05,1)                                         # Compressive forces on top half(kips)
R3=R2                                                                          # Compressive forces on top half(kips)
#Bending Moment
M=2*((R1*7.5)+(R2*4.67))                                                       # Bending moment(kip.in)
#Radius of Curvature
p=round(((7/0.001724)*0.0833),0)                                               # Radius of curvature(ft)        

# Result
print ('Case(a) Bending moment = %lf kip.in' %My)
print ('Case(a) Radius of curvature = %.lf in' %pY)
print ('Case(b) Bending moment = %lf kip.in' %M)
print ('Case(b) Radius of curvature = %.lf ft' %p)
Case(a) Bending moment = 9525.000000 kip.in
Case(a) Radius of curvature = 4640 in
Case(b) Bending moment = 10226.342000 kip.in
Case(b) Radius of curvature = 338 ft

SAMPLE PROBLEM 4.6, Page number 264

In [29]:
import math


# Variable declaration
sY=240                                                                   # Yield strength(MPa)
A1=(0.1*0.02)                                                            # Area of cross section(m**2)
A2=(0.02*0.02)                                                           # Area of cross section(m**2)
A3=(0.02*0.06)                                                           # Area of cross section(m**2)
A4=(0.06*0.02)                                                           # Area of cross section(m**2)

# Calculation
#Neutral Axis
A=(100)*(20) + (80)*(20) + (60)*(20)                                     # Total area(mm**2)
y=(2400-((20)*(100)))/(20)                                               # Distance(mm)
#Plastic Moment
R1=(A1*sY*1000)                                                          # Resultant force(kN) 
R2=(A2*sY*1000)                                                          # Resultant force(kN)
R3=(A3*sY*1000)                                                          # Resultant force(kN) 
R4=(A4*sY*1000)                                                          # Resultant force(kN) 

Mp=(0.030*R1) + (0.010*R2) + (0.030*R3) + (0.070*R4)                     # Plastic moment(kN.m)

# Result
print ('Case(a) Plastic moment = %lf kN.m' %Mp)
Case(a) Plastic moment = 44.160000 kN.m

SAMPLE PROBLEM 4.7, Page number 265

In [30]:
import math

# Variable declaration
y=7                                                                                  # Distance(in)
s=-3.01                                                                              # Stress(ksi)

# Calculation
#Loading
M=10230                                                                              # Couple of moment(kip.in)
#Elastic Unloading
sMl=((10230)*(8))/(1524.0)                                                             # Maximum stress(ksi)
#Permanent Radius of Curvature
p=round(((7)*(29*pow(10,6))*(pow(10,-3)))/(3.01),-2)                                 # Permanent radius of curvature(in)
p=round((p*0.083333),-1)                                                             # Conversion(ft)

# Result
print ('Case(a) Residual stress = %lf ksi' %sMl)
print ('Case(a) Permanent radius of curvature = %lf ft' %p)
Case(a) Residual stress = 53.700787 ksi
Case(a) Permanent radius of curvature = 5620.000000 ft

SAMPLE PROBLEM 4.8, Page number 273

In [38]:
import math
from sympy import symbols

# Variable declaration
M=symbols('M')                                                                    # Variable declaration
P=symbols('P')                                                                    # Variable declaration

# Calculation
#Properties of Cross Section.
A=(3*pow(10,-3))                                                                  # Area of cross section(mm**2)
Y=0.038                                                                           # Distance(m)
I=868*(pow(10,-9))                                                                # Moment of inertia(m**4)
d=(0.038-0.010)                                                                   # Distance(m) 
#Force and Couple at C
M=0.028*P                                                                         # Equivalent force-couple system
s0=(P*(round((1/(3*pow(10,-3))),0)))                                              # Stress compression
s1=(P*round(((0.028)*(0.022))/(868*(pow(10,-9))),0))                              # Stress tension 
s2=(P*round(((0.028)*(0.038))/(868*pow(10,-9))))                                  # Stress compression
#Superposition
sA=-s0+s1                                                                         # Stress tension
sB=-s0-s2                                                                         # Stress compression
#Largest Allowable Force
P1=(30/377.0)*(1000)                                                                # The magnitude of P for which the tensile stress at point A is equal to the allowable tensile stress of 30 MPa is found by writing
P2=(120/1559.0)*(1000)                                                              # Determine the magnitude of P for which the stress at B is equal to the allowable compressive stress of 120 MPa.

# Result
print ('Largest allowable force = %lf kN' %P2)
Largest allowable force = 76.972418 kN

Example 4.09, Page number 285

In [42]:
import math

#Variable declaration
F=4.80                                                    # Couple(lb)
l=80                                                      # Length(mm)
b=120                                                     # Breadth(mm)
n=-1

#Calculation         
# Case(a)
Mx=(F)*(40)                                               # Stress in x(N.m)
Mz=(F)*(60-35)                                            # Stress in z(N.m)
A=(0.080)*(0.120)                                         # Area of cross section(m**2)
Ix=(1/12.0)*(0.120)*(0.080**3)                              # Centroidal moment of inertia(m**4)
Iz=(1/12.0)*(0.080)*(0.120**3)                              # Centroidal moment of inertia(m**4)
S0=n*(P/A)                                                # Stress(MPa) 
S1=((192)*(40))/(5.12*(pow(10,-6)))                       # Stress in x(MPa)
S2=((120)*(60))/(11.52*(pow(10,-6)))                      # Stress in y(MPa)
Sa=-0.5-1.5-0.625                                         # Stress at A(MPa)
Sb=-0.5-1.5+0.625                                         # Stress at B(MPa)
Sc=-0.5+1.5+0.625                                         # Stress at C(MPa) 
Sd=-0.5+1.5-0.625                                         # Stress at D(MPa)


# Case(b)
BG=((1.375)/(1.625+1.375))*80                             # Distance(mm)
HA=((2.625)/(2.625+0.375))*80                             # Distance(mm)

# Result
print ('Stress at A = %lf MPa' %Sa)
print ('Stress at B = %lf MPa' %Sb)
print ('Stress at C = %lf MPa' %Sc)
print ('Stress at D = %lf MPa' %Sd)
print ('Neutral axis BG = %lf mm' %BG)
print ('Neutral axis HA = %lf mm' %HA)
Stress at A = -2.625000 MPa
Stress at B = -1.375000 MPa
Stress at C = 1.625000 MPa
Stress at D = 0.375000 MPa
Neutral axis BG = 36.666667 mm
Neutral axis HA = 70.000000 mm

SAMPLE PROBLEM 4.9, Page number 287

In [43]:
import math
from sympy import symbols

# Variable declaration
P=symbols('P')


#Properties of Cross Section
A=7.46                                                                  # Area of cross section(in**2)
Sx=24.7                                                                 # Section moduli x(in**3)
Sy=2.91                                                                 # Section moduli y(in**3) 
#Force and Couple at C
Mx=4.75*P                                                               # Moment x()
My=1.5*P                                                                # Moment y()
#Normal Stresses
s1=P*round((1/7.46),4)                                                  # Normal stresses 
s2=P*round((4.75/24.7),4)                                               # Normal stresses
s3=P*round((1.5/2.91),4)                                                # Normal stresses
#Superposition
sA=-s1 + s2 +s3                                                         # Stress at A 
sB=-s1+s2-s3                                                            # Stress at B
sD=-s1-s2+s3                                                            # Stress at D
sE=-s1-s2-s3                                                            # Stress at E

P=12/0.842                                                              # Maximum compressive stress(kips) 

# Result
print ('Largest permissible load = %lf kips' %P)
Largest permissible load = 14.251781 kips

SAMPLE PROBLEM 4.10, Page number 288

In [7]:
import math
from sympy import sin,cos

# Variable declaration  
M0=1500                                                                 # Couple of magnitude(kN)
yA=50                                                                   # Distance() 
zA=74
Iy=(3.25*(pow(10,-6)))                                                  # Moment of inertia(m**4)
Iz=(4.18*(pow(10,-6)))                                                  # Moment of inertia(m**4)
Iyz=(2.87*(pow(10,-6)))                                                 # Moment of inertia(m**4)   

# Calculation
# Principal axes
Theta=(80.8)/2.0                                                          # Angle 
R=math.sqrt(pow(0.465,2)+pow(2.87,2))                                        # Radius
R=2.91*(pow(10,-6))                                                     # Converting to meter
Iu=3.72-2.91                                                            # Moment of inertia(m**4)
Iv=3.72+2.91                                                            # Moment of inertia(m**4) 
#Loading
Mu=(M0*sin(40.4))                                                       # Applied couple(N.m)             
Mv=(M0*cos(40.4))                                                       # Applied couple(N.m)
#Case(a) Stress at A
uA=50*cos(40.4*((2*math.pi)/360.0))+74*sin(40.4*((2*math.pi)/360.0))        # Perpendicular distances(mm)
vA=-50*sin(40.4*((2*math.pi)/360.0))+74*cos(40.4*((2*math.pi)/360.0))       # Perpendicular distances(mm)
sA=((972*0.0239)/(0.810*(pow(10,-6))) - ((1142)*(0.0860))/(6.63*pow(10,-6)))/(pow(10,6)) # Stress at A(MPa)
#Case(b) Neutral Axis
phy=81.8                                                                # Angle neutral axis with the v axis(degree) 
B=81.8-40.4                                                             # Angle neutral axis with the horizontal axis(degree)

# Result
print ('Stress at point A  = %lf MPa' %sA)
print ('The angle formed by the neutral axis and the horizontal is  = %lf degree' %B)
Stress at point A  = 13.866727 MPa
The angle formed by the neutral axis and the horizontal is  = 41.400000 degree

Example 4.10, Page number 298

In [11]:
import math
from sympy import integrate,symbols,log

#Variable declaration
r=6                                                    # mean radius(in)
b=2.5                                                  # Breadth(in) 
h=1.5                                                  # Height(in)
n=-1
b=symbols('b')                                         # Variable declaration
h=symbols('h')                                         # Variable declaration
r=symbols('r')                                         # Variable declaration
r1=symbols('r1')                                       # Variable declaration
r2=symbols('r2')                                       # Variable declaration


#Calculation         
# Case(a)
A=b*h                                                  # Area
R=A/(integrate((1/r), (r, r1, r2)))*(b)                # Radius
r1=6-((1/2.0)*(1.5))                                     # Inner radius(in)
r2=6+((1/2.0)*(1.5))                                     # Outer radius(in)
# Case(b)
R=(1.5)/(log(r2/r1))                                   # Distance(in)
e=6-R                                                  # Distance between the centroid and the neutral axis 

# Result
print ('The distance e between the centroid and the neutral axis of the cross section  = %lf in' %e)
The distance e between the centroid and the neutral axis of the cross section  = 0.031381 in

Example 4.11, Page number 299

In [22]:
import math
from sympy import symbols

#Variable declaration
M=8                                                                        # Bending moment(kip.in)
A=(2.5)*(1.5)                                                              # Area(in**2)  
R=5.969                                                                     
e=0.0314                                                                   # Distance(in)

#Calculation         
# Case(a)
Smax=((8)*(6.75-5.969))/((3.75)*(0.0314)*(6.75))                           # Maximum stress(ksi)  
Smin=((8)*(5.25-5.969))/((3.75)*(0.0314)*(5.25))                           # Minimum stress(ksi)

# Result
print ('Maximum stress = %lf ksi' %Smax)
print ('Minimum stress = %lf ksi' %Smin)
Maximum stress = 7.860974 ksi
Minimum stress = -9.304620 ksi

SAMPLE PROBLEM 4.11, Page number 300

In [13]:
import math
from sympy import integrate,symbols,solve

#Variable declaration
P=symbols('P')                                                             # Variable declaration  
s=symbols('s')                                                             # Variable declaration
A=symbols('A')                                                             # Variable declaration
M=symbols('M')                                                             # Variable declaration
R=symbols('R')                                                             # Variable declaration 
e=symbols('e')                                                             # Variable declaration

# Calculation
#Centroid of the Cross Section
r=(120*(pow(10,3)))/(2400.0)                                                 # Distance(m) 
#Force and Couple at D
M=((50+60)/1000.0)*P                                                         # Moment
#Superposition
s=-(P/A) + (M*(r-R))/(A*e*r)                                               # Total stress
#Radius of Neutral Surface
r=symbols('r')                                                             # Variable declaration 
R=(2400)/((integrate((80/r), (r, 30, 50)))+(integrate((20/r), (r, 50, 90))))# Radius of neutral surface(m) 
#Allowable Load
P=solve(-(P/(2.4*pow(10,-3))) + (0.110*P*(0.030-0.04561))/(2.4*(pow(10,-3))*(0.00439)*(0.030)) + (50*pow(10,6)),P) # Allowable load(kN)
P=round((P[0]/1000.0),2)                                                     # Rounding off(kN) 

# Result
print ('Largest force that can applied to the component = %lf kN' %P)                                 
Largest force that can applied to the component = 8.550000 kN