PYTHON Project: Breaking Ice (Air Cushion Vehicle)

Aim
To write a PYTHON program to solve for the minimum value of pressure using the Newton-Raphson Method.

Theory
The Newton-Raphson method is a method for finding a good approximation for the root of a real-valued function `f(x) =0` rapidly [1]. It is premised on the idea that a tangent to a continuous and differentiable function would approach it.

Consider the situation where you need to discover the root of a continuous, differentiable function `f(x)`, and you know the root you're looking for is near `x=x_0`. Then, according to Newton's approach, a better approximation for the root is,
`x_1=x_0-(f(x_0))/(f'(x_0)`

This procedure can be done as many times as necessary to get the appropriate level of precision. In general, the next value for every `x` value `x_n` is,
`x_(n+1)=x_n-(f(x_n))/(f'(x_n)`

Figure 1: Pictorial demonstration of Newton-Raphson Method

The tangent line of the function `f(x)` at the point `x=x_n` has the slope `f'(x_n)` and goes through the point `(x_n,f(x_n))`. Therfore, this line has the equation, `y=f'(x_n)(x-x_n)+f(x_n)`. Setting `y=0` and `x=x_(n+1)` for our new approximation, we can now identify the root of this tangent line. The new approximation is obtained by solving is,
`x_(n+1)=x_n-(f(x_n))/(f'(x_n))`

When the relaxation factor is considered,
`x_(n+1)=x_n-(alpha*f(x_n))/(f'(x_n))`
where, `alpha` is the relaxation factor.

An Air-Cushion Vehicle is a craft with the whole or a significant part of its weight supported, whether at rest or in motion, by a continuously generated cushion of air with its effectiveness dependent on the proximity of the surface over which the craft operates [2]. A heavy ACV traveling at a crucial speed on the float ice sheet might break it into little slabs. Ice-breaking using ACV is a new inland river ice-breaking method that has a significant applicability in preventing ice jam floods in the Yellow River. An ACV can easily break the ice in flexure by forming depressions and sending cracks through the ice, depending on the thickness of the ice. Once a break is created, the ACV's air curtain follows the path of least resistance beneath the ice, reducing the supporting water beneath the ice. After that, the ice collapses under its own weight. At high speeds, an ACV can create a wave in the ice and cause cracks to appear. This would behave similarly to the dynamics of a truck driving on an ice road.

The formula for ice breaking is as follows:

