Project

BREAKING ICE WITH AIR CUSHION VEHICLE-MINIMUM PRESSURE USING NEWTON RAPHSON METHOD

PPraveen Reddy11 Sept 2019

OBJECTIVE:

The main objective of this challenge is to find the minimum pressure using newton-raphson method and script a code for choosing minimum relaxation factor according to the given equation

FORMULA:

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”) in
Mathematical Modeling: Classroom Notes in Applied Mathematics, SIAM 1987) derived
the 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). 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.6
and 4.2 feet

CODE FOR FINDING MINIMUM PRESSURE BY TAKING RELAXATION FACTOR AS 1 AND h=0.6:

# newton raphson method
import matplotlib.pyplot as plt
import math 
#h = 0.6,1.2,1.8,2.4,3.0,3.6,4.2
h =0.6
b =0.5
r =40
s = 150*144
p=[]
def f(p):
	term1=pow(p,3)*(1-pow(b,2))
	term2 = ((0.4*h*b**2)-(s*h**2/r**2))*p**2
	term3 = s**2*pow(h,4)*p/(3*pow(r,4))
	term4 = pow((s*h**2/(3*r**2)),3)
	F = term1+term2+term3-term4
	return F
def fprime(p):
	term1 = 3*p**2*(1-b**2)
	term2 = 2*p*((0.4*h*b**2)-(s*h**2/r**2))
	term3 =  s**2*pow(h,4)/(3*pow(r,4))
	Fp = term1+term2+term3
	return Fp
x_p =10
tol =1e-4
alpha=1
i=1
while(abs(f(x_p))>tol):
	x_p = x_p - (alpha*f(x_p)/fprime(x_p))
	i=i+1
print(i,(x_p))

COMMANDS USED IN THE CODE:

matplotlib -> matplotlib is an amazing visualization library in python for 2D plots of arrays.This command is used to plot graphs using the code

math -> This library is uesd when we use the trigonometric functions

def () -> This keyword is used to define function along with all suitable parameters related to the problem

return() -> The return statement is used when a function is ready to return a value to its caller.

while() -> The while loop in Python is used to iterate over a block of code as long as the test expression (condition) is true

print() -> This statement is used to print the required output

OUTPUT:

8 4.239063886378783
[Finished in 1.3s]

The minimum pressure value for the above equation is 4.239 and it is obtained after 8 iterations using Newton-raphson method

CODE FOR PRINTING PRESSURES,h-VALUES AND ITERATIONS IN A TABLE FORMAT:

import matplotlib.pyplot as plt
import math 
from tabulate import tabulate
import numpy as np
h = [0.6,1.2,1.8,2.4,3.0,3.6,4.2]
b =0.5
r =40
s = 150*144
def f(p,h):
	term1=pow(p,3)*(1-pow(b,2))
	term2 = ((0.4*h*b**2)-(s*h**2/r**2))*p**2
	term3 = s**2*pow(h,4)*p/(3*pow(r,4))
	term4 = pow((s*h**2/(3*r**2)),3)
	F = term1+term2+term3-term4
	return F
def fprime(p,h):
	term1 = 3*p**2*(1-b**2)
	term2 = 2*p*((0.4*h*b**2)-(s*h**2/r**2))
	term3 =  s**2*pow(h,4)/(3*pow(r,4))
	Fp = term1+term2+term3
	return Fp
