Breaking Ice with Air cushion Vehicle - Find minimum pressure with Newton-Raphson method
1. AIM:
To write a python script to solve for the minimum value of pressure using the newton raphson method.
2. GOVERNING EQUATION OF ANY:
- `x_{n ew}=x_{guess}-[frac {f(x_{guess})}{f'(x_{guess})}]`.............(1)
- It is the equation of newton raphson method for finding the root of any given function (in this case, function of x).
3. PROBLEM STATEMENT AND OBJECTIVES:
- To write a python script to solve for the minimum value of pressure using the newton raphson method.
`p^3(1-b^2)+(0.4hb^2-(sigmah^2)/r^2)p^2+(sigma^2h^4)/(3r^4)p -((sigmah^2)/(3r^2))^3=0`
where p denotes the cushion pressure, h the thickness of the ice field, R the size of the air cushion, σ the tensile strength of the ice and β is related to the width of the ice wedge.
Take b = 0.5, r = 40 feet and σ = 150 pounds per square inch (psi).
- Use the newton raphson method to find out the value of pressure for h = 0.6.
- Find the optimal relaxation factor for this problem with the help of a suitable plot.
- Tabulate the results of p for h = [0.6,1.2,1.8,2.4,3,3.6,4.2] assuming a suitable relaxation factor.
- Plot the graph of thickness and pressure.
4. FLOW OF PROJECT:
- Define a function of p having the main function.
- Define new function having the derivative of the main / parent function.
- Define new functions according to the objectives
- Call the functions in main body.
- Plot the graphs and obtain conclusions.
5. PYTHON CODE:
import matplotlib.pyplot as plt
import numpy as np
#Defining the main function of x in a function named f
def f(p,b,h,r,sigma):
return (pow(p,3)*(1-pow(b,2)) + (0.4*h*pow(b,2)-(sigma*pow(h,2))/pow(r,2)*pow(p,2) + (pow(sigma,2)*pow(h,4)*p)/(3*pow(r,4)) - pow((sigma*pow(h,2))/(3*pow(r,2)),3)))
#Defining the derivative of the function f
def fprime(p,b,h,r,sigma):
return 3*pow(p,2)*(1-pow(b,2)) + (0.4*h*pow(b,2)-(sigma*pow(h,2))/pow(r,2)*2*p + (pow(sigma,2)*pow(h,4))/(3*pow(r,4)))
#Defining function named 'only' which is used for CASE - I
def only(f,fprime,p_guess,iter):
while (abs(f(p_guess,b,h,r,sigma))>tol): #Running the loop til function value gets lesser than tolerance value
p_guess = p_guess - alpha*(f(p_guess,b,h,r,sigma)/fprime(p_guess,b,h,r,sigma)) #Newton raphson formula
iter = iter + 1 #Iteration counter
print('Value of pressure at h=',h,'is',p_guess) #Displaying the outputs
print('The iteration number is',iter)
print('The value of function is',f(p_guess,b,h,r,sigma),'\n')
#Defining function named 'listh' which is used for CASE - III
def listh(f,fprime,iter,p_guess):
h = [0.6,1.2,1.8,2.4,3.0,3.6,4.2] #Creating a list of thickness
print(' h ',' Pressure \n')
for i in h: #Loop will run for every i in h i.e i will have values of h
while (abs(f(p_guess,b,i,r,sigma))>tol):
p_guess = p_guess - alpha*(f(p_guess,b,i,r,sigma)/fprime(p_guess,b,i,r,sigma))
iter = iter + 1
pressure.append(p_guess) #Adding values of p_guess to empty array of pressure
print(i,' ','{:0.4f}'.format(p_guess),'\n') #Displaying in command window
plt.figure() #Open a new figure window
plt.plot(h,pressure) #Plotting a graph of thickness and pressure
plt.xlabel('Thickness (H)')
plt.ylabel('Pressure')
plt.grid(True)
plt.show()
#Defining function named 'alphaopti' which is used for CASE - II
def alphaopti(f,fprime,iter,p_guess,alpha):
print('Alpha value ','Iteration #\n')
iteration = [] #Creating empty array for iteration value storing
fiter = 10000 #Basically any value than initial iteration value for lowest alpha
falpha = [] #Empty array for storing critical alpha value
for i in alpha: #i will have values of alpha while looping
iter = 1
p_guess = 20
while (abs(f(p_guess,b,h,r,sigma))>tol):
p_guess = p_guess - i*(f(p_guess,b,h,r,sigma)/fprime(p_guess,b,h,r,sigma))
iter = iter + 1
if fiter<iter: #If previous iteration count is lesser than current previous count
falpha.append(j)
j = i #Then add that alpha value to its empty array
fiter = iter
print(' {:0.2f}'.format(i),' ',iter)
iteration.append(iter)
minimum = min(iteration) #Find minimum value of iteration count in the iteration array
plt.figure()
plt.plot(alpha,iteration) #Plotting values of alpha and respective iteration values
plt.plot(falpha[0],minimum,'o') #Plotting the lowest point
plt.xlabel('Alpha')
plt.ylabel('Iteration Count')
plt.grid(True)
plt.text(falpha[0],minimum,'{:0.2f} {} '.format(falpha[0],minimum))
plt.show()
#Newton raphson iteration
b = 0.5 #Related to width of ice wedge
r = 40 #Size of air cushion
sigma = 150*144 #tensile strength, pounds/feet^2
h = 0.6 #Thickness of ice field
p_guess = 20 #Initial guess of pressure
alpha = 0.4 #initial alpha values
tol = 1e-4 #Tolerance value
iter = 1
pressure = []
z1 = only(f,fprime,p_guess,iter) #Calling function f for CASE - I
iter = 1
p_guess = 20
alpha = np.arange(0.02,1.98,0.02) #Creating equal spaced array for CASE - II
print('\n\n')
z2 = alphaopti(f,fprime,iter,p_guess,alpha) #Calling function alphaopti for CASE - II
iter = 1
p_guess = 20
alpha = 0.4
print('\n\n')
z3 = listh(f,fprime,iter,p_guess) #Calling function listh for CASE - IIICode explaination:
We have defined function 'f' which contains the given equation in problem statement.
A new function 'fprime' is created which contains the derivative of the parent function.
- function 'only'.....CASE - I
- This function is used to find the value of pressure at a given thickness 'h'.
- We are using while loop till the value of the main function becomes less than the tolerance value. We are using newton raphson governing equation to evaluate the value of pressure in the while loop.
- function 'listh'....CASE - III
- This function is used to find values of pressures at different thickness 'h'.
- We first have created a list of thicknesses and then in for loop, for every value of thickness 'h', while the value of main function is greater than tolerance value, i.e until function value gets less than tolerance, the newton raphson will do iterations. Once, it has found a root, i.e value of main function is approximately 0, the next value of thickness 'h' is taken and again for loop will run.
- Then plots of thickness and pressure is plotted.
- function 'alphaopt'......CASE - II
- This function is used to find the minimum value of iteration and its corresponding value of alpha.
- It does so by checking every current value of no. of iterations with previous value of no. of iterations obtained. So basically if current iterations is less than previous iteration count then still scope of improvement is there. Once current iteration count is more than previous one, it means previous iteration count is the lowest and and corresponding to that, alpha value is found by storing only those alpha values by appending it to an empty array.
- Plot of values of alpha and iteration count is plotted and lowest point is marked.
- In the rest of the code, we have defined the given parameters like b, r, sigma (convert psi to psf), initial pressure guess, alpha, tolerance value (tol).
- And in the end, we have called the respective functions according to the case.
6. OBSERVATIONS:
- The first graph shows the parabolic nature of pressure changing as per change in thickness of ice field.
- The second figure shows the point where lowest iteration count was obtained.
7. CONCLUSION:
- We found out the roots of the main equation by using Newton-Raphson Method.
- The alpha value for minimum iteration count (9) is