Analysing ice breaking pressure using a air cushion vehicle by the use of Newton-Raphson method

Objective:

Euler equation:

`(1-beta^2)p^3 + (0.4hbeta^2-(sigmah^2)/r^2)p^2 + ((sigma^2h^4)/(3r^4))p - ((sigmah^2)/(3r^2))^3 =0`            -Equation 1

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.

1. Write a program to calculate the minimum value of pressure required to break the ice field using Newton-Raphson method
2. Write a program to plot the no of iterations vs relaxation factor graph
3. Write a program to calculate the minimum value of pressures corresponding to a range of thickness of ice field

Problem Statement:

1. Usage of the Newton-Raphson method to give the minimum value of pressure
2. Finding the optimal relaxation factor for this problem with the help of a suitable plot
3. Usage of the tabulate module to make a table of pressure, thickness and iterations and see the effect of thickness on pressure

Solution Procedure:

For Part a:

A python program has been made with the values of thickness of the ice field, size of the air cushion, tensile strength of the ice, a factor β related to the width of the ice wedge, relaxation factor of 1 and each coefficient of P in equation 1 predefined and keeping them as constant.

We then defined two functions for f(x) and f’(x) as required in the Newton-Raphson equation

Then, using a while loop, we can calculate the value of f(x) for each iteration of pressure value starting with an initial guess value of 10 until we reach a value of P where f(x)=0

This gives us the minimum value of pressure required to break the ice field

The detailed code is given below:

import math
import matplotlib.pyplot as plt

#factor related to the width of ice wedge
beta=0.5
#thickness of ice field
h=0.6
#tensile strength of the ice in psf (1 psf = 144 psi)
sigma=21600
#size of the air cushion
r=40

#breaking down the calculations of equation 1 into smaller terms
a=1-pow(beta,2)
b=(0.4*h*pow(beta,2))-(sigma*pow(h,2)/pow(r,2))
c=(pow(sigma,2)*pow(h,4))/(3*pow(r,4))
d=pow((sigma*pow(h,2))/(3*pow(r,2)),3)
#Tolerance for f(x)
tol=1e-4
#Initial guess for the minimum value of pressure
pguess=10
#defining a iteration variable
itr=0
#Relaxation factor
alpha=1

#function calculating f(x)
def f1(pguess,a,b,c,d):
	return (a*pow(pguess,3))+(b*pow(pguess,2))+(c*pguess)-d

#function calculating f'(x)
def f2(pguess,a,b,c):
	return (3*a*pow(pguess,2))+(2*b*pguess)+c

#condition to run the loop until we reach the required value of f(x)
while(abs(f1(pguess,a,b,c,d))>tol):

	#newton-raphson formula
	pguess=pguess-(alpha*(f1(pguess,a,b,c,d)/f2(pguess,a,b,c)))
	itr=itr+1

#printing minimum value of pressure and corresponding iterations
print("Min Pressure:",pguess)
print("Iterations:",itr)

	

For Part b:

The process is the same as in part one only with the additional use of a for loop to calculate the values of minimum pressure corresponding to each value of relaxation factor. And later, used the logic that the best optimised relaxation factor would be the one that gives the minimum value of pressure in the required tolerance zone in the minimum no of iterations.

For that we have plotted the iterations vs relaxation value graph.

Detailed code given below:

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

#factor related to the width of ice wedge
beta=0.5
#thickness of ice field
h=0.6
#tensile strength of the ice (1 psf = 144 psi)
sigma=21600
#size of the air cushion
r=40

#breaking down the calculations of equation 1 into smaller terms
a=1-pow(beta,2)
b=(0.4*h*pow(beta,2))-(sigma*pow(h,2)/pow(r,2))
c=(pow(sigma,2)*pow(h,4))/(3*pow(r,4))
d=pow((sigma*pow(h,2))/(3*pow(r,2)),3)
#Tolerance for f(x)
tol=1e-4

#defining 30 relaxation factors from 0.1 to 1.9
alpha=np.linspace(0.1,1.9,30)
#an empty list for storing no of iterations
iter=[]

