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

AIM : To write a program in PYTHON which calculates the minimum pressure required to break the given ice using an air cusion vehicle , by NEWTON - RAPHSON method

INTRODUCTION:

Air cushion vehicles(ACV) are used to break ice in cold places, where ice is formed on the surface of lakes, rivers and on the roads too. Depending on the ice thickness , an ACV can readily break ice in flexure by making a depression and send cracks throgh the ice. Once a crack is made , the path of least resistance for the ACV's air curtain gets under the ice and reduces the supporting water under the ice. The ice then fails under it's own weight. Thus the ice breaks while the vehicle moves forward by applying the pressure radially on the ice.

Various parameters in modelling ice breaking operation include thickness of the ice layer (h), pressure (p), tensile strength of ice(σ),size of air cushion (r) and (β) corelated to the width of the ice wedge.

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

In this program we are finding the minimum pressure required to break the ice for varying thickness, tensile strength air cushion size and β by using Newton-Raphson method which is a iterative process.

import math
import matplotlib.pyplot as plt

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



def fprime(p,h):                         # derivation of a function f'(x)
	return (3*(pow(p,2)*(1-pow(beta,2)))) + (2*p*((0.4*h*pow(beta,2))-(s*pow(h,2)/pow(r,2)))) + ((pow(s,2)*pow(h,4))/(3*pow(r,4)))

# INPUTS 

beta = 0.5                               # width of ice wedge
r = 40                                   # size of air cushion 
s = 150                                  # tensile strength of the ice (SIGMA)
height = [0.6,1.2,1.8,2.4,3.0,3.6,4.2]   # thickness of ice field 
p_guess = 0.1                            # cushion pressure
alpha = 1                                # relaxation factor
t = 1e-8                                 # tolerance
pressure = []
c = 1                                    

for h in height :
	while(abs(fun(p_guess,h))>t):        # absolute value of f(x guess) should be less than tolerance provided 
		p_guess = p_guess - (alpha*(fun(p_guess,h)/fprime(p_guess,h)))
		c = c + 1
	pressure.append(p_guess)
	

for i in range (0,7) :                  # creating a table
	if i == 0 :
		print ("Sl.no Height  Pressure ")

	print(i+1 ,"    ", height[i] , "=", pressure[i])


# GRAPH PLOTS

plt.plot(height,pressure)
plt.xlabel('Height')
plt.ylabel('Pressure ')
plt.xlim([0,5])
plt.ylim([0,0.5])
plt.show()

OUTPUT:

For relaxation faction Alpha = 1 , and for seven iterations time = 20.7s

Hence the minimum pressure required to break the ice is 0.00305547 psi , and needs to be done at a height of 0.6 feet.

Optimum relaxation factor : Lets consider h = 0.6 and then when the program was executed with an alpha value of 3,5,8 etc i.e alpha> 1 the output was not shown for a very long time as its an over relaxed condition. When alpha <= 1 we get the optimum results . The under-relaxation factor for all variables are set to values that are near optimal for the largest possible number of cases. These values are suitable for many problems, but for some particularly nonlinear problems.

GRAPHS :

for alpha= 1