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

AIM: Find minimum pressure of Breaking Ice with Air cushion Vehicle with Newton-Raphson method.

OBJECTIVES:

• Use the newton Raphson method to find out the value of pressure for h = 0.6 ft.
• Find the optimal relaxation factor for this problem with the help of a suitable plot.
• 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.

BREAKING ICE WITH AIR CUSHION VEHICLE

Equation of breaking ice with air cushion vehicle is as follow:

`p^3 (1-β^2 )+(0.4hβ^2-(σh^2)/r^2 ) p^2+(σ^2 h^4)/(3r^4 ) p-((σh^2)/(3r^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.

We have to take β = 0.5, r = 40 feet and σ = 150 pounds per square inch (psi) for the programming.

NEWTON-RAPHSON METHOD

Newton-Raphson method is a root-finding algorithm which produces successively better approximation to the roots (or zeroes) of a real -valued function. The most basic version starts with a single-variable function f defined for a real variable x, the function's derivative f ′, and an initial guess x₀ for a root of f. If the function satisfies sufficient assumptions and the initial guess is close. Then,

`x_1=x_0-f(x_0 )/(f'(x_0 ))`

is a better approximation of the root than x₀. Geometrically, (x₁, 0) is the intersection of the x-axis and the tangent of the graph of f at (x₀, f (x₀)): that is, the improved guess is the unique root of the linear approximation at the initial point. The process is repeated as

`x_(n+1)=x_n-α f(x_n )/(f'(x_n ))` where α relaxation factor

So, we use this equation for finding the minimum pressure.

1.At h=0.6

PROCEDURE

1. We have to import important module like math, scipy, numpy and tabulate.
2. I give initial values of β, σ, r and h values.
3. Now I define the function which have variable p and h as above equation. For simplicity I break this equation in 4 terms.
4. Then I define derivative function of these equation.
5. Now I guess the pressure value as 30.
6. For newton Raphson method we assume the relaxation factor equal to 0.6
7. I define tolerance value for stop the newton Raphson method.
8. I define iteration value equal to 1
9. Then I use while loop to find the minimum pressure vale with condition of absolute value of function is greater than tolerance value.
10. I define the header name for the table.
11. Then give the value to table variable which I want to print in the table and print the table.

PROGRAM

import math
import matplotlib.pyplot as plt 
import numpy as np 
from tabulate import tabulate   #for otput in tabular form we have to install this module

#inputs
beta=0.5             #related to width of the ice wedge
r=40                 #size of the air cushion in feet
sigma=150*144        #tensile strength of the ice in pound per square feet
h=0.6                #thickness of the ice field in feet

#define the function
def f(p,h):
	
	term_1=pow(p,3)*(1-pow(beta,2))
	term_2=(0.4*h*pow(beta,2)-(sigma*pow(h,2)/pow(r,2)))*pow(p,2)
	term_3=p*((pow(sigma,2)*pow(h,4))/(3*pow(r,4)))
	term_4=pow((sigma*pow(h,2))/(3*pow(r,2)),3)

	equ=term_1+term_2+term_3-term_4
	return equ

#define the derivative function
def fprime(p,h):
	term_1=3*pow(p,2)*(1-pow(beta,2))
	term_2=(0.4*h*pow(beta,2)-(sigma*pow(h,2)/pow(r,2)))*2*p
	term_3=((pow(sigma,2)*pow(h,4))/(3*pow(r,4)))

	equ_prime=term_1+term_2+term_3
	return equ_prime

#guess the value of pressure
p_guess=30

#relaxation factor
alpha=0.6

#tolerance value
tol=1e-4

iteration=1

#newton raphson method
while(abs(f(p_guess,h))>tol):
	p_guess=p_guess - alpha*(f(p_guess,h)/fprime(p_guess,h))
	iteration=iteration+1

#for output in table form
headers=['Thickness of ice field','cushion pressure(psf)','No.of iteration','relaxation factor']
table=[[h,p_guess,iteration,alpha]]
print(tabulate(table, headers, tablefmt="pipe",numalign="center"))

OUTPUT

2.finding optimum relaxation factor

PROCEDURE

1. We have to import important module like math, scipy, numpy and tabulate.
2. I give initial values of β, σ, r and h values.
3. Now I define the function which have variable p and h as above equation. For simplicity I break this equation in 4 terms.
4. Then I define derivative function of these equation.
5. For newton Raphson method we assume 20 values of relaxation factor between 0.1 to 1.9.
6. Now I create the 2-empty array for the pressure and no. of iteration.
7. I define tolerance value for stop the newton Raphson method.
8. I use for loop for finding pressure value and no. of iteration for each value of relaxation factor.
9. I define the p_guess and iteration value for over writing after each time while loop run for finding the no. of iteration value required for get minimum pressure from 30 at particular relaxation factor.
10. Then I use while loop in for loop to find the minimum pressure value with condition of absolute value of function is greater than tolerance value.
11. I append the p_guess value and no. of iteration value in empty array.
12. I define the header name for the table.
13. Then give the value to table variable which I want to print in the table for these I am using for loop because variable is in the form of array and print the table.
14. I plot the graph of no. of iteration vs relaxation factor with x and y label and save the graph.

PROGRAM

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

#inputs
beta=0.5             #related to width of the ice wedge
r=40                 #size of the air cushion in feet
sigma=150*144        #tensile strength of the ice in pound per square feet
h=0.6                #thickness of the ice field in feet

#define the function
def f(p,h):
	
	term_1=pow(p,3)*(1-pow(beta,2))
	term_2=(0.4*h*pow(beta,2)-(sigma*pow(h,2)/pow(r,2)))*pow(p,2)
	term_3=p*((pow(sigma,2)*pow(h,4))/(3*pow(r,4)))
	term_4=pow((sigma*pow(h,2))/(3*pow(r,2)),3)

	equ=term_1+term_2+term_3-term_4
	return equ

#define the derivative function
def fprime(p,h):
	term_1=3*pow(p,2)*(1-pow(beta,2))
	term_2=(0.4*h*pow(beta,2)-(sigma*pow(h,2)/pow(r,2)))*2*p
	term_3=((pow(sigma,2)*pow(h,4))/(3*pow(r,4)))

	equ_prime=term_1+term_2+term_3
	return equ_prime

#relaxation factor
alpha=np.linspace(0.1,1.9,20)

#empty array for store pressure value
v_p=[]

#empty array for store num. of iteration
n_i=[]

#tolerance value
tol=1e-4


#for loop for newton raphson method
for i in alpha:
	iteration=1
	p_guess=30
	while(abs(f(p_guess,h))>tol):
		p_guess=p_guess - i*(f(p_guess,h)/fprime(p_guess,h))
		iteration=iteration+1
	v_p.append(p_guess)
	n_i.append(iteration)

#for output in table form
headers=['Thickness of ice field','cushion pressure(psf)','No.of iteration','relaxation factor']
table=[]
for i in range(0,20):
	table_tmp=[h,v_p[i],n_i[i],alpha[i]]
	table.append(table_tmp)
print(tabulate(table, headers, tablefmt="pipe",numalign="center"))

#plotting graph
plt.plot(alpha,n_i,'r-')
plt.ylabel('No. of iteration')
plt.xlabel('Relaxation Factor')
plt.savefig('optimum relaxation factor.png')
plt.show()

OUTPUT

GRAPH OF NO. OF ITERATION VS RELAXATION FACTOR

3. P FOR h = [0.6,1.2,1.8,2.4,3,3.6,4.2]

PROCEDURE

1. We have to import important module like math, scipy, numpy and tabulate.
2. I give initial values of β, σ, r and h values.
3. Now I define the function which have variable p and h as above equation. For simplicity I break this equation in 4 terms.
4. Then I define derivative function of these equation.
5. For newton Raphson method we assume relaxation factor equal to 0.6.
6. Now I create the 2-empty array for the pressure and no. of iteration.
7. I define tolerance value for stop the newton Raphson method.
8. I use for loop for finding pressure value and no. of iteration for each value of h.
9. I define the p_guess and iteration value for over writing after each time while loop run for finding the no. of iteration value required for get minimum pressure from 30 at particular h.
10. Then I use while loop in for loop to find the minimum pressure value with condition of absolute value of function is greater than tolerance value.
11. I append the p_guess value and no. of iteration value in empty array.
12. I define the header name for the table.
13. Then give the value to table variable which I want to print in the table for these I am using for loop because variable is in the form of array and print the table.
14. I plot the graph of no. of iteration vs relaxation factor with x and y label and save the graph.

PROGRAM

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

#inputs
beta=0.5                                  #related to width of the ice wedge
r=40                                      #size of the air cushion in feet
sigma=150*144                             #tensile strength of the ice in pound per square feet
h=[0.6,1.2,1.8,2.4,3,3.6,4.2]             #thickness of the ice field in feet

#define the function
def f(p,h):
	
	term_1=pow(p,3)*(1-pow(beta,2))
	term_2=(0.4*h*pow(beta,2)-(sigma*pow(h,2)/pow(r,2)))*pow(p,2)
	term_3=p*((pow(sigma,2)*pow(h,4))/(3*pow(r,4)))
	term_4=pow((sigma*pow(h,2))/(3*pow(r,2)),3)

	equ=term_1+term_2+term_3-term_4
	return equ

#define the derivative function
def fprime(p,h):
	term_1=3*pow(p,2)*(1-pow(beta,2))
	term_2=(0.4*h*pow(beta,2)-(sigma*pow(h,2)/pow(r,2)))*2*p
	term_3=((pow(sigma,2)*pow(h,4))/(3*pow(r,4)))

	equ_prime=term_1+term_2+term_3
	return equ_prime

#relaxation factor
alpha=0.6

#empty array for store pressure value
v_p=[]

#empty array for store num. of iteration
n_i=[]

#tolerance value
tol=1e-4

#for loop for newton raphson method
for i in h:
	iteration=1
	p_guess=30            #guess the value of pressure
	while(abs(f(p_guess,i))>tol):
		p_guess=p_guess - alpha*(f(p_guess,i)/fprime(p_guess,i))
		iteration=iteration+1
	v_p.append(p_guess)
	n_i.append(iteration)

#for output in table form
headers=['Thickness of ice field','cushion pressure(psf)','No.of iteration','relaxation factor']
table=[]
for i in range(0,7):
	table_tmp=[h[i],v_p[i],n_i[i],alpha]
	table.append(table_tmp)
print(tabulate(table, headers, tablefmt="pipe",numalign="center"))

#plotting graph
plt.plot(h,v_p,'ro-')
plt.xlabel('Thickness of ice field in feet')
plt.ylabel('Cushion Pressure(PSF)')
plt.savefig('P vs h.png')
plt.show()

OUTPUT

GRAPH OF P vs h

CONCLUSION

• As we increase the value of relaxation factor upto 1, we get results very rapidly so require less no.of iteration. But for further increment of relaxation factor generate over-shooting which require more time. For this problem optimum relaxation factor is 1.
• In this problem relaxation factor value must be less than 2 because at relaxation factor 2 while loop runs infinite times.
• As we increase h, pressure value increase rapidly.