Ice Breaking Challenge
Objective: The objective of this project is to determine the minimum pressure required to break a given thickness of ice using an air cushion vehicle through Newton Raphson Method
Equation:
`p^3(1-beta^2)+(0.4hbeta^2-(sigmah^2)/r^2)p^2+((sigma^2h^4)/(3r^4))p-((sigmah^2)/(3r^2))^3=0`
The following parameters shall be studied to obtain the final results
1. Part 1: Using `alpha=1` obtain the results as a test run
2. Part 2: Optimization of `alpha` for a particular thickness of ice i.e, finding the Optimum Relaxation Factor for h=0.6
3. Part 3: Obtaining the appropriate relaxation factor
Theory: Newton Raphson method is an effective way to solve non linear equations which cannot be solved usign conventional means. The method finds the roots of a given equation iteratively to achive the final value at which the function reduces to minimum.
For performing Newton Raphson the following equation was used
`x=x+alpha(f(x)/(f'(x)))`
where `alpha` is the relaxation factor (if `alpha`<1: Under Relaxation & `alpha`>1: Over Relaxation)
x is the quantity at which the function becomes minimum
For the target equation:
`f(p)=p^3(1-beta^2)+(0.4hbeta^2-(sigmah^2)/r^2)p^2+((sigma^2h^4)/(3r^4))p-((sigmah^2)/(3r^2))^3`
Now the 1st order derivative of the funtion: `f'(p)=3p^2(1-beta^2)+(0.4hbeta^2-(sigmah^2)/r^2)2p+((sigma^2h^4)/(3r^4))`
Now for performing Newton Raphson
`p_(guess)=p_(guess)+alpha(f(p_(guess))/(f'(p_(guess))))`
Program:
Part 1:
import matplotlib.pyplot as plt
b = 0.5 #Beta: 0.5
r = 40 #Size of Air Cushion in feet
s = 150*144 #sigma: tensile strength of Ice in psi (converted to pound per square foot by multiplying 144)
h = [0.6,1.2,1.8,2.4,3.0,3.6,4.2] # Multiple thickness of Ice field given in Feet (the multiple values are clubbed in an array for ease)
tol = 1e-5 #Considering the convergence criteria for the pressure equation
i = 1
p_guess = 100 #Providing an initial guess value for pressure
def f(p,h): #Defining the pressure equation for Air Cushion Vehicle
return ((pow(p,3)*(1-pow(b,2)))+ ((0.4*h*pow(b,2) - (s*pow(h,2)/pow(r,2)))*pow(p,2)) + (pow(s,2)*pow(h,4)*p/(3*pow(r,4))) - (pow((s*pow(h,2)/(3*pow(r,2))),3)))
def f_prime(p,h): #Defining the 1st order derivative (w.r.t pressure) function of the equation
return ((3*pow(p,2)*(1-pow(b,2)))+ 2*p*(0.4*h*pow(b,2) - (s*pow(h,2)/pow(r,2)))+ (pow(s,2)*pow(h,4)/(3*pow(r,4))))
min_pressure = [] #Defining a pressure array to store the mininum air cushion pressure to be obtained for respective thicknees of ice field
iter = [] #Defining an iteration array to store the values of iteration at which colvergence is obtained
for x in h: #Using a loop to feed in different values of ice field thickness to the pressure equation
alpha = 1 #calling the function optalpha for calculating the optimum value of alpha for each values of ice field thickness
i=1
while(abs(f(p_guess,x))>tol): #Using an exit controlled loop which would terminate when the convergence criteria is met by performing Newton Raphson method
p_guess = p_guess - alpha * (f(p_guess,x)/f_prime(p_guess,x))
i=i+1
iter.append(i)
min_pressure.append(p_guess)
p_guess=100
#Printing the results in a tabular format
for j in range(0,len(h)):
if j == 0:
print('Ice Thickness (feet) | ', 'Min. Air Cushion Pressure (psf) | ', 'No. of Iteration for Convergence |')
print(' ', h[j], ' ', min_pressure[j], ' ', iter[j])
#Plotting the graph for min. air cushion pressure vs thickness of Ice Field
plt.plot(h,min_pressure)
plt.xlabel('Ice Feild Thickness (feet)')
plt.ylabel('Min. Air Cushion Pressure (psf)')
plt.show()
#END of ProgramResults:
Output Table: Minimum Pressure by Ice Thickness (for `alpha`=1)
Inference: The above values depict say for ice thickeness of h=0.6 feet, the minimum pressure of air cushion vehicle required to break it is 4.239064120503728 psf and the program took 14 iterations for achieving this minimum pressure. The same inference can be derived for other ice thickness values as well. The below plot of minimum air cushion pressure required for varing ice thickness
Fig: Minimum Pressure by Ice Thickness (for `alpha`=1)
Part 2: Optimization of `alpha` for a particular thickness of ice i.e, finding the Optimum Relaxation Factor for h=0.6.
This program is a test program for finiding the optimum `alpha` value for a given thickness
import matplotlib.pyplot as plt
import numpy as np
b = 0.5
r = 40
s = 150*144
h = 0.6
tol = 1e-5
def f(p):
return ((pow(p,3)*(1-pow(b,2)))+ ((0.4*h*pow(b,2) - (s*pow(h,2)/pow(r,2)))*pow(p,2)) + (pow(s,2)*pow(h,4)*p/(3*pow(r,4))) - (pow((s*pow(h,2)/(3*pow(r,2))),3)))
def f_prime(p): #Defining the 1st order derivative (w.r.t pressure) function of the equation
return ((3*pow(p,2)*(1-pow(b,2)))+ 2*p*(0.4*h*pow(b,2) - (s*pow(h,2)/pow(r,2)))+ (pow(s,2)*pow(h,4)/(3*pow(r,4))))
a = np.linspace(0.1,1.9,20)
pressure = []
count = []
p_guess = 100
for x in a:
i=1
while(abs(f(p_guess))>tol): #Using an exit controlled loop which would terminate when the convergence criteria is met by performing Newton Raphson method
p_guess = p_guess - x * (f(p_guess)/f_prime(p_guess))
i=i+1
pressure.append(p_guess)
p_guess = 100
count.append(i)
plt.plot(a,count)
plt.xlabel('Alpha')
plt.ylabel('No. of Interations')
plt.show()
plt.plot(a,pressure)
plt.xlabel('Alpha')
plt.ylabel('Pressure (psf)')
plt.show()
Results
Fig: No. of Iterations by `alpha`(for h=0.6)
Inference: the fig above denotes the number of iterations for individual values of `alpha` ranging from 0.1 to 1.9 considering ice thickness constant(h=0.6). It is observed that the number of iterations reduce with increase in `alpha` and after a certian point it increses again.
Fig: Test Pressure by `alpha` (for h=0.6)
Inference: the fig above denotes the test pressure values for each values of `alpha` and it is observed that the values fluctuate with change in alpha and after point (`alpha=1.6`) we get the minimum pressure (4.23906267904548 psf) and after that the value increases again till the end `alpha` value and there is no significant drop below this point again.
In order to the check whether the value of `alpha` changes with change in ice thickness value the program was build again but with h=1.2 and the following results were obtained
Fig: No. of Iterations by `alpha` (for h=1.2)
Inference: similar phenomenon was observed in the no. of iterations vs `alpha` values as in case of h=0.6
Inference: the graph obtained for h=1.2 for pressure vs `alpha` is different than that observed for =0.6. So it can be inferred that the value of `alpha` for minimum pressure changes with change in ice thickness as well.
In order to capture more accurate results the optimum `alpha`value has to be calculated for each ice thickness value
Part 3: Obtaining the appropriate relaxation factor
import matplotlib.pyplot as plt
import numpy as np
b = 0.5 #Beta: 0.5
r = 40 #Size of Air Cushion in feet
s = 150*144 #sigma: tensile strength of Ice in psi (converted to pound per square foot by multiplying 144)
h = [0.6,1.2,1.8,2.4,3.0,3.6,4.2] # Multiple thickness of Ice field given in Feet (the multiple values are clubbed in an array for ease)
tol = 1e-5 #Considering the convergence criteria for the pressure equation
i = 1
p_guess = 100 #Providing an initial guess value for pressure
def f(p,h): #Defining the pressure equation for Air Cushion Vehicle
return ((pow(p,3)*(1-pow(b,2)))+ ((0.4*h*pow(b,2) - (s*pow(h,2)/pow(r,2)))*pow(p,2)) + (pow(s,2)*pow(h,4)*p/(3*pow(r,4))) - (pow((s*pow(h,2)/(3*pow(r,2))),3)))
def f_prime(p,h): #Defining the 1st order derivative (w.r.t pressure) function of the equation
return ((3*pow(p,2)*(1-pow(b,2)))+ 2*p*(0.4*h*pow(b,2) - (s*pow(h,2)/pow(r,2)))+ (pow(s,2)*pow(h,4)/(3*pow(r,4))))
def optalpha(H): # defining a function which would return the optimum values of alpha for invidual values of h (ice field thickness)
a = np.linspace(0.1,1.9,20)
test_pressure = []
p_guess=100
for x in a:
i=1
while(abs(f(p_guess,H))>tol): #Using an exit controlled loop which would terminate when the convergence criteria is met by performing Newton Raphson method
p_guess = p_guess - x * (f(p_guess,H)/f_prime(p_guess,H))
i=i+1
test_pressure.append(p_guess)
p_guess=100
ind=test_pressure.index(min(test_pressure))
return(a[ind])
min_pressure = [] #Defining a pressure array to store the mininum air cushion pressure to be obtained for respective thicknees of ice field
iter = [] #Defining an iteration array to store the values of iteration at which colvergence is obtained
opt_alpha = []
for x in h: #Using a loop to feed in different values of ice field thickness to the pressure equation
alpha = optalpha(x) #calling the function optalpha for calculating the optimum value of alpha for each values of ice field thickness
i=1
while(abs(f(p_guess,x))>tol): #Using an exit controlled loop which would terminate when the convergence criteria is met by performing Newton Raphson method
p_guess = p_guess - alpha * (f(p_guess,x)/f_prime(p_guess,x))
i=i+1
iter.append(i)
min_pressure.append(p_guess)
opt_alpha.append(alpha)
p_guess=100
#Printing the results in a tabular format
for j in range(0,len(h)):
if j == 0:
print('Ice Thickness (feet) | ', 'Min. Air Cushion Pressure (psf) | ', 'No. of Iteration for Convergence | ', 'Optimum Alpha')
print(' ', h[j], ' ', min_pressure[j], ' ', iter[j], ' ' ,round(opt_alpha[j],2))
#Plotting the graph for min. air cushion pressure vs thickness of Ice Field
plt.plot(h,min_pressure)
plt.xlabel('Ice Feild Thickness (feet)')
plt.ylabel('Min. Air Cushion Pressure (psf)')
plt.show()
#END of ProgramOutput Table: Minimum Pressure by Ice Thickness
Inference: The above values depict say for ice thickeness of h=0.6 feet, the minimum pressure of air cushion vehicle required to break it is 4.2398626790548 psf and the program took 31 iterations for achieving this minimum pressure with `alpha` value of 1.62. The same inference can be derived for other ice thickness values as well.
The below plot of minimum air cushion pressure required for varing ice thickness
Fig: Minimum Pressure by Ice Thickness
Conclusions:
1. The number of iterations for convergence of the results (min. pressure) decreases with increase in the value of `alpha` and after a certain point it starts increasing again
2. The minimum pressure required to break the a particular thickness of ice using air cushion vehicle increases with the increase in the ice feild thickness
3. The value of `alpha` changes with change in ice feild thickness and in order to capture accurate results the optimum `alpha` needs to be calculated separately for individual ice feild thickness values.