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

Objective: In determining the minimum cushion pressure needed to break a given thickness of ice
using an air cushion vehicle, Muller derived the equation

P³(1-β²)+(0.4hβ²-σh²/r²)p²+ σ²h⁴p/(3r⁴)+( σh²/3r²)³ =0

where;

β = 0.5, r = 40 feet and σ = 150 pounds per square inch (psi)

Also,h = 0.6, 1.2, 1.8, 2.4, 3.0, 3.6 and 4.2 feet.

Concept: This equation is to be resolved using Newton-Raphson method i.e.

`dy/dx= (f(xguess)-f(xnew))/(xguess-xnew)`

or, xnew=xguess- f(xguess)/f’(xguess)

but for more accuracy a relaxation factor, i.e. `α`

i.e. xnew=xguess-`α`*f(xguess)/f’(xguess)

Therefore, by using the equation we can determine the value of minimum cushion pressure at different thickness of ice.

Methodology: Let's start with declaration of two functions i.e f & d_f,

1. The function f takes pressure(p) and thickness(h) as input and return the value above equation.

2. The function d_f returns the derivative of the above equation.

3. In main function, the first task is to define all the parameters and use a nested a loop to find the value of pressure (p) until the return of f becomes less than tol i.e. 1e-4. and this pressure is to be calculated for all the thickness values and the output is to be shown.

4.Next step will be to calculate the optimum relaxation factor, the factor which finds the solution in the least number of iterations. So, to calculate it we declare a relaxation factor as 0.15 and has calculated no. of iterations for each factor. After calculation, we plot a graph of Iterations vs Relaxation factor.

5. The last step would be to plot the graph between the Cuishion pressure vs thickness of ice for the relaxation factor of 1.1.

Output:

1.

2.

3.

4.

Graphs:

1.Relaxation Factor vs No. of iterations

2. Cushion pressure vs thickness of ice field

Python Program:

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



b= 0.5;
r= 40;
al= 150*144;


# FUNCTION-1: To define the function of minimum cushion pressure 
def f(p,h):

	global b,r,al
	t1= pow(p,3)*(1-pow(b,2))
	t2= (0.4*h*pow(b,2)-al*pow(h,2)/pow(r,2))*pow(p,2)
	t3= pow(al,2)*pow(h,4)/3/pow(r,4)*p
	t4= pow(al*pow(h,2)/3/pow(r,2),3)
	term= t1+t2+t3-t4
	return term

# FUNCTION-1: The derivative of the function of minimum cushion pressure 
def d_f(p, h):
	global b,r,al
	dt1= 3*pow(p,2)*(1-pow(b,2))
	dt2= (0.4*h*pow(b,2)-al*pow(h,2)/pow(r,2))*2*p
	dt3= pow(al,2)*pow(h,4)/3/pow(r,4)
	d_term= dt1+dt2+dt3
	return d_term

# Let's come to the main function now.

h= np.linspace(0.6,4.2,7)
tol= 1e-4
p= 100 #psi
itr=1
print("      Pressure        | Thickness",'\n')
# Using the Relaxation factor AS 1
for i in h:

	while(abs(f(p,i)) > tol):
		p=p-f(p,i)/d_f(p,i)
		itr= itr+1
	print(p,'psf |',i,'\n')

# Now lets find the optimum relaxation factor.
rel=[]
ITR=[]
rel1=0.15
itr=1
p= 100 #psi
print("Relaxation Factor  | Iterations",'\n')

while (rel1<2):
	rel.append(rel1)
	while(abs(f(p,0.6)) > tol):
		p=p-rel1*f(p,0.6)/d_f(p,0.6)
		itr= itr+1
	ITR.append(itr)

	print(rel1,'|',itr,'\n')
	rel1=rel1+0.05
	p=100
	itr=1
min= np.amin(ITR)
plt.xlabel('Relaxation factor')
plt.ylabel('Iterations')
plt.plot(rel,ITR,'-mo',min)
plt.show()

#Using the Optimum Relaxation factor AS 1.1	shown in graph
Pre=[]
p=100

for j in h:
	while(abs(f(p,j)) > tol):
		p=p-1.1*f(p,j)/d_f(p,j)
	Pre.append(p)


plt.ylabel('Cushion Pressure- psf')
plt.xlabel('thickness of ice field-feet')
plt.plot(h,Pre)
plt.show()