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

Aim: Determining minimum pressure to break ice using an air cushion ice breaker - Newton Raphson method

Objective:

1. Use the newton raphson method to find out the value of pressure for h = 0.6.

2. Find the optimal relaxation factor for this problem with the help of a suitable plot.

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

Background:

Air Cushion Vehicle:

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

f=p3(1−β2)+p2(0.4hβ2−σh2r2)+p(σ2h43r4)−(σh23r2)3=0

∂f∂p=3p2(1−β2)+2p(0.4hβ2−σh2r2)+(σ2h43r4)

Numerical solution of equations:

"Numerical methods are often, of a repetitive nature. These consist in repeated execution of the same process where at each step the results of the preceding step is used. This is known as iteration process and is repeated till the result is obtained to a desired degree of accuracy".[1]

Newton Raphson Method: [1]

Let approximate root be x0

for a function f(x)=0

If x1=x0+h

is the exact solution then f(x1)=0

⇒f(x0+h)=0

Expading using Taylor Series and neglectic ^2 and higher order terms we get:

f(x0)+hf'(x0)=0

⇒h=−f(x0)f'(x0)

⇒x1=x0−f(x0)f'(x0)

Newton-Raphson formula

This is also used with a relaxation factor based on the function to be optimized

x1=x0−rf⋅f(x0)f'(x0)

Code:

# Breaking Ice with Air cushion Vehicle - Find minimum pressure with Newton-Raphson method
# By : Anubrat Padhee

# Importing libraries
import math
import matplotlib.pyplot as plt
import numpy as np

#p denotes the cushion pressure 
#h the thickness of the ice field 
h=[0.6,1.2,1.8,2.4,3.0,3.6,4.2]
#r the size of the air cushion unit=feet
r=40
#σ the tensile strength of the ice unit=psf (pound per square foot)
sigma=150*144
#β is related to the width of the ice wedge
beta=0.5

#function
def f(p,h):
	return (pow(p,3)*(1-pow(beta,2)))+(((0.4*h*pow(beta,2))-(sigma*pow(h,2)/pow(r,2)))*pow(p,2))+(((pow(sigma,2)*pow(h,4))/(3*pow(r,4)))*p)-(pow((sigma*pow(h,2)/(3*pow(r,2))),3))

#derivative of function
def fprime(p,h):
	return (3*(pow(p,2)*(1-pow(beta,2))))+(2*p*((0.4*h*pow(beta,2))-(sigma*pow(h,2)/pow(r,2))))+((pow(sigma,2)*pow(h,4))/(3*pow(r,4)))

# Performing newton iteration
p_guess=10

pressure=[]

# Tolerance
tol=1e-4

# Relaxation factor
alpha=0.5

for i in h:
	while(abs(f(p_guess,i))>tol):
		p_guess=p_guess-alpha*(f(p_guess,i)/fprime(p_guess,i))
	pressure.append(p_guess)

print("The minimum value of p for the value of pressure for h = 0.6 at relaxation factor 0.5 is : ",pressure[0])
print("Tabulation of the results of p for h = [0.6,1.2,1.8,2.4,3,3.6,4.2] for a  relaxation factor of 0.5.")
print("height : pressure")
for i in range(len(h)):
	print(h[i],"   :",pressure[i])

# To find the most suitable relaxation factor
h1=0.6
alpha_all=np.linspace(0.1,1.5,50)

#iterations
iterations=[]

for alpha in alpha_all:
	iter=1
	p_guess=10
	while(abs(f(p_guess,h1))>tol):
		p_guess=p_guess-alpha*(f(p_guess,h1)/fprime(p_guess,h1))
		iter=iter+1	
	iterations.append(iter)		

# Plotting
plt.plot(alpha_all,iterations)
plt.xlabel("Relaxation factor")
plt.ylabel("iterations")
plt.show()	

Explanation of the code:

• First the important libraries were imported.
• An array of different thickness of the ice field were taken.
• The size of the air cushion, the tensile strength of the ice unit, β is related to the width of the ice wedge variable were taken.
• Functions were created for the minimum cushion pressure needed to break a given thickness of ice using an air and the function derivative.
• Performing newton iterations taking p_guess=10, taking tolerance value as 1e-4 and taking relaxation factor as 0.5.
• Using newton raphson method p_guess value was appeneded into an empty pressure array for different values of h.
• Then the minimum pressure was tabulated for different values for h.
• Then, taking h=0.6 number of iterations were appended into an empty array of iterations for different values for relaxation factor.
• Then, we plot the graph between relaxation factor and iterations.
• So, now seeing the graph we can see that number of iterations is minimum for relaxation factor nearly equal to one.

Output:

Conclusions:

• The minimum value of pressure for different value of thickness of ice shield was tabulated.
• So, now seeing the graph we can see that number of iterations is minimum for relaxation factor nearly equal to one.
• So, we can say optimal relaxation factor for h=0.6 is nearly equal to one.