x_p =10
tol =1e-4
alpha=1
headers = [\"iteration\",\"h-value\",\"p-vaue\"]

it =[]
P=[]

for a in h:
	i=1
	while(abs(f(x_p,a))>tol):
		x_p = x_p - (alpha*f(x_p,a)/fprime(x_p,a))
		i=i+1
	it.append(i)
	P.append(x_p)

	#print(\"%-15s %-15s %-15s\"%(i,a,x_p))
	#mydata = print(\"%-15s %-15s %-15s\"%(i,a,x_p))
mydata =tuple(zip(it,h,P))

print(tabulate(mydata,headers=headers,tablefmt=\"grid\"))

COMMANDS USED IN THE ABOVE CODE:

matplotlib -> matplotlib is an amazing visualization library in python for 2D plots of arrays.This command is used to plot graphs using the code

math -> This library is uesd when we use the trigonometric functions

def () -> This keyword is used to define function along with all suitable parameters related to the problem

return() -> The return statement is used when a function is ready to return a value to its caller.

while() -> The while loop in Python is used to iterate over a block of code as long as the test expression (condition) is true

print() -> This statement is used to print the required output

tabulate -> This module is used for drawing tables in different form

numpy -> This module is used for assigning n vaues in between two values with same common difference

tuple() -> This command is used to arrange the values in ordered pairs and shown as different rows in array

append() -> This function is used to append the calculated values in an empty array

To install tabulate module go to command prompt window and type \"python -m pip install tabulate\"

OUTPUT:

+-------------+-----------+-----------+
|   iteration |   h-value |    p-vaue |
+=============+===========+===========+
|           8 |       0.6 |   4.23906 |
+-------------+-----------+-----------+
|           7 |       1.2 |  17.2367  |
+-------------+-----------+-----------+
|          24 |       1.8 |  38.9901  |
+-------------+-----------+-----------+
|          18 |       2.4 |  69.4993  |
+-------------+-----------+-----------+
|          40 |       3   | 108.764   |
+-------------+-----------+-----------+
|          21 |       3.6 | 156.785   |
+-------------+-----------+-----------+
|          60 |       4.2 | 213.562   |
+-------------+-----------+-----------+
[Finished in 1.0s]

This is the table related to Iterations ,pressure values and h-values

CODE FOR FINDING MINIMUM RELAXATION FACTOR USING GRAPH:

import matplotlib.pyplot as plt
import math 
from tabulate import tabulate
import numpy as np
h = [0.6,1.2,1.8,2.4,3.0,3.6,4.2]
b =0.5
r =40
s = 150*144
def f(p,h):
	term1=pow(p,3)*(1-pow(b,2))
	term2 = ((0.4*h*b**2)-(s*h**2/r**2))*p**2
	term3 = s**2*pow(h,4)*p/(3*pow(r,4))
	term4 = pow((s*h**2/(3*r**2)),3)
	F = term1+term2+term3-term4
	return F
def fprime(p,h):
	term1 = 3*p**2*(1-b**2)
	term2 = 2*p*((0.4*h*b**2)-(s*h**2/r**2))
	term3 =  s**2*pow(h,4)/(3*pow(r,4))
	Fp = term1+term2+term3
	return Fp
x_p =10
tol =1e-4
alpha=1
headers = [\"iteration\",\"h-value\",\"p-vaue\"]

it =[]
P=[]

for a in h:
	i=1
	while(abs(f(x_p,a))>tol):
		x_p = x_p - (alpha*f(x_p,a)/fprime(x_p,a))
		i=i+1
	it.append(i)
	P.append(x_p)

	#print(\"%-15s %-15s %-15s\"%(i,a,x_p))
	#mydata = print(\"%-15s %-15s %-15s\"%(i,a,x_p))
mydata =tuple(zip(it,h,P))

print(tabulate(mydata,headers=headers,tablefmt=\"grid\"))
Alpha =np.linspace(0.01,1.99,200)
J=[]
CT =[]
L=[]
for j in Alpha:
	x_p=10
	ct=1
	while(abs(f(x_p,0.6))>tol):
		x_p = x_p - (j*f(x_p,0.6)/fprime(x_p,0.6))
		ct=ct+1
	J.append(j)
	CT.append(ct)
	L.append(ct)
W=tuple(zip(L,Alpha))
w=sorted(W)

mydata1 =tuple(zip(J,CT))
headers1 =[\"relaxationfctor\",\"iteration\"]
print(tabulate(mydata1,headers=headers1,tablefmt=\"grid\"))
print(w[0])
plt.figure()
plt.xlabel(\'Thickness\')
plt.ylabel(\'Pressure\')
plt.plot(h,P)
plt.show()
plt.title(\'minimum relaxation factor is 1.0547 at 7th iteration \')
plt.xlabel(\'Relaxation factor\')
plt.ylabel(\'Iteration\')
plt.plot(J,CT)
plt.show()