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

AIM:

In determining the minimum cushion pressure needed to break a given thickness of iceusing an air cushion vehicle, Muller (“Ice Breaking with an Air Cushion Vehicle”) inMathematical 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 aircushion, σ 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 tablewith 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

Solve by using Newton Raphson method and Write a Python code.

Use relaxation factor of 1 and explain your results.

How will you choose the optimum relaxation factor for this case? and what is the optimum relaxation factor for this case? Show with the help of a graph, assume h = 0.6.

Explain your results with respect to the appropriate relaxation factor.

Also, create a small report summarizing your observations.

When solving the problem DO NOT convert to the metric system, use the given system.

SOLUTION:

1. Finding Minimum Pressure Using Relaxation Factor As 1:

CODE:

"""
Finding minimum pressure using relaxation factor as 1 to break ice by air cushion vehicle using Newton Raphson method
"""

import math
import matplotlib.pyplot as plt


#inputs
beta=0.5   #width of ice wedge(no units)
r=40       #size of air cushion(feet)
segma=150   #tensile strength of ice(psi)
segma=150*144   #converting psi into psf
t=[0.6,1.2,1.8,2.4,3.0,3.6,4.2]    #thickness of ice(feet)

#Here
alpha=1   #relaxation factor
tol=1e-6    #tolerance

#simplify the equation
a=1-pow(beta,2)
b=0.4*pow(beta,2)
c=segma/pow(r,2)
d=pow(segma,2)/(3*pow(r,4))
e=pow(segma,3)/(9*pow(r,6))

#defining function
def f(p,h):
	return (a*pow(p,3))+((b*h-c*pow(h,2))*pow(p,2))+(d*pow(h,4)*p)-(e*pow(h,6))

def fprime(p,h):
	return (a*3*pow(p,2))+((b*h-c*pow(h,2))*2*p)+(d*pow(h,4))


#finding minimum pressure to break ice
pressure=[]
iterations=[]

for h in t:
	ct=1
	pg=500
	while(abs(f(pg,h))>tol):
		pg=pg-alpha*(f(pg,h)/fprime(pg,h))
		ct=ct+1
	pressure.append(pg)
	iterations.append(ct)


#printing result
print('Thickness',':','minimum pressure',':','no.of iterations')
for i in range(0,7):
	print(t[i],':',pressure[i],':',iterations[i])

plt.plot(t,pressure)
plt.xlabel('Thickness')
plt.ylabel('Pressure')
plt.show()

OUTPUT:

Pressure VS Thickness

EXPLANATION:

In this question we find the minimum pressure required to break ice with different thickness using relaxation factor as 1.

First of all we get all inputs of beta, segma, r, thickness, relaxation fcator and tolerance.

Then we define the function and its derivative as shown in the coding.

By using for loop, calculate the minimum pressure and no. of iterations takes place.

The next iteration is calculated by the below formula: p_guess=p_guess-alpha*(f(p_guess,h)/fprime(p_guess,h))

Repeat the above procedure and calculate the minimum pressure and no. of iterations takes place for different thickness of ice.

Finally, print the results and make the plot between thickness and pressure.

2. Finding Optimum Relxation Factor Using Thickness as 0.6

CODE:

"""
Finding optimum Relaxation factor using thickness as 0.6
"""

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


#inputs
beta=0.5   #width of ice wedge(no units)
r=40       #size of air cushion(feet)
segma=150   #tensile strength of ice(psi)
segma=150*144   #converting psi into psf
t=0.6    #thickness of ice(feet)

#Here
alpha=np.linspace(0.1,1.8,100)   #relaxation factor
tol=1e-6    #tolerance

#simplify the equation
a=1-pow(beta,2)
b=0.4*pow(beta,2)
c=segma/pow(r,2)
d=pow(segma,2)/(3*pow(r,4))
e=pow(segma,3)/(9*pow(r,6))

#defining function
def f(p):
	h=0.6
	return (a*pow(p,3))+((b*h-c*pow(h,2))*pow(p,2))+(d*pow(h,4)*p)-(e*pow(h,6))

def fprime(p):
	h=0.6
	return (a*3*pow(p,2))+((b*h-c*pow(h,2))*2*p)+(d*pow(h,4))


#finding minimum pressure to break ice
pressure=[]
iterations=[]

for i in alpha:
	ct=1
	pg=500
	while(abs(f(pg))>tol):
		pg=pg-i*(f(pg)/fprime(pg))
		ct=ct+1
	pressure.append(pg)
	iterations.append(ct)


#printing result
print('Relaxation factor',':','minimum pressure',':','no.of iterations')
for j in range(0,100):
	print(alpha[j],':',pressure[j],':',iterations[j])

plt.plot(alpha,iterations)
plt.xlabel('Relaxation factor')
plt.ylabel('No. of iterations')
plt.show()

OUTPUT:

No. of Iterations VS Relaxation Factor:

EXPLANATION:

In this question, we find the optimum relaxation factor for thickness 0.6.

First of all we get all the input values like beta, segma, r, thickness and tolerance.

Then we get 100 relaxation factor values in the range of 0.1 to 1.8 by using linspace command as shown in the coding

By using for loop, we calculate the minimum pressure for every relaxation factor values and thickness as 0.6 as shown in the coding.

Finally, print the results and make plot between relaxation factor and no. of iterations.

Now we see the above result relaxation factor 1.0 has minimum number of iterations as 18.

So, the optimum relaxation factor for thickness 0.6 is 1.0.

OBSERVATION:

From this concept, Newton Raphson method is used to find the exact pressure required to break the different thickness of ice.

Here, the relaxation factor play a major role in finding minimum presssure with less number of iterations.

From the above plots, If you give high relaxation factor then it takes more number of iterations as well as if you give very less factor then it also takes more number of iterations for calculating minimum pressure.

So, for ervey thickness there must be a optimum relaxation factor value.

From the above plots, the exact pressure required is increased with increase in thickness of ice.