Breaking Ice with Air cushion Vehicle - Find minimum pressure with Newton-Raphson method

Objective: Writing a python script to solve for the minimum value of pressure for Breaking Ice with Air cushion Vehicle using the newton raphson method.

Problem Statement:

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

β = 0.5,

r = 40 feet and

σ = 150 pounds per square inch (psi).

we have to,

1. Use the newton raphson method to find out the value of pressure for h = 0.6 ft.
2. Find the optimal relaxation factor for this problem with the help of a suitable plot.
3. 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.

Solution Procedure: As we are using the newton raphson method.

The Newton-Raphson method is a way to quickly find a good approximation for the root of a real-valued function f ( x ) = 0 . It uses the idea that a continuous and differentiable function can be approximated by a straight line tangent to it.

According to Newton Raphson method,

X_new = X_guess - α(f(x_guess)/f'(x_guess)

As we have to find the minimum pressure , in above equation the X will be replaced by the pressure i.e. P_guess

"""Breaking Ice with Air cushion Vehicle - Find minimum pressure with Newton-Raphson method """    
import matplotlib.pyplot as plt  
import math  
import numpy as np      

#function to determine the minimum cushion pressure needed to break a given thickness of ice using an air cushion vehicle

  
def f(p,h):  	
	return (pow(p,3)*(1-pow(beta,2)))+(((0.4*h*pow(beta,2))-(sigma*pow(h,2)/pow(r,2)))*pow(p,2))+(((pow(sigma,2)*pow(h,4))/(3*pow(r,4)))*p)-(pow((sigma*pow(h,2)/(3*pow(r,2))),3))    

#derivative of function  
def fprime(p,h):  	
	return(3*(pow(p,2)*(1-pow(beta,2))))+(2*p*((0.4*h*pow(beta,2))-(sigma*pow(h,2)/pow(r,2))))+((pow(sigma,2)*pow(h,4))/(3*pow(r,4)))    

#inputs available

#width of the ice edge  
beta=0.5    
#size of the air cushion(feet)  
r=40    
#Tensile strength of the ice(psi)  
sigma=150    

height=[0.6,1.2,1.8,2.4,3.0,3.6,4.2]    

#cushion pressure  
p_guess=0.2       

#assuming a suitable relaxation factor (alpha) 
alpha=1

tolerance=1e-4   
pressure=[]    
iter=1  

for h in height:  	
	while(abs(f(p_guess,h))>tolerance):  		
		p_guess=p_guess-((f(p_guess,h)/fprime(p_guess,h))*alpha)  		
		iter=iter+1

	
	pressure.append(p_guess)

    

for i in range (0,7):  	
	if i==0: 
		print("      TABLE        ") 		
		print("height--->pressure")    	


	print(height[i],"--->" ,pressure[i])	  
	
print("Total numbers of iterations is=",iter)    
	
plt.plot(height,pressure)  
plt.xlabel('height')  
plt.ylabel('pressure')
plt.figure()
#by performing program for each alpha we get iter value (manually)
alpha=[0.2,0.4,0.6,0.8,1]
iter=[117,55,33,23,14]
plt.plot(alpha , iter)  
plt.xlabel('alpha')  
plt.ylabel('iter')
plt.show()

Results and Discussion:

output:

plots:

plot 2:

so,

As we increase the value of alpha the no. of iterations are decreasing.

and, so for get faster and desired results i have choosen the alpha value as 1.

i have also plotted the graph pressure vs. height and observed that we got minimum value of pressuree at 0.6. which was our objective.