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

Objective: To determine the minimum cushion pressure needed to break a given thickness of ice using

an air cushion vehicle, Muller (“Ice Breaking with an Air Cushion Vehicle”) using the following equation

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).

Given that:

1. To construct a table with h in the first column, and p in the second column for h = 0.6, 1.2, 1.8, 2.4, 3.0, 3.6and 4.2 feet. (relaxation factor = 1)
   
   
2. To find the optimum relaxation factor by plotting a graph. (h = 0.6)
   
   
3. Explaining the results with the help of appropriate relaxation factor.

Air Cushion Vehicle:

In determining the minimum cushion pressure needed to break a given thickness of ice using an air

cushion vehicle, Muller (“Ice Breaking with an Air Cushion Vehicle”)

The equations:

Iteration: Numerical methods are often, of a repetitive nature. These consist in repeated execution of

thesame process where at each step the results of the preceding step is used. This is known as

iteration process and is repeated till the result is obtained to a desired degree of accuracy.

Relaxation factor: - it represents the value of reduction of the current load step for the restart

providing the solution fails to converge. A new value of the assumed load step is found from the

express.

Newton Raphson method:

The Newton-Raphson method uses an iterative process to approach one root of a function. The specific

root that the process locates depends on the initial, arbitrarily chosen x-value.

- Depending on the conditions under which we are attempting to solve this equation, several of the variables may be changing. So, it may be necessary to use partial derivatives.

- For the purposes of this example, we are assuming that height is the only thing changing, and that

these values are all functions of time. This avoids the use of a partial derivative, we simply differentiate

all variables with respect to pressure, as shown above.

Procedure:

1. Find f'(x) and define NR equation.

2. Guess an initial value of x for the first iteration.

3. substitute x in NR equation and calculate x*.

4. If |x-x*|>tolerance, stop iteration and outpur the root:x*

5. If number of iterations reaches an assumed maximum, stop.

6. Else let x=x* and repeat the steps from 3 until a condition in 4 or 5 is satisfied.

Coding:

import matplotlib.pyplot as plt 
import math
import numpy as np
from prettytable import PrettyTable

#input values
beta = 0.5              #width of ice berg
r = 40                  #size of air cushion
sigma = 150*144         #tensile strengh
h = [0.6, 1.2, 1.8, 2.4, 3, 3, 3.6, 4.2]         #feet
p = [] 

def f(p,h):                                                                                            
	x1 = pow(p,3)*(1-pow(beta,2))
	x2 = (0.4*h*pow(beta,2) -(sigma*pow(h,2)/pow(r,2)))*p**2
	x3 = pow(sigma,2)*pow(h,4)*p/(3*pow(r,4))
	x4 = pow((sigma*h*h/(3*r*r)),3)
	f1= x1+x2+x3-x4
	return f1


def fp(p,h):
	y1 = 3*pow(p,2)*(1-beta*beta) 
	y2 = 2*p*((0.4*h*beta*beta)-(sigma*h*h/(r*r)))
	y3 = pow(sigma,2)*pow(h,4)/(3*pow(r,4))
	fp = y1+y2+y3
	return fp

#performing newton iteration

x_guess = 10
alpha = 1
tol = 1e-4
it=[]

for i in h:
	iter = 1	
	while(abs(f(x_guess,i))>tol): 
		x_guess = x_guess-alpha*(f(x_guess, i)/fp(x_guess, i))
		iter = iter+1
		print(x_guess, iter, h) 
	p.append(x_guess)
	it.append(iter)





table = PrettyTable(['iteration', 'Thickness', 'Pressure'])

for y in range (0,8):
	table.add_row([it[y],h[y],p[y]])
print(table)

alpha = np.linspace(0.1,1.9,100)

iterations = []

tol = 1e-4

for i in alpha:
	iter=1	
	x_guess = 10				
	while(abs(f(x_guess,h[0]))>tol):
		x_guess = x_guess-i*(f(x_guess,h[0])/fp(x_guess,h[0]))
		iter = iter+1
	print(x_guess,iter,h)
	print(i)
	iterations.append(iter)
print(iterations)

print('the number of iterations('+str(min(iterations))+') is least at :')
for i in range(len(iterations)):

	if iterations [i]==min(iterations):
		print(alpha[i])


plt.plot(h,p, color='tab:blue', marker='o', markersize=5)
plt.xlabel('height')
plt.ylabel('pressure')
plt.title('Pressure vs Thickness')
plt.savefig('pressure vs Thickness plot')
plt.show()


plt.plot(alpha,iterations, '.-')
plt.xlabel('aplha')
plt.ylabel('iterations')
plt.title('Relaxation factor')
plt.savefig('plot of Relaxation factor')
plt.show()

output:

Ploting:

Limitations of Newton Raphson method:

1. If f(x) = 0, this method will not work.

2. If f''(x) = 0, this causes the solution to diverge.

3. Sometimes get trapped with local minima and maxima.

Conclusion:

The plots and results show that by using a suitable relaxation factor we can obtain results with as much accuracy, in a shorter processing time. We can also observe that the RF is not

directly related to the number of iteration. In fact the relation is highly dependent on the function being optimized.