Breaking Ice using the Newton-Raphson Method on an Air Cushion Vehicle by calculating the Minimum-Pressure required to break the icy surface.

OBJECTIVE: BREAKING ICE USING THE NEWTON RAPHSON METHOD ON AN AIR CUSHION VEHICLE BY CALCULATING THE MINIMUM PRESSURE REQUIRED TO BREAK THE ICY SURFACE. 

THEORY: 

AIR CUSHION VEHICLE:

Recent experiments, using air cushion vehicles have demonstrated that they can be used successfully to break up ice on rivers. This is particularly useful at the ice-break-up time, when the ice may pack in one section of the river hindering the flow of water and eventually resulting in causing floods above the packed ice. The total weight of the platform is carried by the cushion pressure p(psi). At rest or slow speeds, the platform floats on the water and displaces its weight of water. If the pressure is sufficiently heavy, it causes radial cracks in the surface of the ice which emanates from the position where the applied pressure exceeds the tensile strength limit of ice at that point based on its thickness.

When the pressure is being applied to the surface of the ice, it results in causing stress to the icy surface which can give as:

`sigma = (3pr^2)/(h^2)[1-beta^(2/3)(1-(0.4h)/p)^(1/2)]`

The value of stress changes with respect to the applied pressure. When the stress level inside the ice exceeds the maximum tensile strength of the ice, it breaks. This pressure at which the ice starts breaking is the minimum pressure required to break the ice of that thickness.

To calculate the minimum pressure value for the given equation, we're converting the equation into pressure form as shown below: 

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

Where `beta = 0.5`, r = 40ft and `sigma = 150` pounds per square inch(psi) is given. 

This is a Non-linear equation and we have to find the roots of the equation to calculate the minimum pressure value required to break the ice. We know that root is synonymous with x-intercept and hence it is asked to find the point at which the function crosses the x-axis. Now, we need to set our equation to zero and solve for x. 

NEWTON RAPHSON METHOD: 

This is an iterative process that can approximate solutions to an equation with incredible accuracy and hence this method is used to approximate numerical solutions to equations that are too hard for someone to solve by hand. 

As we are needed to find the root of a continuous, differentiable function `f(x)`, and we know that the root value we're looking for is near the point `x = x_0`. Then, the Newton Raphson Method tells us a better approximation for the root as:

`x_1 = x_0-(f(x_0))/(f'(x_0))`

This process is repeated several times depending on the necessity for the degree of accuracy. For any x value $x_n$ the next value is given as: 

`x_(n-1) = x_n-f(x_n)/(f'(x_n))`

The trick of Newton's Method is to draw a tangent line to the graph, y = f(x) at the point(x1,y1). This tangent line is a good linear approximation to f(x) near x1, so our next guess x2 is the point where the tangent line intersects the x-axis. x2 is closer to the original solution. 

Here, the equation for f(x) is given as:

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

Where, 

`beta` = Width of ice wedge,

h = Thickness of the ice,

r = Size of the air cushion

p = Cushion pressure

`sigma` = Tensile strength of ice

PYTHON CODE:

FINDING THE RELATION BETWEEN PRESSURE THAT IS NEEDED TO BREAK THE ICE ACROSS VARYING THICKNESSES OF ICE: 

# Breaking Ice using Newton Raphson Method for an Air Cushion Vehicle
# By Sagar_Biswas
# Minimum Pressure required across Varying Thickness to Break Ice
# Required Libraries are Imported
import matplotlib.pyplot as plt
import numpy as np

# Inputs Parameters
h = [0.6, 1.2, 1.8, 2.4, 3.0, 3.6, 4.2]  # Array containing multiple values for thicknesses of ice
s = 150 * 144  # Tensile Strength of ice in psi (converted to Pound Per Square Foot by Multiplying 144)
b = 0.5  # 'Beta' stands for Width of Ice-Wedge
r = 40  # Size of Air-Cushion section in ft
tol = 1e-5  # Considering the Permissible Tolerance as the Convergence Criteria for the Pressure-Equation
p_guess = 150  # Providing an initial guess value for pressure
i = 1  # Iteration Counter Variable


# Pressure-Equation for the Air Cushion Vehicle
def f(p, h):
    return ((pow(p, 3) * (1 - pow(b, 2))) + ((0.4 * h * pow(b, 2) - (s * pow(h, 2) / pow(r, 2))) * pow(p, 2)) +
            (pow(s, 2) * pow(h, 4) * p / (3 * pow(r, 4))) - (pow((s * pow(h, 2) / (3 * pow(r, 2))), 3)))


# 1st Order Derivative for the Pressure-Equation
def f_prime(p, h):
    return ((3 * pow(p, 2) * (1 - pow(b, 2))) + 2 * p * (0.4 * h * pow(b, 2) - (s * pow(h, 2) / pow(r, 2))) + (
            pow(s, 2) * pow(h, 4) / (3 * pow(r, 4))))


min_pressure = []  # Creating a list to store the values of minimum pressure
iter = []  # Creating a list to store the number of iterations being created before convergence


