Finding the minimum pressure for breaking ice by 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 .

THEORY : To determine 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 below equation f(x) =

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

To calculate the minimum pressure we do an approximation method , where certain number of iteration are carried out and the closest result is considered. By Newton - Raphson's techinque we get the following equation

p_new = p_guess - (`alpha` * (f(x) / f'(x)))

`here f'(x) = 3p^2[1-beta^2] + (0.4h*beta^2)2p - (sigma*h^2)/r^2 + (sigma^2*h^4)/3r^`

Relaxation Factor (`alpha`) : In numerical mathematics , relaxation methods are iterative methods for solving system of equations , including nonlinear systems. Relaxation factors play a vital role in determining the rate of convergence for simulation. When the relaxation factor is less than one it is said to be under relaxed and it is called over relaxed when the relaxation factor is more than one. Relaxation factor is always a constant and it is multiplied to the algebraic equations to alter the path of the iteration.

PYTHON CODE :

""" PROGRAM TO FIND THE MINIMUM PRESSURE TO BREAK THE ICE WITH AIR CUSHION VEHICLE ,WITH NEWTON RAPHSON METHOD  """

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()

For relaxation faction (`alpha`) = 1 , and for the seven iterations it took 20.7s to finish .

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.

for `alpha`= 3 , no output .

GRAPHS :

for `alpha`= 1