Breaking Ice with Air Cushion Vehicle and Finding minimum Air Cushion pressure using Newton Raphson Method

Breaking Ice with Air Cushion Vehicle and

Finding minimum Air Cushion pressure using Newton Raphson Method

• Concept of Breaking Ice with Air cushion vehicle:

• • In the region of low temperature river covers with thick layer of Ice causing the difficulty in waterways and movement of ships. To remove ice a special vehicle called Air cushion vehicle loaded with the Ice breaking system using Air cushioning are used. These vehicles the pressure of the air cushion is sufficient to depress the water and force air below the ice field.

• Objective:
  
  • Calculate minimum cushion pressure required to break the ice by using one of the Iterative methods namely Newton Raphson method using the equation as below:

    `p^3 (1-β^2 )+(0.4hβ^2-(σh^2)/r^2 ) p^2+((σ^2 h^4)/(3r^3 ))p-((σh^2)/r^2 )^3 = 0`

    Where, p is the minimum air cushion pressure in psf
    β is the term related to ice wedge in
    h is the thickness of the ice field in feet
    r is the size of air cushion in feet
    σ is the tensile strength of the ice in psf

  • To understand how minimum pressure changes with the values of thickness of Ice wedge. Also, the effect of air cushion size in both the parameters.
  • Plot the graph depicting minimum cushion as a function of thickness of Ice wedge and size of Air cushion.
  • Plot the graph of Number of iterations versus relaxation factor.
  • Determining the optimum value of relaxation factor used in the iterative method and understanding the behaviour of process when relaxation factor is 1.

• Code: (To find Minimum cushion pressure)

• • Defining the inputs and importing the various modules like math for mathematical calculation and matplotlib for plotting graph.
  • Nested list for p is created.
  • r and h as defined as list since they take wide range of values.
  • r takes values from 20 to 80 with step of 10 while h takes values from 1.6 to 4.2 with the step of 0.6.
  • σ and ß takes single value
  • list of color is created which will then be used to specify the different plots as a legends.
  • iteration values are stored in the nested list created name iter_values.

# Program to solve a problem of ice cushion breaking using Newton raphson problem
# By: Namit

import matplotlib.pyplot as plt
import math
import matplotlib.patches as mpatches

# Inputs

# thickness of the ice feild (feet)
H = [0.6, 1.2, 1.8, 2.4, 3.0, 3.6, 4.2]

# size of air cushion (feet)
R =[20, 30, 40, 50, 60, 70, 80]

r20 = []
r30 = []
r40 = []
r50 = []
r60 = []
r70 = []
r80 = []