#for loop for getting outputs for every value of relation factor
for rfac in alpha:
	#Initial guess for the minimum value of pressure
	pguess=10
	#defining a iteration variable
	itr=0

	#function calculating f(x)
	def f1(pguess,a,b,c,d):
		return (a*pow(pguess,3))+(b*pow(pguess,2))+(c*pguess)-d


	#function calculating f'(x)
	def f2(pguess,a,b,c):
		return (3*a*pow(pguess,2))+(2*b*pguess)+c

	#condition to run the loop until we reach the required value of f(x)
	while(abs(f1(pguess,a,b,c,d))>tol):

		#newton-raphson formula
		pguess=pguess-(rfac*(f1(pguess,a,b,c,d)/f2(pguess,a,b,c)))
		itr=itr+1
	
	#storing the iteration values in iter list
	iter.append(itr)

#minimum value of iterations amongst all values of relaxation factors
min_y=min(iter)

#plotting iteration vs relaxation factor graph
plt.plot(alpha,iter,marker='.',markerfacecolor='red',markersize='10')
plt.xlabel('Relaxation factor')
plt.ylabel('No of iterations')
plt.title('The relaxation factor corresponding to the minimum no of iterations i.e '+str(min_y)+' is '+str(alpha[iter.index(min_y)])+'')
plt.show()

	

For Part c:

The process is the same as in part one only with the additional use of a for loop to calculate the values of minimum pressure corresponding to each value of thickness of ice field. We have then used the tabulate module to present our output values in an organised tabular format. We have also plotted the pressure vs thickness graph to analyse how pressure varies with thickness.

Detailed code given below:

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

#an empty list for storing minimum pressures corresponding to each h
min_p=[]
#an empty list for storing no of iterations
iteration=[]
#thickness of ice field h
height=[0.6,1.2,1.8,2.4,3.0,3.6,4.2]

#for loop for getting outputs for every value of thickness
for h in height:
	#factor related to the width of ice wedge
	beta=0.5
	#tensile strength of the ice (1 psf = 144 psi)
	sigma=21600
	#size of the air cushion
	r=40

	#breaking down the calculations of equation 1 into smaller terms
	a=1-pow(beta,2)
	b=(0.4*h*pow(beta,2))-(sigma*pow(h,2)/pow(r,2))
	c=(pow(sigma,2)*pow(h,4))/(3*pow(r,4))
	d=pow((sigma*pow(h,2))/(3*pow(r,2)),3)
	#Tolerance for f(x)
	tol=1e-4
	#Initial guess for the minimum value of pressure
	pguess=10
	#defining a iteration variable
	itr=0
	#Relaxation factor
	alpha=1

	#function calculating f(x)
	def f1(pguess,a,b,c,d):
		return (a*pow(pguess,3))+(b*pow(pguess,2))+(c*pguess)-d


	#function calculating f'(x)
	def f2(pguess,a,b,c):
		return (3*a*pow(pguess,2))+(2*b*pguess)+c

	#condition to run the loop until we reach the required value of f(x)
	while(abs(f1(pguess,a,b,c,d))>tol):

		#newton-raphson formula
		pguess=pguess-(alpha*(f1(pguess,a,b,c,d)/f2(pguess,a,b,c)))
		itr=itr+1
	
	#storing the iteration values in iter list
	iteration.append(itr)
	#storing the min pressure values in min_p list
	min_p.append(pguess)

#organising and displaying the thickness, min pressure and iterations in the form of a table using tabulate module
data = zip(height,min_p,iteration)
print(tabulate(data,headers=["Thickness","Min Pressure","Iterations"],tablefmt='fancy_grid'))

#plotting pressure vs thickness graph
plt.plot(height,min_p,marker='.',markerfacecolor='red',markersize='10')
plt.xlabel('Thickness')
plt.ylabel('Pressure')
plt.title('Pressure vs Thickness graph')
plt.show()

Results and Discussion:

For Part a:

As per the Newton Raphson method, the minimum value of pressure required to break the ice field of 0.6ft is 4.239063886378783 which is achieved in 7 iterations.

For Part b:

As per the graph of iteration vs relaxation factor, it has been found that the minimum value of iteration is 7 corresponding to a relaxation factor of ~1.03

For Part c:

The table for the thickness of ice field, minimum cushion pressures and no. of iterations is given below:

Also, it can be seen from the Pressure vs Thickness graph that as thickness increases pressure also increase linearly with it.