for x in h:  # Using for-loop to calculate the results for each value of thickness of the ice-sheet
    alpha = 1  # Relaxation Factor used to Under-relax or Over-relax the operation to control the number of iterations
    i = 1  # Loop-Counter Variable Initial Value
    # Using while loop to keep the operation going till the tolerance value is met or the solution has converged
    while abs(f(p_guess, x)) > tol:  # Calculating the absolute value
        p_guess = p_guess - alpha * (f(p_guess, x) / f_prime(p_guess, x)) # Calculating guess values from the equations
        i = i + 1  # counter variable is set to increase by +1 each time the loop runs
    iter.append(i)  # Storing all the values of i in the list iter
    min_pressure.append(p_guess)  # minimum pressure values to break ice

print('Various Thickness of Ice:',h)
print('Minimum Pressure required to break the Ice respective to the thickness:',min_pressure)

plt.plot(h,min_pressure,'--',color='red', linewidth=3)
plt.xlabel('Varying Thickness of Ice (feets)')
plt.ylabel("Minimum Air Cushion's Pressure (psi)")
plt.title('Air Cushion Pressure VS Varying Thicknesses of Ice')
plt.savefig("minpressure.-thickness_ice.png")
plt.show()

Code is explained Step-by-Step inside the Python Code itself using Comments. 

PLOTTING RESULTS: 

It is evident from the plot that as the thickness of the ice increases, the pressure required to break its surface increases.

FINDING OUT THE OPTIMUM RELAXATION FACTOR WHEN THE THICKNESS OF ICE IS 0.6:

# Breaking Ice using Newton Raphson Method for an Air Cushion Vehicle
# By Sagar_Biswas
# Finding Relation between No. of Iterations and Successive Relaxation(Alpha)
# Finding Relation between Pressure being applied & Successive Relaxation
# Required Libraries are Imported
import matplotlib.pyplot as plt
import numpy as np
import math

# Input Parameters
h = 0.6  # We are going to find out the optimum
r = 40  # Size of Air-Cushion section in ft
s = 150 * 144  # Tensile Strength of ice in psi (converted to Pound Per Square Foot by Multiplying 144)
b = 0.5  # 'Beta' stands for Width of Ice-Wedge


# Pressure-Equation for the Air Cushion Vehicle
def f(p, h):
    return ((pow(p, 3) * (1 - pow(b, 2))) + ((0.4 * h * pow(b, 2) - (s * pow(h, 2) / pow(r, 2))) * pow(p, 2)) +
            (pow(s, 2) * pow(h, 4) * p / (3 * pow(r, 4))) - (pow((s * pow(h, 2) / (3 * pow(r, 2))), 3)))


# 1st Order Derivative for the Pressure-Equation
def f_prime(p, h):
    return ((3 * pow(p, 2) * (1 - pow(b, 2))) + 2 * p * (0.4 * h * pow(b, 2) - (s * pow(h, 2) / pow(r, 2))) + (
            pow(s, 2) * pow(h, 4) / (3 * pow(r, 4))))


alp = np.linspace(0.1, 1.9, 15)  # Creating multiple values of Alpha using linspace

iter = []  # Creating a list to store the number of iterations
P_val = []  # Creating a list to store the varying values of Pressure with respect to Varying values of Alpha

for x in alp:
    p_guess = 50
    tol = 1e-5
    i = 1
    while abs(f(p_guess, h)) > tol:
        p_guess = p_guess - (x * f(p_guess, h)) / f_prime(p_guess, h)
        i = i + 1

    P_val.append(abs(f(p_guess, h)))
    iter.append(i)

Lowest_Number_of_Iterations=iter.index(min(iter))
Optimum_Alpha=alp[Lowest_Number_of_Iterations]
print("Optimized Value of Alpha:->",Optimum_Alpha)
print("Minimum_Required_Iterations:->",min(iter))

plt.figure(1)
plt.plot(alp, iter, color='red', linewidth=2)
plt.xlabel('Successive Relaxation(Alpha)--->')
plt.ylabel('No. of iterations--->')
plt.title('Successive Relaxation(Alpha) VS No. of Iterations')
plt.show()

plt.figure(2)
plt.plot(alp, P_val, color='green', linewidth=2)
plt.xlabel('Successive Relaxation(Alpha)--->')
plt.ylabel('Varying Pressure Values--->')
plt.title('Successive Relaxation(Alpha) VS Varying Pressure Values ')
plt.show()

 Code is explained Step-by-Step inside the Python Code itself using Comments.

PLOTTING RESULTS:

From the above plot, we can understand that when the relaxation factor(alpha) is near the value of '1' then we can obtain the results using the least number of iterations. If the value of the relaxation factor(alpha) increases or decreases from this value then the number of iterations shoots up and hence increases the time to reach convergence.

In the above plot, we have plotted the relationship between the Successive Relaxation(Alpha) and Optimum Pressure values and we can observe that when we're using the value of 0.9 as the relaxation factor, only then we're getting the value of the pressure equation nearing the value of '0'.

We have also calculated the exact Optimal-Value of Alpha for our Program which leads to requiring the least number of iterations before convergence.

RESULTS & OBSERVATIONS: 

1. No. of Iterations required to reach the point where Minimum Pressure is the least decreases with increment in the successive relaxation factor. Even that relaxation factor after reaching a certain point known as 'Optimum Relaxation-Factor, the no. of iterations required starts increasing again.

2. For increment in the thicknesses of the ice, the Minimum-Pressure required to break them increases.