`p^3(1-beta^2)+(0.4*hbeta^2-(sigmah^2)/r^2)p^2+(sigma^2h^4)/(3r^4)p-(sigmah^2)/(3r^2)=0`
where, `p` is the cushion [ressure, `h` is the thickness of the ice field, `r` is the size of the air cushion, `sigma` is the tensile strength of the ice, and `beta` is related to the width of the ice wedge.

Here, we are considering certain values: `beta=0.5`, `r=40` feet, and `sigma = 150` pounds per square inch (psi).

Objectives

1. Understanding the Newton-Raphson Method.
2. Understanding the concept of Air Cushion Vehicle.
3. Develop a PYTHON program to solve by Newton-Raphson Method.
4. Find the value of pressure for h=0.6 feet using the developed program.
5. Construct a PYTHON program to find the optimal relaxation factor with the help of a suitable plot.
6. Develop a PYTHON program to 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.

Program to Calculate the Cushion Pressure for h=0.6 feet
Programming is done in Python

# Breaking Ice with Air Cushion Vehicle: Newton-Raphson Method
# To calculate the pressure at h = 0.6 feet

import math
import matplotlib.pyplot as plt 
import numpy

# Inputs: given
beta = 0.5	#is related to width of ice edge
r = 40		#is the size of air cushion in feets
sigma = 150*144	#is the tensile strength of the ice (pounds/sq.inch) // 1psi=144psf
h = 0.6		#is the thickness of ice in feet

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

# Defining derivative of the function
def func_der(p):
	term5 = 3*pow(p,2)*(1-pow(beta,2))
	term6 = 2*p*((0.4*h*pow(beta,2))-(sigma*pow((h/r),2)))
	term7 = (pow(sigma,2)*pow((h/r),4))/3
	return (term5+term6+term7)

# Newton_Raphson Method
p_guess = 30	#initial guess to start Newton-Raphson Method
alpha = 0.5	#relaxation factor, random choice
tol = 1e-5	#setting tolerance value

# Interation count
ite = 1		#setting initial value


# Initiating Iteration
while(abs(func(p_guess))>tol):
	p_guess = p_guess - (alpha*(func(p_guess)/func_der(p_guess)))
	ite = ite+1
	
print('The value of pressure for h = 0.6 feet is',p_guess,'psf')

Each step of the program is explained with comments.

Further Explanation
Here, the relaxation factor (alpha) is assumed to be 0.5. The more the relaxation value closer to zero, the lesser will be difference for the new p_guess formed after each iteration and more will be the total iteration required. The tensile strenth value is multiplied with 144 to convert it from psi to psf.

Output of the Program

The Cushion Pressure value obtained to be 4.239 psf for thickness of ice of 0.6 ft.

Program to Find the Optimal Relaxation Factor with the help of a Suitable Plot
The programming is done in PYTHON

# Breaking Ice with Air Cushion Vehicle: Newton-Raphson Method
# To find appropriate relaxation factor

import math
import matplotlib.pyplot as plt 
import numpy as np

# Inputs: given
beta = 0.5	#is related to width of ice edge
r = 40		#is the size of air cushion in feets
sigma = 150*144	#is the tensile strength of the ice (pounds/sq.inch) // 1psi=144psf
h = 0.6		#is the thickness of ice in feet

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

# Defining derivative of the function
def func_der(p):
	term5 = 3*pow(p,2)*(1-pow(beta,2))
	term6 = 2*p*((0.4*h*pow(beta,2))-(sigma*pow((h/r),2)))
	term7 = (pow(sigma,2)*pow((h/r),4))/3
	return (term5+term6+term7)


# Newton_Raphson Method
alpha = np.linspace(0.01,1.99,500)#relaxation factor varied from 0.01 to 1.99
tol = 1e-5			  #setting tolerance value

# Interation count
tot_itr = []			  #for counting total iteration took for each alpha


# Initiating Iteration
for i in range(0,len(alpha)):
	p_guess = 10		  #initial guess to start Newton-Raphson Method for each iteration
	itr = 0			  #setting initial value for counting iteration for each alpha
	
	while(abs(func(p_guess))>tol):
		p_guess = p_guess - (alpha[i]*(func(p_guess)/func_der(p_guess)))
		itr = itr+1
	
	tot_itr.append(itr) 

# Finding Optimum Relaxation Factor
low_itr = min(tot_itr)			#finding the loop with min iteration		
pos_low_itr = tot_itr.index(low_itr)	#finding position of loop with min iteration
opt_alpha = alpha[pos_low_itr]

print("The optimum relation factor is ",opt_alpha, "with total iteration of", low_itr)

plt.figure()
plt.plot(alpha,tot_itr)
plt.plot(opt_alpha,low_itr,'*r', label= 'Position of Optimal Relaxation Factor')
plt.legend(loc='best')
plt.xlabel('Relaxation Factor')
plt.ylabel('No. of Iterations')
plt.show()

Each step of the program is explained with comments.

The relaxation parameter is considered from 0.01 to 1.99 using the linspace command. For these various relaxation factor, the number of iteration required is found. The relaxation factor with minimum iterations taken is the optimum value. In this program, first the minimum value of iteration is found using the min command, then the position of that respective iteration is found from the array/list, from which the optimum relaxation factor is identified.

Output of the Program

The optimum relaxation factor was found to be 0.998016 which requires only seven iteration for the solution. This optimum value is found for the thickness of ice as 0.6 feet.

The above diagram shows the number of iterations required for each relaxation factor. From the graph, the optimum value found is marked on the graph using the star symbol as shown in figure. The graph depicts that the number of iteration is observed to be high near 0.01 and 1.99, whereas it requires only less number of iterations for in the range of 0.25 to 1.75. The number of iterations is seemed to be overshooting when the relaxation factor value approaches two, and the loop is seemed to run infinitely without meeting the required tolerance.

Program to 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

The programming is done in PYTHON

# Breaking Ice with Air Cushion Vehicle: Newton-Raphson Method
# To 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

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

# Inputs: given
beta = 0.5	#is related to width of ice edge
r = 40		#is the size of air cushion in feets
sigma = 150*144	#is the tensile strength of the ice (pounds/sq.inch) // 1psi=144psf
h = [0.6, 1.2, 1.8, 2.4, 3, 3.6, 4.2]	#is the thickness of ice in feet

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

# Defining derivative of the function
def func_der(p,h):
	term5 = 3*pow(p,2)*(1-pow(beta,2))
	term6 = 2*p*((0.4*h*pow(beta,2))-(sigma*pow((h/r),2)))
	term7 = (pow(sigma,2)*pow((h/r),4))/3
	return (term5+term6+term7)


# Newton_Raphson Method
p_guess = 30	#initial guess to start Newton-Raphson Method
alpha = 0.5	#relaxation factor, random choice
tol = 1e-5	#setting tolerance value

# Interation count
itr = 1		#setting initial value

# Setting up array to store values
pressure = []	#to store min cushion pressure for each h
tot_itr = []	#to store the total iteration required for finding min p for each h

# Initiating Iteration
for i in range(0,len(h)):
	
	while(abs(func(p_guess,h[i]))>tol):
		p_guess = p_guess - (alpha*(func(p_guess,h[i])/func_der(p_guess,h[i])))
		itr = itr+1
	
	tot_itr.append(itr)
	pressure.append(p_guess)


# Tabulation 
contents = tuple(zip(h,tot_itr,pressure))
heading = ['Thickness of Ice (ft)', 'Iterations Taken','Minimum Cushion Pressure (psf)']
Table = tabulate(contents, headers=heading, tablefmt='fancy_grid')
print(Table)

Each step of the program is explained with comments. The table is created using the zip [3] and tuple [4] command, where the tabular format is set as fancy grid. All the calculations is done for the relaxation factor of 0.5.

Output of the Program

The table is created in PYTHON as shown in above figure, which gives the information of minimum cushion pressure and iterations required for various thickness.

It is found that, as thickness of the ice increases, the program requires more iterations with the same relaxation factor. Also, the cushion pressure increases with the ice thickness.

Conclusion
The PYTHON Program to solve the Newton-Raphson method is constructed successfully. The minimum cushion pressure required ice thickness of 0.6 feet is found. An optimum relaxation factor is found using the plot, later which the minimum cushion pressure required for various ice thickness is found and tabulated in PYTHON.

Reference
[1] https://brilliant.org/wiki/newton-raphson-method/
[2] https://www.wartsila.com/encyclopedia/term/air-cushion-vehicle-(acv)-hovercraft
[3] https://www.programiz.com/python-programming/methods/built-in/zip
[4] https://www.programiz.com/python-programming/methods/built-in/tuple