Breaking Ice with Air cushion Vehicle - Python script to find minimum Pressure with Newton-Raphson method and the optimum Relaxation Factor.

Problem Statement:

The problem is to write a python script solving the minimum value of pressure needed to break the ice of different thickness with air cushion vehicle using the Newton-Raphson method. The function given below defines the process and its influencing parameters.

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.

Assuming:

β = 0.5, r = 40 feet and σ = 150 pounds per square inch (psi).

Objectives:

1. To find out the value of pressure for h = 0.6.

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

3. To find the optimal relaxation factor for this problem with the help of a suitable plot.

Program:

Part 1:

# Newton Raphson Method
# Breaking Ice with Air Cushion Vehicle
# By Arjun Damodar

# 1. To mimimise the cushion pressure (p) required to break ice using Newton Raphson Method with relaxation factor 1

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

#Inputs
beta = 0.5                                    # Constant related to the width of the ice wedge
r = 40                                        # Size of the air cushion (in feet)
sigma = 150*144                               # Tensile strength of ice in psi; (1psi = 144 psf)
h = [0.6, 1.2, 1.8, 2.4, 3.0, 3.6, 4.2]       # Thickness of the ice (in feet)

def f(p,h):
	term1 = pow(p,3)*(1-pow(beta,2))
	term2 = pow(p,2)*((0.4*h*pow(beta,2))-(sigma*pow((h/r),2)))
	term3 = p*pow(sigma,2)*pow((h/r),4)/3
	term4 = pow((sigma*pow((h/r),2)/3),3)
	term5 = term1+term2+term3-term4
	return term5

def f_dash(p,h):
	term_1 = 3*pow(p,2)*(1-pow(beta,2))
	term_2 = 2*p*((0.4*h*pow(beta,2))-(sigma*pow((h/r),2)))
	term_3 = (pow(sigma,2)*pow((h/r),4))/3
	term_4 = term_1+term_2+term_3
	return term_4

# Newton Iteration
p_guess = 20  # Initial guess
alpha = 1     # Relaxation Factor
tol = 1e-5   # Tolerance
ite = 1
pressures = []           # Array to store minimum cushion pressure values
iterations = []          # Array to store no. of ierations 

for i in h:
	while(abs(f(p_guess,i))>tol):
		p_guess = p_guess - alpha*(f(p_guess,i)/f_dash(p_guess,i))
		ite = ite + 1
	pressures.append(p_guess)
	iterations.append(ite)

# Tabulation
data = tuple(zip(h,pressures,iterations))
headers = ['Thickness of Ice (feet)', 'Minimum Cushion Pressure (psf)', 'Number of Iterations']
print(tabulate(data,headers=headers,tablefmt='grid'))

# Plotting
plt.plot(h,pressures)
plt.xlabel('Thickness of ice (feet)')
plt.ylabel('Minimum Cushion Pressure Required to Break the Ice (psf)')
plt.show()

Part 2:

# Newton Raphson Method
# Breaking Ice with Air Cushion Vehicle
# By Arjun Damodar

# 2. To find the optimum relaxation factor

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

#Inputs
beta = 0.5                                    # Constant related to the width of the ice wedge
r = 40                                        # Size of the air cushion (in feet)
sigma = 150*144                               # Tensile strength of ice in psi; (1psi = 144 psf)
h = 0.6                                       # Thickness of the ice (in feet)

def f(p):
	term1 = pow(p,3)*(1-pow(beta,2))
	term2 = pow(p,2)*((0.4*h*pow(beta,2))-(sigma*pow((h/r),2)))
	term3 = p*pow(sigma,2)*pow((h/r),4)/3
	term4 = pow((sigma*pow((h/r),2)/3),3)
	term5 = term1+term2+term3-term4
	return term5

def f_dash(p):
	term_1 = 3*pow(p,2)*(1-pow(beta,2))
	term_2 = 2*p*((0.4*h*pow(beta,2))-(sigma*pow((h/r),2)))
	term_3 = (pow(sigma,2)*pow((h/r),4))/3
	term_4 = term_1+term_2+term_3
	return term_4

# Newton Iteration

alpha = np.linspace(0.01,1.99,500)     # Relaxation Factor
tol = 1e-5  # Tolerance
iterations = []

for i in alpha:	
	itr=0
	p_guess = 10
	h = 0.6  
	while(abs(f(p_guess))>tol):
		itr = itr + 1
		p_guess = p_guess - i*(f(p_guess)/f_dash(p_guess))
	iterations.append(itr)

# Finding optimum alpha

lowest_iteration = min(iterations)
position_of_the_lowest_iteration = iterations.index(lowest_iteration)
optimum_aplha = alpha[position_of_the_lowest_iteration]

# Results

print('Lowest Itertion is:',lowest_iteration)
print('Alpha value causing lowest iteration is',optimum_aplha)
plt.plot(alpha,iterations)
plt.xlabel('Relaxation Factor')
plt.ylabel('No. of Iterations')
plt.show()

Explanation Of the Code:

Part1:

1. Given values are assigned.
2. Tensile strength value given in psi is converted to psf by multiplying with 144.
3. The function 'f' is defined with the given formula.
4. Another function 'f_dash' is defined which returns the derivative of the function 'f'.
5. Factors required for the newton iteration are assumed and defined.
6. The initial guess for pressure is given as 20 and the relaxation factor is assumed to be 1.
7. Arrays to store pressure and iteration values are created.
8. The tolerance value is given as 1e-5 which is used as a condition in the 'while' loop defined inside the 'for' loop.
9. The 'for' loop runs for each value of h and 'while' loop performs the iteration in order to find the minimum pressure.
10. Corresponding pressure and iteration values are stored inside the arrays created.
11. Data required for the tabulation is stored and sorted using 'tuple(zip())' command.
12. Headers were defined according to the order and tabulation is created using the 'tabulate' command.
13. Plotting is done to get P vs h.

Part 2:

1. Instead of the 'h array' given in the first part, h is given a constant value of 0.6. Rest everything is same as of the first part till defining the function. The function is defined with a single variable.
2. A set of relaxation factors from 0.01 to 1.99 for 500 intervals were created using 'np.linspace' command and stored in the array 'alpha[]'.
3. An array for storing iteration value is created.
4. The 'for ' loop is assigned in order to run the loop for each alpha value and while loop which defined inside the for loop performs the iteration for each value of alpha and stores the no. of iteration taken for each.
5. To find the optimum relaxation factor or alpha, the lowest iteration value is found from the 'iterations array' using 'min()' command.
6. Then the position or index of that lowest value (element) is found from the list ('iterations array') using 'index()' command. The syntax is as follows: "listname.index(element)"
7. Optimum relaxation factor is found by using the index value in the array of relaxation factor(alpha[]).
8. Plotting is done for no.of iterations vs alpha.

Results:

Part 1:

Screenshot of the output:

Tabulation:

P V/s h:

Part 2:

Screenshot of the Output:

Optimum Alpha:

Optimum alpha is 0.998 which is causing the lowest iteration of 7.

No. of iterations V/s Alpha:

Observation and Conclusions:

1. The minimum pressured required to break the ice is increasing when the thickness of ice is increased. For h=0.6, the pressure required is equal to 4.23psf which is equal to 0.029psi.
2. The optimum relaxation factor is found to be 0.998 for the given set of alpha values ranging from 0.01 to 1.99. No. of iteration is observed to be high near 0.01 and 1.99 and observed to be less in the range of 0.25 to 1.75.
3. The iteration overshoots when the relaxation factor value is set to 2 and the loop runs infinitely without meeting the required tolerance.