Find Minimum Cushion Pressure with Newton-Raphson method of an Air Cushion Vehicle to break Ice in Python

Objective:

Write a Python program to find minimum cushion pressure of an Air Cushion Vehicle to Break Ice with Newton-Raphson method.

Governing Equation:

The factors involved are:

p - cushion pressure

h - thickness of the ice field

r - size of the air cushion

σ - tensile strength of the ice 

β - related to the width of the ice wedge

`alpha` - relaxation factor.

The solution will be obtained for the following cases:

Case-1: Find p for different values of h assuming a suitable relaxation factor.

Case-2: Find the optimal relaxation factor at h=0.6.

Air Cushion Vehicle (ACV): 

A hovercraft, also known as an air-cushion vehicle or ACV, is an amphibious craft capable of travelling over land, water, mud, ice, and other surfaces.

Hovercraft use blowers to produce a large volume of air below the hull, or air cushion, that is slightly above atmospheric pressure. The pressure difference between the higher pressure air below the hull and lower pressure ambient air above it produces lift, which causes the hull to float above the running surface. For stability reasons, the air is typically blown through slots or holes around the outside of a disk- or oval-shaped platform, giving most hovercraft a characteristic rounded-rectangle shape.

Newton-Raphson Method:

The Newton-Raphson method is one of the most widely used methods for root finding. It can be easily generalized to the problem of finding solutions of a system of non-linear equations, which is referred to as Newton's technique. Moreover, it can be shown that the technique is quadratically convergent as we approach the root.

Unlike the bisection and false position methods, the Newton-Raphson (N-R) technique requires only one inital value x₀, which we will refer to as the initial guess for the root. To see how the N-R method works, we can rewrite the function f(x) using a Taylor series expansion in (x-x₀):

where f'(x) denotes the first derivative of f(x) with respect to x, f''(x) is the second derivative, and so forth. Now, suppose the initial guess is pretty close to the real root. Then (x-x₀) is small, and only the first few terms in the series are important to get an accurate estimate of the true root, given x₀. By truncating the series at the second term (linear in x), we obtain the N-R iteration formula for getting a better estimate of the true root:

 Thus the N-R method finds the tangent to the function f(x) at x=x₀ and extrapolates it to intersect the x axis to get x₁. This point of intersection is taken as the new approximation to the root and the procedure is repeated until convergence is obtained whenever possible. Mathematically, given the value of x = x_i at the end of the ith iteration, we obtain x_i+1 as

We assume that the derivative does not vanish for any of the x_{k, k=0,1,..., i+1}. The result obtained from this method with x₀ = 0.1 for the equation of Example 1, x*sin(pi x)-exp(-x)=0, is graphically shown in Figure 2.

Case-1: Find p for different values of h assuming a suitable relaxation factor

Inputs:

Width of ice wedge (beta) = 0.5

Size of air cushion (r) = 40

Tensile strength In PSF = PSI*144 (sigma) = 21600

Thickness of ice (h) = [0.6,1.2,1.8,2.4,3,3.6,4.2]

Relaxation factor (alpha) = 1

Python Code:

"""
Find Minimum Pressure using Newton Raphson method by Giridhar
Case 1
"""

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

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


# f'(x)		
def fprime(p,beta,r,sigma,h):
	term7 = 3*pow(p,2)*(1-pow(beta,2))
	term8 = 0.8*p*h*pow(beta,2)
	term9 = 2*sigma*p*pow(h,2)/pow(r,2)
	term10 = (pow(sigma,2)*pow(h,4))/(3*pow(r,4))
	term = term7+term8-term9+term10				
	return term									

# Inputs

# Width of ice wedge
beta = 0.5 	
# Size of air cushion
r = 40	
# Tensile strength In PSF = PSI*144	
sigma = 21600 	
# Thickness of ice
h = [0.6,1.2,1.8,2.4,3,3.6,4.2] 
# relaxation factor
alpha = 1

p_guess = 20
tol = 1e-4
itr = 1
p = []
iter = []

for i in h:
	while(abs(f(p_guess,beta,r,sigma,i))>tol):
		p_guess = p_guess-alpha*(f(p_guess,beta,r,sigma,i)/fprime(p_guess,beta,r,sigma,i))
		itr = itr+1
	iter.append(itr)
	p.append(p_guess)

print(iter)
print(p)