c_size = [\'r = 20\', \'r = 30\', \'r = 40\', \'r = 50\', \'r = 60\', \'r = 70\', \'r = 80\']

# Nested List for pressure range at different height for various sizes of cushion
p_range = [r20, r30, r40, r50, r60, r70, r80]


iter_values = [[],[],[],[],[],[],[]]

color = [\'green\', \'blue\', \'red\', \'purple\', \'magenta\', \'orange\', \'pink\']

i = 0

• • Defining the function which calculates the value of equation as a function of p i.e. f(p)
  • Separating equation in separate terms to make calculation easy
    t1 = `p^3 (1-β^2 )`
    t2 = `(0.4hβ^2-(σh^2)/r^2 ) p^2`
    t3 = `((σ^2 h^4)/(3r^3 ))p`
    t4 = `( (σh^2)/r^2 )^3`
  • Defining the function which calculates the value of derivative of equation as a function of p i.e. fprime(p)
  • Finding derivative of each term with respect to p
    t1_prime = `3p^2 (1-β^2 )`
    • t2_prime = `(0.4hβ^2-(σh^2)/r^2 )2p`
      t3_prime = `((σ^2 h^4)/(3r^3 ))`
      t4_prime = 0
  • A nested for loop runs for various values of R and H, also while loop runs to do the iterative calculations using Newton Raphson for finding the minimum cushion pressure p.
  • First for loop access the values of r , nested for loop access the values of h and lastly the while loop helps in iterative process for p using the two functions f(p,h,r) and fprime(p,h,r).
    
    `p=p-a( f(p,h,r)/(fprime(p,h,r)) )`
    
    Here, a is the relaxation factor alpha. also we take a =1
    
  • Initial guess value of p is taken as 200 psf

for r in R:
    for h in H:
        beta = 0.5                           #Related to width of Ice wedge
        sigma = 150*144                      #Tensile Strength of Ice wedge (psf)
        p = 200                              #Cushion pressure (psf)
        a = 1                                #Relaxation factor
        def f(p, h, r):

            # formula is seperated in different terms denoted as t1, t2, t3, t4, t5

            t1 = pow(p, 3) * (1 - pow(beta, 2))
            t2 = ((0.4 * h * pow(beta, 2)) - ((sigma * pow(h, 2)) / pow(r, 2))) * pow(p, 2)
            t3 = ((pow(sigma, 2) * pow(h, 4)) / (3 * pow(r, 4))) * p
            t4 = pow((sigma * pow(h, 2)) / (3 * pow(r, 2)), 3)

            return t1 + t2 + t3 - t4

        def fprime(p, h, r):

            # first derivative of terms with repect to p is

            t1_prime = 3 * pow(p, 2) * (1 - pow(beta, 2))
            t2_prime = 2 * p * ((0.4 * h * pow(beta, 2)) - ((sigma * pow(h, 2)) / pow(r, 2)))
            t3_prime = ((pow(sigma, 2) * pow(h, 4)) / (3 * pow(r, 4)))
            t4_prime = 0

            return t1_prime + t2_prime + t3_prime - t4_prime

        tol = 1e-4
        iter = 1
        #Iterative Process (Newton Raphson Method)
        while abs(f(p, h, r)) > tol:
            p = p - a * (f(p, h, r)/fprime(p, h, r))
            iter += 1
        p_range[i].append(abs(p))
        iter_values[i].append(iter)
    i = i+1

• • creating Table of h and p for a = 1 and r = 40 using for loop
  • We have 7 values for h and p both so the loop runs for 7 times executing the print statement where values of h and p are printed side by side and concatenated with space for better alingnment

# Creating table for h and p at r = 40 feet
print(\'\\n\"Values of Minimum cushion pressure corresponding to the different values of thickness of Ice feild at r = 40 feet are:\"  \')

print(\'\\nThickness of Ice Field (h) : Minimum Cushion Pressure(p)\')

for j in range(0,7):
    print(\'         \' + str(H[j]) + \'              \' + \' : \' + \'       \' + str(r40[j]))

• • Graph of p vs h is plotted for different values of r using the for loop using the counter variable y.

# Plotting the graph
for y in range(0, 7):
    print(\'\\nvalue of cushion pressure at different heights for \' + str(c_size[y]) + \' is \' + \'\\n\',str(p_range[y]))
    print(iter_values[y])
    plt.plot(H, p_range[y], color = color[y])


#Customizing the plots
green_patch = mpatches.Patch(color=\'green\', label = \'r = 20\')
blue_patch = mpatches.Patch(color=\'blue\', label = \'r = 30\')
red_patch = mpatches.Patch(color=\'red\', label = \'r = 40\')
purple_patch = mpatches.Patch(color=\'purple\', label = \'r = 50\')
magenta_patch = mpatches.Patch(color=\'magenta\', label = \'r = 60\')
orange_patch = mpatches.Patch(color=\'orange\', label = \'r = 70\')
pink_patch = mpatches.Patch(color=\'pink\', label = \'r = 80\')
plt.legend(handles=[green_patch, blue_patch, red_patch, purple_patch, magenta_patch, orange_patch, pink_patch])


# Adding Label
plt.xlabel(\'Thickness of Ice Feild(h) (feet)\')
plt.ylabel(\'Minimum Cushion Pressure(p) (psi)\')

#Labelling and title of the graph
plt.suptitle(\" Minimum Cushion Pressure(p) vs Thickness of Ice Feild(h)\")
plt.title(\'p as a function of h and r\')
plt.show()

• Code: (To find optimal value of relaxation factor)

• • For this program the values of a varies from 0.1 to 1.8.
  • while h is taken as 0.6 and initial p guess is 20.
  • The code consists two functions f(p) and fprime(p) which helps in calculating the minimum pressure value using iterative method of newton raphson and the number of iterations in which required to reach the solution is stored in the iter_values. This process is repeated for every value of alpha defined using a for loop.

\'\'\'
Program to find the optimal value of relaxation factor alpha (a) for iterative method of air cushion problem
and plot alpha vs number of iterations
\'\'\'
# By: Namit


import matplotlib.pyplot as plt
import math
import numpy as np
import matplotlib.patches as mpatches



alpha = np.linspace(0.1, 1.8, 100)                           # Relaxation factor


# Inputs
h = 0.6  # Thickness of the Ice feild (feet)
beta = 0.5  # Related to width of Ice wedge
sigma = 150*144  # Tensile Strength of Ice wedge (psf)

r = 40  # Size of the air cushion (feet)


def f(p):
    # formula is seperated in different terms denoted as t1, t2, t3, t4, t5

    t1 = pow(p, 3) * (1 - pow(beta, 2))
    t2 = ((0.4 * h * pow(beta, 2)) - ((sigma * pow(h, 2)) / pow(r, 2))) * pow(p, 2)
    t3 = ((pow(sigma, 2) * pow(h, 4)) / (3 * pow(r, 4))) * p
    t4 = pow((sigma * pow(h, 2)) / (3 * pow(r, 2)), 3)

    return t1 + t2 + t3 - t4


def fprime(p):
    # first derivative of terms with repect to p is

    t1_prime = 3*pow(p, 2)*(1 - pow(beta, 2))
    t2_prime = 2*p*((0.4*h*pow(beta, 2)) - ((sigma * pow(h, 2)) / pow(r, 2)))
    t3_prime = ((pow(sigma, 2)*pow(h, 4))/(3*pow(r, 4)))
    t4_prime = 0

    return t1_prime + t2_prime + t3_prime - t4_prime

p_range = []
iter_values = []
tol = 1e-6

for a in alpha:
    iter = 1
    p = 200   # Cushion pressure (psf)
    # Iterative Process (Newton Raphson Method)
    while abs(f(p)) > tol:
        p = p - a * (f(p) / fprime(p))
        iter = iter + 1
    p_range.append(p)
    iter_values.append(iter)

• • After the values of Iterations are stored for every alpha value we need the optimal value of alpha which gives the minimum cushion pressure in minimum number of iterations.
  • To achieve this in-built function namely min() function is used which help us to find the minimum value of the iteration from the list of values stored in iter_values.
  • Index value of this minimum iteration is stored in variable which is then used to access the corresponding value of alpha to find the optimal value.

#Finding minima of graph (optimal value of alpha)

iter_min = min(iter_values)
alpha_pos = iter_values.index(iter_min)
opt_alpha = alpha[alpha_pos]
iter_min = min(iter_values)
alpha_pos = iter_values.index(iter_min)
opt_alpha = alpha[alpha_pos]
print(\'Optimal value of Relaxation Factor is \' + str(round(opt_alpha)))
print(\'Minimum number of iterations required to find minimum cushion pressure is  \' + str(round(iter_min)))

• • Graph of Number of iterations vs alpha is plotted which depicts the variations in number of iterations according to the value of alpha
  • Local minima of the graph is also plotted showing the minimum value of iterations and corresponding value of alpha is its optimal value.

#plotting graph of No. of iteration vs alpha

plt.title(\'Number of Iterations vs Relaxation Factor(alpha)\')
plt.plot(alpha, iter_values)
plt.xlabel(\'Relaxation factor (alpha)\')
plt.ylabel(\'Number of iterations\')

plt.show()

#plotting the optimum value of alpha (Local Minima of graph)

plt.title(\'Local Minima of the f(alpha, No. of iterations) \')
plt.plot(alpha, iter_values)
plt.xlabel(\'Relaxation factor (alpha)\')
plt.ylabel(\'Number of iterations\')
plt.plot(opt_alpha, iter_min, marker = \'*\' ,color= \'red\', markersize = 6)
plt.xlim(0.8, 1.2)
plt.ylim(0, 75)

plt.show()

• Output:
  • From plot of minimum cushion pressure versus thickness of ice wedge we can conclude that the pressure varies as a function of thickness of ice wedge and size of air cushion forming a parabolic shape.

• • Below is the table of showing the corresponding p values for various h values.

• • So, it is clear that for fixed value of r (cushion size) with increase in ice thickness the minimum cushion pressure required also increases parabolically. Also, for the same of value of h pressure required for breaking ice decreases.
  • The plot of p vs h for r = 20 is much steeper or has a higher slope which means with small change in thickness of ice wedge the minimum cushion pressure increases rapidly/noticeably while for r = 70 and 80 this plot almost has zero slope which indicates that p is almost constant for all values of h.

• • The plot of Number of Iterations vs alpha has a U-shaped curve which depicts that the iterations required to achieve minimum cushion pressure reduces with the increase in value of alpha but upto certain limit i.e. upto alpha = 1, after which number of iterations starts increasing.

• • Local minima of the plot is found which helps us determine the optimal value of alpha for minimum value of iterations.This is shown as a red mark in the plot of Local minima.

• • Minimum number of iterations comes out to be 17 and optimal value of alpha is 1.