Breaking Ice with Air cushion Vehicle using the Newton-Raphson method in Python

Objective:

1) Writing a python script to solve for the minimum value of pressure using the newton raphson method.

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

2) Using the newton raphson method to find out the value of pressure for h = 0.6 ft.

3) Finding the optimal relaxation factor for this problem with the help of a suitable plot.

4) Tabulating the results of p for h = [0.6,1.2,1.8,2.4,3,3.6,4.2] assuming a suitable relaxation factor.

Theory:

Newton Raphson Method:

The Newton-Raphson method (also known as Newton's 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. Suppose you need to find the root of a continuous, differentiable function f(x), and you know the root you are looking for is near the point (`x=x_0`). Then Newton's method tells us that a better approximation for the root is

`x _1 =x _0− ( f ′ (x_0))/ f(x_0)`

This process may be repeated as many times as necessary to get the desired accuracy. In general, for any x-value `x_n`​, the next value is given by

`x_(n+1) ​ =x_n −(f ′(x_n))/ f(x_n)`

Ice Breaking with an Air Cushion Vehicle:

A hovercraft, also known as an air-cushion vehicle or ACV, is an amphibious craft capable of travelling over land, water, mud, ice, and other surfaces. Hovercraft use blowers to produce a large volume of air below the hull, or air cushion, that is slightly above atmospheric pressure. The pressure difference between the higher pressure air below the hull and lower pressure ambient air above it produces lift, which causes the hull to float above the running surface. For stability reasons, the air is typically blown through slots or holes around the outside of a disk- or oval-shaped platform, giving most hovercraft a characteristic rounded-rectangle shape. Typically this cushion is contained within a flexible "skirt", which allows the vehicle to travel over small obstructions without damage.

The Formula for ice breaking is as follows

`p^3 ( 1 − β^2 ) + ( 0.4*h*β^2 − (σ*h^2)/ (r^2) ) p^2 + ((σ^2*h^4)/ (3*r^4))*p − ((σ*h^2)/(3*r^2))^3 = 0`

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

Procedure:

Full Code:

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

sigma = 150*144
beta = 0.5
r = 40



def f(p,h):

	f = (1-beta**2)*p**3 + (0.4*h*beta**2 - sigma*h**2/r**2)*p**2 + (sigma**2*h**4)*p/(3*r**4) - ((sigma*h**2)/(3*r**2))**3
	return f

def fprime(p,h):

	fprime = (3*p**2*(1-beta**2) + 2*p*(0.4*h*beta**2 - sigma*h**2/r**2) + sigma**2*h**4/(3*r**4))

	return fprime

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


p_guess = 10
alpha = 1.0
tol = 1e-4
iter = 1
iterations = []
P = []

for h in h_range:

	while(abs(f(p_guess,h)) > tol):
		p_guess = p_guess - alpha*(f(p_guess,h)/fprime(p_guess,h))
		iter = iter + 1

	iterations.append(iter)
	P.append(p_guess)

	
print(P)
print(iterations)

print('-----------------------------------------------------------------')
print('Height(ft)','|','iterations','|','Cushion pressure(psft)')
print('-----------------------------------------------------------------')
for i in range(0,7):
    print(' ',h_range[i],'              ',iterations[i],'         ',P[i])
print('-----------------------------------------------------------------')

plt.plot(h_range,P)
plt.xlabel('Thickness of ice (ft)')
plt.ylabel(' Minimum cushion pressure (psft)')
plt.title('Cushion pressure V/S ice thickness')
plt.grid('on')
plt.savefig('plot1')
plt.show()


print('-----------------------------------------------------------------')
print('CASE 2: find the optimum relaxation factor for ice thickness = 0.6 ft')
print('-----------------------------------------------------------------')

alpha_range = np.arange(0.02,2.0,0.02)
iterations =[]
i_min =10000
alpha_best = []

for alpha in alpha_range:
    i=1
    p_guess=10
    while (abs(f(p_guess,0.6))>tol):
        p_guess = p_guess-alpha*f(p_guess,0.6)/fprime(p_guess,0.6)
        i = i + 1
    if (i_min > i):
        i_min = i
        alpha_opt = alpha
    iterations.append(i)

#print(alpha_opt)
#print(iterations)
print('minimum number of iterations = ',i_min, 'at alpha = ','{:.2f}'.format(alpha_opt))

#printing the table
print('-----------------------------------------------------------------')
print('Relaxation factor(alpha) ','|','iterations')
print('-----------------------------------------------------------------')
for i in range(0,99):
    print('    ','{:.2f}'.format(alpha_range[i]),'                       ',iterations[i])
print('-----------------------------------------------------------------')
plt.plot(alpha_range,iterations)
plt.plot(alpha_opt,i_min,marker = '*',color = 'red',label = 'point corresponding to minimum number of iterations')
plt.legend()
plt.xlabel('relaxation factor(alpha)')
plt.ylabel('iterations to converge')
plt.grid('on')
plt.savefig('optimum alpha')
plt.show()

Stepwise Explanation:

1) We first, import the libraries like math library, plotting library and numpy library for making necessary calculation and plotting the graphs.

2) We then define the functions for the formula given for ice breaking i.e.

p3(1−β2)+(0.4⋅h⋅β2−σ⋅h2r2)p2+(σ2⋅h43⋅r4)⋅p−(σ⋅h23⋅r2)3=0p3(1−β2)+(0.4⋅h⋅β2−σ⋅h2r2)p2+(σ2⋅h43⋅r4)⋅p−(σ⋅h23⋅r2)3=0

3) We then use the while loop for perfoming iterations for pressure of the cushion vehicle according to the thickness of the ice.

4) Then, we plot the graph for different height of the ice as provided in the aim.

5) We then, proceed to find the optimal relaxation factor again by performing iterations.

6) We then, plot the graph for optimal relaxation factor and also mark the optimal alpha value according to the iterations.

7) We also print the table for relaxation factor and iterations.

Case 1:

In case-1 we get the value of pressure for height of 0.6 ft as 4.239 psft as shown in the table in the image given below.

We also get the graph for Cushion Pressure vs Thickness Plot.

Case 2:

In case 2 we get the optimal relaxation factor for ice thickness = 0.6 ft. We also plot the graph between iterations vs relaxation factor (alpha)

Conclusion:

1) We have written a python script to solve for the minimum value of pressure using the newton raphson method.

2) We have find out the value of pressure for h = 0.6 ft i.e. 4.239 psft.

3) We have Tabulated the results of p for h = [0.6,1.2,1.8,2.4,3,3.6,4.2] assuming a suitable relaxation factor.

4) We have plotted graphs for Cushion Pressure vs Thickness Plot and iterations vs relaxation factor (alpha) plot.