for i in range(0,7):
	if i==0:
		print("height : Pressure")
	print(h[i],":",p[i])

#Table
headers = ['Thickness (h)','Iterations','Pressure (p in PSF)']			
table = tuple(zip(h,iter,p))
print(tabulate(table,headers=headers,tablefmt="fancy_grid"))

#Plotting 
plt.grid()
plt.plot(h,p)
plt.title("Min. Cushion pressure (vs) Ice Wedge Thickness")
plt.xlabel("Ice Wedge Thickness (ft)")
plt.ylabel("Min. Cushion Pressure (PSF)")
plt.savefig("Case_1")
plt.show()

Output:

Tabulated results of Pressure wrto h:

For 0.6ft thickness of ice 4.23 PSF minimum pressure is required.

Pressure vs Thickness Plot:

As the thickness increases the minimum cushion pressure required to break the ice wedge increases in a slight non-linear path. This can be assumed as a direct relationship between the minimum cushion pressure and thickness of the ice. As the thickness increases, more pressure is required to break the ice.

Python Code Explanation:

• Imported necessary modules required for code execution.
• Functions f and fprime(derivative) are defined according to governing equation and the value is returned.
• All required inputs are defined.
• h is defined in array.
• p and iter are defined with an empty array to find the values wrto h.
• A for loop is created to calculate pressure and number of iterations for each value of thickness of ice.
• Another for loop is created to print pressure for each value of thickness.
• Table properties are defined to print table with height, number of iterations and pressure values.
• The graph Minimum Cushion Pressure vs Thickness of Ice Wedge is defined using plt command.

Case-2: Find the Optimal Relaxation factor at h = 0.6

Python Code:

"""
Find Minimum Pressure using Newton Raphson method by Giridhar
Case 2
"""

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

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


# f'(x)		
def fprime(p,beta,r,sigma,h):
	term7 = 3*pow(p,2)*(1-pow(beta,2))
	term8 = 0.8*p*h*pow(beta,2)
	term9 = 2*sigma*p*pow(h,2)/pow(r,2)
	term10 = (pow(sigma,2)*pow(h,4))/(3*pow(r,4))
	term = term7+term8-term9+term10				
	return term									

# Inputs

# Width of ice wedge
beta = 0.5 	
# Size of air cushion
r = 40	
# Tensile strength In PSF = PSI*144	
sigma = 21600
# Thickness of ice
h = 0.6

# relaxation factor
alpha=np.linspace(0.1,1.9,20) 
tol = 1e-4
p = []
a = []
iter = []

for i in alpha:
	p_guess=20
	itr=1
	while((abs(f(p_guess,beta,r,sigma,h)))>tol):
		p_guess=p_guess-i*(f(p_guess,beta,r,sigma,h)/fprime(p_guess,beta,r,sigma,h))
		itr=itr+1
	iter.append(itr)
	p.append(p_guess)
	a.append(i)

# Optimum relaxation factor
print(iter)
print(p_guess)
print("The minimum iterations are:",min(iter))

index_min_iteration = iter.index(min(iter))
alpha_optimum=alpha[index_min_iteration]

print("Optimum Relaxation Factor=",alpha_optimum)

# Plotting
plt.grid()
plt.plot(a,iter)
plt.title('Relaxation Factor (vs) Number of iterations')
plt.xlabel("Relaxation Factor")
plt.ylabel("Number of iterations")
plt.savefig("Case_2")
plt.show()

Output:

Optimum Relaxation Factor obtained is 1.047 and the minimum number of iterations required for solution are 10.

Relaxation Factor vs Number of Iterations Plot:

Python Code Explanation:

• Imported necessary modules required for code execution.
• Functions f and fprime(derivative) are defined according to governing equation and the value is returned.
• All required inputs are defined.
• alpha is defined using linspace command and the value of alpha should be in between 1 and 1.9 with 20 values.
• p, a and iter are defined with an empty array.
• A for loop is created to calculate pressure, alpha and number of iterations
• Print command is used for minimum number of iterations and optimum relaxation factor output.
• The graph Relaxation Factor vs Number of Iterations is defined using plt command.

Conclusion:

All the required results are obtained using Newton-Raphson method which is mostly used for non-linear functions and the solution can be obtained in less number of iterations and the code is executed succesfully.

Reference:

1. https://en.wikipedia.org/

2. https://web.mit.edu/

Pic Credits: https://www.griffonhoverwork.com/