Breaking ice with air cushion vehicle - finding minimum cushion pressure using newton-raphson method

Aim: To calculate minimum pressure required to break ice with air cushion vehicle using Python

Theory:

The relation between cushion pressure and thickness of ice field is given by the equation,

where,

p = cushion pressure,

h = thickness of ice field,

β is related to the width of the ice wedge = 0.5,

r = size of the air cushion = 40 feet and

σ = tensile strength of the ice = 150 psi

From the above equation value of p can be found for different values of h using Newton-Raphson method.

Newton-Raphson Method:

The Newton-Raphson method is a way to quickly find a good approximation for the root of a real-valued function f(x) = 0. It uses the idea that a continuous and differentiable function can be approximated by a straight line tangent to it. Tgis is

Suppose you need to find the root of a continuous, differentiable function f(x), and you know the root you are looking for is near the point x=x0​. Then Newton's method tells us that a better approximation for the root is,

`x_1 = x_0 -alpha \frac{f(x_0)}{f'(x_0)}`

where, `alpha` = Relaxation factor

The new value gives a better approximation for function and this value can be used again as an input to newton's formulae to find next better approximation. After repeating this for few iterations, solution converges. The relaxation factor determines how change in value of x takes places. If it is less than 1 (under relaxation), the value of the x change by lesser value in each iteration and if it is more than 1, the value of x changes by higher value in each iteration. The value of relaxation factor helps in reaching convergence in least iterations.

Python Program:

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

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

#Define function to calculate equation terms
def equ_terms(beta,r,h,sigma):
    term = []
    term.append(1 - pow(beta,2))
    term.append(0.4*h*pow(beta,2) - sigma*pow(h,2)/pow(r,2))
    term.append(pow(sigma,2)*pow(h,4)/(3*pow(r,4)))
    term.append(sigma*pow(h,2)/(3*pow(r,2)))
    return term

#function
def f(beta,r,h,sigma,p):
    term = equ_terms(beta,r,h,sigma)
    return term[0]*pow(p,3) + term[1]*pow(p,2) + term[2]*p - pow(term[3],3)

#first derivative of function
def fprime(beta,r,h,sigma,p):
    term = equ_terms(beta,r,h,sigma)
    return 3*term[0]*pow(p,2) + 2*term[1]*p + term[2]

# input parameters
beta = 0.5
r = 40
sigma_psi = 150
h = 0.6

# converting psi to pound per sqaure foot(psf)
sigma_psf = 150*144

# Newton Raphson parameters
p_guess = 25
tol = 1e-6

# Starting iterations....
while(abs(f(beta,r,h,sigma_psf,p_guess))>tol):
    p_guess = p_guess - f(beta,r,h,sigma_psf,p_guess)/fprime(beta,r,h,sigma_psf,p_guess)

print('Minimum value of p for h = 6ft is ' + str(p_guess) + ' psf')

# Optimal relaxation factor
alpha_ = np.linspace(0.5,1.5,200)
iterations = []
for alpha in alpha_:
    p_guess = 25
    i = 1
    while(abs(f(beta,r,h,sigma_psf,p_guess))>tol):
        p_guess = p_guess - alpha*(f(beta,r,h,sigma_psf,p_guess)/fprime(beta,r,h,sigma_psf,p_guess))
        i = i + 1
    iterations.append(i)    

plt.figure()
plt.plot(alpha_,iterations)
plt.xlabel('Relaxation Factor')
plt.ylabel('iterations')
plt.title('Optimization of Relaxation factor')
plt.savefig('Optimization of Relaxation factor')

min_iter = min(iterations)
opt_alpha = alpha_[iterations.index(min_iter)]
print('Optimat value of Relaxation factor :'+str(opt_alpha))
print('iterations at optimal Relaxation factor:'+str(min_iter))

# Using optimal relaxation factor to find p
h_ = [0.6,1.2,1.8,2.4,3,3.6,4.2]
p_calculated = []
for h in h_:
    p_guess = 25
    while(abs(f(beta,r,h,sigma_psf,p_guess))>tol):
        p_guess = p_guess - opt_alpha*(f(beta,r,h,sigma_psf,p_guess)/fprime(beta,r,h,sigma_psf,p_guess))
    p_calculated.append(p_guess)

plt.figure()
plt.plot(h_,p_calculated)
plt.xlabel('Thickness of the ice field(ft)')
plt.ylabel('Cushion pressure(psf)')
plt.title('Results')
plt.savefig('Results')

#vertically stack two 1D arrays into one 2D array
results = np.vstack((h_,p_calculated)).T

# convert numpy array to pandas dataframe
result = pd.DataFrame(data=results, columns=['Thickness of ice field(ft)','Cushion pressure(psf)'])
# start index with 1
result.index += 1
print(result)

A function equ_terms is defined to calculate coefficients of the equation. It return coefficients as an array.

A function f is defined to calculate value of function. A function fprime is defined to calculate first derivative of the function.

To have uniformity of dimension sigma_psi is converted into sigma_psf by multiplying it with 144.

Initial guess is set to 25 and tolerance value is 1e-6 and value of p is calculated using newton-raphson method for h=0.6 using while loop. The while loop continues till the value to function is greater than tolerance value.

To find optimal value of alpha, an array of value is created between 0.5 and 1.5 using np.linspace for 200 values. The iterations for different value of alpha are calculated using for loop and a plot is made for alpha vs iterations. The value of alpha for min value of itertions if optimal value of relaxation factor. Min value of iteration is found using min command and value of alpha is found using index of min iteration using index command.

The value of p is found for given values of h using optimal relaxation factor and plot is made for h vs p. The two 1D arrays,h_ and p_calculated, are stacked together in a single array using 2D array. The numpy array is converted using pandas DataFrame using pd.DataFrame to tabulet the data with suitable column title. The result.index+=1 is used to increases index value of each row by 1, this makes idex start with 1.The result is printed using print command.

Output:

Conclusion:

The optimal value of relaxation factor is 0.99749. The result chart shows the value of cushion pressure for different values of ice field thickness.