Finding Minimum Air Cushion Pressure With Newton-Raphson Method in Python
1) Air Cushion Vehicle
[source: https://bit.ly/2z5LgIl]
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 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. Typically this cushion is contained within a flexible "skirt", which allows the vehicle to travel over small obstructions without damage. [1]
[source: https://en.wikipedia.org/wiki/File:Hovercraft_-_scheme.svg]
2) Ice Breaking with an ACV and Muller's Equation
Ice breaking can be done conventionally by ice-breaking ships or by an ACV. A conventional icebreaker uses mechanical advantage and the vessels mass to break ice downward in flexural stress. Dependent on ice thickness, a hovercraft can readily break the ice in flexure by making a depression and send cracks through the ice. Once a crack is made, the path of least resistance for the hovercraft's air curtain gets under the ice and reduces the supporting water under the ice. The ice then fails under its own weight. At speed, an ACV can set up a wave in the ice and send cracks through it. [2]
E.R Muller, in his paper titled Ice Breaking with an Air Cushion Vehicle, explored the ice breaking mechanism under a slow moving air cushion vehicle and derived a formula that determines the minimum air cushion pressure that an ACV needs to develop in order to break an ice of particular thickness. [3]
In order to find the minimum cushion pressure, the following equation must be satisfied:
`p^3(1-beta^2) + (0.4beta^2 - sigmah^2/r^2)p^2 +sigma^2h^4/(3r^4)p -(sigmah^2/(3r^2))^3 =0`
Where,
`p` is the cushion pressure in psi
`beta` is a constant
`sigma` is the tensile strength of the ice
`r` is a term related to the dimension of the ACV
`h` is the thickness of the ice (in some cases the cushion height)
Note:
To keep the dimensionality consistent, we can either convert the sigma into psf (per square feet; 1 psi = 12^2 psf = 144 psf) or we can simply multiply 'r' and 'h' by 12 to convert them into inches. However, note that in the Muller's equation, 'r' and 'h' always appear in a ratio (h/r). Therefore, it doesn't matter in what unit you keep 'r' and 'h' as long as they are in the same units. However, changing sigma to psf will always give the outputs in psf and not in psi.
3) Newton-Raphson Method
In numerical analysis, the Newton-Raphson method is a method for finding successively better approximations to the roots (or zeroes) of a real-valued function. In other words, it is an algorithm to find the roots of an equation.
The Newton–Raphson method in one variable is implemented as follows:
The method starts with a function f defined over the real numbers x, the function's derivative f ′, and an initial guess x0 for a root of the function f. If the function satisfies the assumptions made in the derivation of the formula and the initial guess is close, then a better approximation x1 is:
`x_1 = x_0 - f(x_0)/(f'(x_0)`
Geometrically, (x1, 0) is the intersection of the x-axis and the tangent of the graph of f(x) at (x0, f (x0)). The process is repeated as:
`x_(n+1) = x_(n) - f(x_(n))/(f'(x_(n))`
until a sufficiently accurate value is reached. [4]
4) Using Python To Find The Minimum Cushion Pressure
In Python, we can use the Newton-Raphson method to find the minimum cushion pressure to break the ice. For this, we use the Muller's equation and write it in the form:
`y(p) = a*p^3 + b* p^2 +c.p -d`
Where:
`a = 3(1-beta^2)`
`b = 0.4 beta^2 - sigmah^2/r^2`
`c = sigma^2h^4/(3r^4)`
`d = (sigmah^2/(3r^2))^3`
Keeping a, b, c and d constant, the equation primarily has one variable, 'p', pressure which has real values. Thus, the equation is real-valued and the Newton-Raphson method can be applied to find the values of the cushion pressure where the equation becomes (approximates) zero.
import numpy as N
import matplotlib.pyplot as M
beta = 0.5
r = 40*12 #converting to inches
sigma = 150 #psi
h = 4.2*12 #converting to inches
''' To keep the dimensionality consistent, we can either convert the sigma into psf (per sqaure feet; 1 psi = 12^2 psf = 144psf)
or we can simply multiply 'r' and 'h' by 12 to convert them into inches. However, note that in the Muller's equation, r and h always
appear in ratio (h/r). Therefore, it doesn't matter in what unit you keep r and h as long as they are in the sam,e units.
However, changing sigma to psf will always give the outputs in psf and not in psi
'''
A = 1 - pow(beta,2)
B = (0.4*pow(beta,2)) - (sigma*pow(h/r,2))
C = (pow(sigma,2)*pow(h/r,4))/3
D = pow((sigma*(pow(h/r,2))/3),3)
p = N.linspace(0,1.3,100)
y = A*pow(p,3) + B*pow(p,2) + C*p - D
M.plot(p,y)
M.title('Visualising The Root of The Function')
M.ylabel('f(p) @h = 4.2')
M.xlabel('Pressure (psi)')
M.grid()
M.show()If a plot between `f(p)` at h = 4.2 and `p = (0,1.4)` is drawn, we see that the graph crosses the X-axis somewhere between 1.2 and 1.4 psi.
Fig. 4.1: Visualising The Root of Muller's Equation
Newton-Raphson root finding method is an iteration technique. As such, we can either use a 'for' loop: where we specify the number of iteration the code shall run. However, using for-loop may not always give better results. This is because of the reason that the loop runs for a specified number of iterations which can lead to the following three problems:
- The approximation might not converge when the number of iterations is sufficiently less.
- The roots (approximation) might have converged way before the last iteration. This leads to unnecessary time consumption.
- Lastly, there is a chance that the approximation might oscillate around the root and the last iteration overshoots the best approximation. This is one of the main failures of the method.
To avoid such a situation, we can set a tolerance level (a minimum value), which when achieved stops the program. For this, a 'while' loop is ideal. For example, in a particular case, a minimum of 0.0001 may be an acceptable value of minima. Thus setting a tolerance of 0.0001 a while loop will iterate only till this required condition is obtained. (see the main program)
Program:
1) The derivative of the function can either be written in the exact form or, it can be achieved by the ab-initio derivative method i.e.,
#ab-initio differentiation
def f_prime(x):
dh = 0.0001
f_prime = (f(x+dh) - f(x))/dh
return f_prime
''' where f(x) is the Muller's function at all values.'''
2) The Newton-Raphson equation can have an additional factor, called the under-relaxation factor.
`x_(n+1) = x_(n) - alpha( f(x_(n))/(f'(x_(n))))`
Suppose, in a function, that there are two roots close enough. Using the general Newton's method to find the root might not provide both the roots as in the process, it might skip the points where the other root lies. Thus, in order to lower this 'amount of skipping points' process, the under-relaxation factor is introduced. This factor basically makes sure that in the iterative process, the method progresses (converges) slowly so that no roots are skipped. If the multiplicity of the function is known, then it is equal to alpha, otherwise, it can be calculated after few iterations. [5]
Main Program:
import matplotlib.pyplot as M #a module to plot like Matlab
import numpy as np
#1a)Defining the constant values
beta = 0.5
r = 40*12 # ACV dimension in inches
sigma = 150 #tensile strength (psi)
alpha = 0.8 #relaxation factor
#1b) Height Array
Heights = [0.6,1.2,1.8,2.4,3.0,3.6,4.2] #list of heights in feet
H = np.array(Heights)*12
#making a numpy 1D array and multiplying elementwise by 12. Elements in inches
#1c)Output Heading
print ('%s %30s %20s' % ('Thickness,h(ft)', 'Minimum Pressure,p(psi)', 'Iterations'))
'''%Ns creates num number of spaces in the output, where N is an integer'''
#2) Defining Various Functions:
#2a) Muller's f(x)
def f(x):
function = A*pow(x,3) + B*pow(x,2) + C*x - D
return function
#2b) derivative, f_prime
def f_prime(x):
f_prime = 3*A*pow(x,2) + 2*B*x + C
#else use ab-initio method
return f_prime
#2c) newton-raphon formula
def newton(x):
return (x - alpha*(f(x)/f_prime(x)))
Min_Pressure = [] #preallocation
#3) Iterations
i = 0
for h in H:
#3a) Condesing constant terms
A = 1 - pow(beta,2)
B = (0.4*pow(beta,2)) - (sigma*pow(h/r,2))
C = (pow(sigma,2)*pow(h/r,4))/3
D = pow((sigma*(pow(h/r,2))/3),3)
guess= 200 #initial guess
tol = 1e-4 #tolerance
counter = 1 #counts the number of iterations
#3b) program loop
while (abs(f(guess)) > tol): #abs(f) because in our case f(p) is always > 0
guess = newton(guess) #evaluates the new value for x
counter = counter + 1
Min_Pressure.append(guess) #appends the value of min pressures at different heights
#3c) Printing the values
print('%10s'% Heights[i],'%20s' %'','%.5f' %guess, '%20s' %counter)
i = i+1
'''
Prints the values of height, guess and iterations at a partiuclar height, and places them below the respective headings.
Note the different method of syntax as compared to Output Heading
'''
#4) Plotting
Alp = 'alpha = ' + str(alpha)
print('\n',Alp)
M.plot(Heights,Min_Pressure, 'r')
M.title('Min Cushion Pressure Vs Ice Height')
M.xlabel('Height (feet)')
M.ylabel("Min Cushion Pressure (psi)")
M.grid()
M.show()
Output in Sublime Window:
Output in Command Prompt:
The plot of Min Cushion Pressure Vs Ice Thickness/Cushion Height:
Fig. 4.2 Minimum Cushion Pressure
From the above plot, we can clearly see that as the Ice thickness or cushion height increases, the ice breaking pressure also increases.
Changing alpha from under-relaxation(0.1, 0.5), 1 to over-relaxation(1.5, 2):
`alpha` 0.1
`alpha` 0.5
`alpha` 1
`alpha` 1.5
`alpha` 2
From the above outputs, it is clear that the number of iterations increases when alpha deviates from 1.5. Hence, the ideal value of alpha for the said function lies around 1.5.
5) Optimum Value of Relaxation Factor
To calculate the optimum value of `alpha`, we can fix a particular ice thickness, and vary the values of `alpha` from 0 to 2 (Strictly speaking, ` 0 < alpha <=2`). This is because, we have already seen from the hit-and-trials for the values of alpha, that at `alpha` equal to 2, the iterations for heights 0.6 and 1.2 feet are very small, but for other heights, the number of iterations is awfully large.
Selecting the average of the heights, i.e., 2.4 feet, the program to calculate the minimum number of iterations, which will give the optimum value of `alpha` is given below:
Code:
import numpy as np
import matplotlib.pyplot as plt
beta = 0.5
r = 40*12
sigma = 150
h = 2.4*12 #average value of heights
Alpha = np.linspace(0.01,2,200) #Alpha should have non zero values. A zero value will yeild infiite iterations
# Muller's function, f(x)
def f(x):
function = A*pow(x,3) + B*pow(x,2) + C*x - D
return function
# Derivative, f_prime
def f_prime(x):
#f_prime = 3*A*pow(x,2) + 2*B*x + C # exact form
dh = 0.000001
f_prime = (f(x+dh)-f(x))/dh #ab-initio method
return f_prime
# Newton-Raphon formula
def newton(x,alpha):
return (x - alpha*(f(x)/f_prime(x)))
# Iterations
iterations=[] #preallocating, defining data type
for i in range(len(Alpha)):
#
# Condesing constant terms in Muller's equation
A = 1 - pow(beta,2)
B = (0.4*pow(beta,2)) - (sigma*pow(h/r,2))
C = (pow(sigma,2)*pow(h/r,4))/3
D = pow((sigma*(pow(h/r,2))/3),3)
guess= 200 #initial guess
ig = guess # storing initial guess value as 'guess' changes in the next steps
tol = 1e-4 #tolerance
counter = 1 #counts the number of iterations
# Program loop
while (abs(f(guess)) > tol): #abs(f) because in our case f(p) is always > 0
guess = newton(guess,Alpha[i]) #evaluates the new value for x
counter = counter + 1
iterations.append(counter)
# Plotting, no. of iterations vs alpha values)
plt.plot(Alpha,iterations)
plt.title('Varaition of No. of Iterations With Alpha')
plt.xlabel('Alpha')
plt.ylabel('No. of Iterations')
plt.grid()
plt.show()
#Finding the minimum value in the 'iterations' list.
minimum = min(iterations) #min() fucntion calculates the minimum value in a list
alpha_optimum=[]
for j in range(len(Alpha)):
if iterations[j] == minimum: #execute only if the value is minimum
alpha_optimum.append(Alpha[j]) #append the values at the same index
alpha_optimum = np.array(alpha_optimum) #converting into a numpy array for proper decimal numbers
'''
#confirming the values
alpha_optimum1 = (alpha_optimum-0.01)*100
#as Alpha starts ahead with 0.01 and not 0.00. alpha_optimum1 stores the indices where min(iterations) occur.
#printing the itertion values
for i in range(len(alpha_optimum1)):
print(iterations[int(alpha_optimum1[i])])
'''
print('The minimum number of iterations to converge to the solution, when height = %2.1f, initial guess = %s and tolerance = %s is %s which occurs when alpha is = ' %(h,ig, tol,minimum),alpha_optimum)
'''
>>>
The minimum number of iterations to converge to the solution, when height = 28.8, initial guess = 200 and tolerance = 0.0001 is 13 which occurs when alpha is = [1.46 1.64 1.67 1.71 1.76]
'''Output:
Fig. 5.1: Variation of No. Iterations With Alpha
It can be seen from the above figure and the output results that the number of iterations to converge to the solution is extremely large (infinite) at 0 and decreases abruptly till 0.25 and then decreases gradually up to `alpha` equal to 1.46, where it attains its first minima. The curve hovers around this value until the last minima at `alpha` equal to 1.76. The number of iterations again increases largely towards `alpha` equal to 2 and above.
Thus, the optimum value for alpha `alpha` when height is 2.4 feet, initial guess for Newton-Raphson method is 200 and tolerance is 1.46, 1.64, 1.67, 1.71 or 1.76.
***End***