Newton Raphson Method- Minimum Pressure to Break Ice with an Air cushion Vehicle using PYTHON

Abstract - The ice breaking mechanism under a slow moving air cushion vehicle is explored (E R Muller - Mathematical modelling : Classroom notes in applied mathematics -page(29)) . Determining the minimum cushion pressure needed to break a given thickness of ice . Solving by using the newton-raphson method in python . Iterating towards the solution for different thickness of ice and understanding optimum relaxation factor for a particular case. 

 Introduction :

Air Cushion vehicle : 

Recent experiments using self-powered air cushion vehicles have demonstrated that they can be used successfully to break up ice on rivers. This is particularly useful at ice break-up time, when the ice may pack in one section of the river hindering the flow of water and causing floods above the packed ice. Other experiments have used air cushion platforms tied to the front of ships. Such an experiment near Thunder Bay, Ontario, Canada, has demonstrated that an air cushion platform can be a relatively efficient means of ice breaking compared to conventional displacement vessels with their large power and crew requirements. 

                               (ref - https://dfdickins.com/our-projects/air-cushion-vehicles)

Newton-Raphson method :

It is a root finding algorithm which produces successively better approximations to the roots of real valued function. Generally we start with an initial educated guess value and the newtons method tells us the better approximation for the root. This process maybe repeated as many times as necessary to get the desired accuracy. In our case, we set a certain tolerance as condition to stop the iterations.

Objectives of the study : 

1. Implementing the study using the mentioned iterative technique in python.

2. Finding Minimum pressure to break certain bunch of given thickness of ice with relaxation factor 1. 

3. Understanding the  choice of an appropriate optimum relaxation factor for the method. 

4. Explaining results w.r.t appropriate relaxation factor. 

Mathematical equation : 

(E R Muller - Mathematical modelling : Classroom notes in applied mathematics -page(34))

`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`

where, 

p = cushion pressure 

h = thickness of ice field

r = size of air cushion

sigma = tensile strength of ice

beta = related to width of ice wedge

Applying Newton-raphson method to the equation to find minimum Pressure :  

Initial conditions and given parameters : 

1. p = 10 (initial guess value)

2. tolerance = 1e-4. 

3. alpha = 1 (for first case)

Given Parameters : 

1. beta = 0.5 

2. r = 40 feet

3. sigma = 150 psi

4. h = [0.6, 1.2, 1.8, 2.4, 3.0, 3.6, 4.2]

Newton-raphson algorithm: 

`p_(n ew) = p_(guess ) - alpha*(f(p_(guess ))/(f'(p_(guess ))))`

we take derivative of the equation and feed it in the algorithm where f' is mentioned. 

This algorithm iterates until the set tolerance limit is achieved. 

Python code: 

''' 
code by Aadil shaikh
10/9/19

Newton raphson method - Finding minimum pressure for
breaking ice with air cushion vehicle 

'''

import matplotlib.pyplot as plt 
from tabulate import tabulate
import numpy as np 


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

def fprime(p):
	term1 = (0.4*h*pow(beta,2) - ((sigma*pow(h,2))/(pow(r,2))))
	term2 = ((pow(sigma,2)*pow(h,4))/(3*pow(r,4)))
	term3 = pow(((sigma*pow(h,2))/(3*pow(r,2))),3)
	return 3*pow(p,2)*(1-pow(beta,2)) + 2*p*term1 + term2 


# Performing newton iteraion 

# given values 
beta = 0.5
r = 40 # feet
sigma = 150*144 
h1 = [0.6,1.2,1.8,2.4,3.0,3.6,4.2]


# iteration values
p_guess = 10
alpha = 1  #  relaxation parameter



tol = 1e-4
iter = 1

Data = []
h_h = []
iteration = []

for i in range(len(h1)):
	h = h1[i]
	while (abs(f(p_guess)) > tol):
		p_guess = p_guess - alpha*(f(p_guess)/fprime(p_guess))
		iter = iter + 1
	#print(iter)

	Data.append(p_guess)
	h_h.append(h)
	iteration.append(iter)
	

fig = plt.figure(dpi=200)
ax = fig.add_subplot(1,1,1)
table_data=[
	["Thickness of Ice field- h (feet)", "Pressure (psft)" , "Iteration"],
	[h_h[0],Data[0],iteration[0]],
    [h_h[1] ,Data[1],iteration[1]],
    [h_h[2] ,Data[2],iteration[2]],
    [h_h[3] ,Data[3],iteration[3]],
    [h_h[4] ,Data[4],iteration[4]],
    [h_h[5] ,Data[5],iteration[5]],
    [h_h[6] ,Data[6],iteration[6]]]
table = ax.table(cellText=table_data, loc='center')
table.set_fontsize(15)
table.scale(1,2)
plt.title('Pressure values for increasing thickness of ice field & no. of iterations')
ax.axis('off')

plt.figure(2)
plt.plot(h_h,Data,color="red",linewidth=3)
plt.xlabel('Thickness of Ice field- h (feet)',fontsize = 'large' , fontweight = 'bold')
plt.ylabel('Presure (psft)',fontsize = 'large' , fontweight = 'bold')
plt.title('Relation between Pressure and Thickness of ice field')
plt.grid()
plt.show()	

The equation is broken into 3 terms and typed in the code. 

For the values of given thickness of ice (h) we are calculating the minimum pressure required.  The newton-raphson algorithm is inside the while loop iterating till the tolerance is achieved. 

A table is constructed to display results and the graph shows the variation of minimum pressure for the given thickness array of ice.

Output : -

Results : 

1. Using relaxation factor 1 for above case, we find the minimum pressure values for each increasing thickness of ice. 

2. we observe as the thickness of ice increases, the minimum pressure required to break the ice also increases exponentially.

3. From the plot, we can observe how the curve is increasing in a non linear fashion demonstrating the values obtained above. 

Python code for objective 3 : 

import matplotlib.pyplot as plt 
import numpy as np 


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

def fprime(p):
	term1 = (0.4*h*pow(beta,2) - ((sigma*pow(h,2))/(pow(r,2))))
	term2 = ((pow(sigma,2)*pow(h,4))/(3*pow(r,4)))
	term3 = pow(((sigma*pow(h,2))/(3*pow(r,2))),3)
	return 3*pow(p,2)*(1-pow(beta,2)) + 2*p*term1 + term2 


beta = 0.5
r = 40
h = 0.6
sigma = 150*144
p_guess = 10
tol = 1e-4

iter = 1 

a1 = []
itera = []
for a in range(1,200):
	a = a*0.01
	
	iter = 1
	p_guess = 10
	while (abs(f(p_guess)) > tol):
		p_guess = p_guess - a*(f(p_guess)/fprime(p_guess))
		iter = iter + 1
	
	itera.append(iter)
	a1.append(a)


plt.plot(a1,itera,color="green",linewidth=3)
plt.xlabel('Alpha',fontsize = 'large' , fontweight = 'bold')
plt.ylabel('No. of Iterations',fontsize = 'large' , fontweight = 'bold')
plt.savefig('alpha.png')
#plt.show()

 This code is written with thickness as 0.6 (assumption) to study what could be the optimum relaxation factor and which one should be chosen to converge quickly and accurately. Its scripted for alpha (relaxation factor) from 0-2 . 

Output : 

Results : -

1. From observation of plot and code output, we can see that it takes almost same number of iterations(least) to converge for a relaxation factor of 0.78 - 1.35 .

2. For a more precise observation, the iteration values are very closeby to each other varying with just 1-2 number of iterations for a range of 0.82-1.15. 

3. The newton-raphson technique, when using with under-relaxed condition, will deliver more accurate results, even if it takes more iteration or simulation time.

4. Whereas over relation can fasten the process but iterations could increase.

5. But observing on here the appropriate range of alpha comes out to be 0.82-1.15 to get the best result. 

6. On observing the iteration values, one iteration converged at 7 which was the least. It was for relaxation factor 1.16. The for loop was removed and solution was iterated for a = 1.16. 

7. 

8. so 1.16 came out to be the most optimum relaxation factor for this case. 

References : 

1. http://www.math.utah.edu/~cherk/teach/5740MathModeling/12mathmodel/sources/Mathematical%20Modelling,%20Classroom%20Notes%20in%20Applied%20Mathematics%20(Murray%20S.%20Klamkin)%200898712041.pdf

2. https://en.wikipedia.org/wiki/Hovercraft

3. https://dfdickins.com/our-projects (pic credit)

4. https://brilliant.org/wiki/newton-raphson-method/

                                               THE END 

 keywords - PYTHON, NUMERICAL-ANALYSIS, CFD, PYTHON-BASICS