Breaking Ice with Air cushion Vehicle - Find minimum pressure with Newton-Raphson method (Python)

AIM: To write a python script to solve for the minimum value of pressure using the Newton Raphson method.

OBJECTIVES :

1. To use the Newton Raphson method to find out the value of pressure for h = 0.6.

2. To find the optimal relaxation factor for this problem with the help of a suitable plot.

3. To tabulate the results of p for h = [0.6,1.2,1.8,2.4,3,3.6,4.2] assuming a suitable relaxation factor.

THEORY: Ice sheet can be broken into small slabs with a heavy air-cushion vehicle (ACV) traveling on the float ice sheet at the critical speed. Ice-breaking with ACV is a new ice-breaking method in the inland river which has an important application in the prevention of ice jam flood in the Yellow River. Dependent on ice thickness, an ACV can readily break the ice in flexure by making depression and send cracks through the ice. Once a crack is made, the path of least resistance for the ACV'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 setup a wave in the ice and send cracks through it. This would act similar to the dynamics of a transport truck on an ice road.

NEWTON RAPHSON METHOD : 

Newton's method, also called the Newton-Raphson method, is a root-finding algorithm that uses the first few terms of the Taylor series of a function $f\left(x\right)$ in the vicinity of a suspected root. Newton's method is sometimes also known as Newton's iteration, although in this work the latter term is reserved for the application of Newton's method for computing square roots.

The Taylor series of $f\left(x\right)$ about the point $x=x_0+\epsilon$ is given by-

$f(x_0+\epsilon )=f\left(x_0\right)+f^{\prime }\left(x_0\right)\epsilon +\frac{1}{2}{f}^{\text{'}\text{'}}\left(x_0\right){\epsilon }^2+....$ ,                                      -$\left(1\right)$

Keeping terms only to first order,

$f(x_0+\epsilon )\approx f\left(x_0\right)+f^{\prime }\left(x_0\right)\epsilon .$ ,                                                                               - $\left(2\right)$

Equation (2) is the equation of the tangent line to the curve at $(x_0,f\left(x_0\right))$, so $(x_1,0)$ is the place where that tangent line intersects the x-axis. A graph can, therefore, give a good intuitive idea of why Newton's method works at a well-chosen starting point and why it might diverge with a poorly-chosen starting point. This expression above can be used to estimate the amount of offset epsilon needed to land closer to the root starting from an initial guess $x_0$. Setting $f(x_0+\epsilon )=0$ and solving (2) for $\epsilon ={\epsilon }_0$ gives,

${\epsilon }_0=-\frac{f\left(x_0\right)}{f^{\prime }\left(x_0\right)}$ ,                                                                                                               - $\left(3\right)$

which is the first-order adjustment to the root's position. By letting $x_1=x_0+{\epsilon }_0$, calculating a new ${\epsilon }_1$ , and so on, the process can be repeated until it converges to a fixed point (which is precisely a root) using,

${\epsilon }_n=-\frac{f\left(x_n\right)}{f^{\prime }\left(x_n\right)}$ ,                                                                                                                -$\left(4\right)$

Unfortunately, this procedure can be unstable near a horizontal asymptote or a local extremum. However, with a good initial choice of the root's position, the algorithm can be applied iteratively to obtain,

 $x_{n+1}=x_n-\frac{f\left(x_n\right)}{f^{\prime }\left(x_n\right)}$                                                                                                    - $\left(5\right)$

for n=1, 2, 3, .... An initial point x_0 that provides safe convergence of Newton's method is called an approximate zero.

GOVERNING EQUATION USED :

$\Rightarrow f\left(x\right)=p^3(1-{\beta }^2)+(0.4h{\beta }^2-\frac{\sigma h^2}{r^2})p^2+\left(\frac{{\sigma }^2{h}^4}{3r^4}\right)p-{\left(\frac{\sigma {}^{}}{3{}^{}}\right)}^3$ ,

$\Rightarrow f\text{'}\left(x\right)=3p^2(1-{\beta }^2)+2p(0.4h{\beta }^2-\frac{\sigma h^2}{r^2})+\frac{{\sigma }^2{h}^4}{3r^4}$

where $p$ denotes the cushion pressure,

$h$ the thickness of the ice field,

$r$ the size of the air cushion,

$\sigma$ the tensile strength of the ice and

$\beta$ is related to the width of the ice wedge.

 SOLUTION STEPS :

1. At first, necessary modules such as math, matplotlib, NumPy are imported and required inputs are given. For the first part of code, an array of $h$ and an assumed value of $\alpha$ is considered according to the challenge statement.
2. Two different functions named $f\left(x\right)$ and $f\text{'}\left(x\right)$ are defined according to the formula given, and an empty array of a number of interactions and pressure(in psi and psft) are created.
3. For calculating minimum pressure for different values of ice thickness a for loop is executed taking an initial guessed value of pressure, tolerance, and a loop counter =1 which will increment after every iteration. Then a while loop is initiated inside the same for loop for checking whether all inputs in the formula if is greater than tolerance provided, a new value of pressure will be calculated according to the Newton Raphson method.
4. The calculated values of pressure and number of iterations are stored in their respective arrays using 'append' command and necessary statements are printed in the command window.
5. Graphs for pressure versus ice thickness are plotted using the calculated data set.
6. For calculating optimal relaxation factor($\alpha$) for a given thickness of ice same steps as above are repeated for except that now an array of $\alpha$ from 0.1 to 1.9 is taken with 50 points in between them and $h$ value is taken. Same for loop and while loop is again executed for checking whether all inputs in the formula if is greater than tolerance provided, a new value of pressure will be calculated according to the Newton Raphson method.
7. The number of iterations calculated is stored as an array and a graph is plotted between the same and relaxation factors.

PYTHON CODE : 

# A program to calculate pressure for given thickness of ice using Newton Raphson Method 
import math
import numpy as np
import matplotlib.pyplot as plt

# Inputs
# beta  is related to the width of the ice wedge
beta = 0.5
# r is the size of the air cushion (feet)
r = 40
# sigma the tensile strength of the ice (pounds per square inch (psi))
sigma = 150
# sigma the tensile strength of the ice (pounds per square foot)| 1foot = 12inches | 1 square feet = 12*12 square inches  
sigma = 21600
# h the thickness of the ice field
h = [0.6, 1.2, 1.8, 2.4, 3.0, 3.6, 4.2]
# Relaxation factor assumed
alpha = 1.1

def f(p,h,r,beta,sigma):
    
    term_1  = pow(p, 3)*(1-pow(beta, 2))
    term_2a = 0.4*h*pow(beta, 2)
    term_2b = (sigma*(pow(h,2))/pow(r,2))
    term_3  = pow(sigma,2)*(pow(h,4))/(pow(r,4))
    term_4  = (sigma/3)*(pow(h,2))/(pow(r,2))
      
    return term_1 + (term_2a - term_2b)*pow(p,2) + term_3*(p/3) - pow(term_4, 3)

def fprime(p,h,r,beta,sigma):
    
    term_5  = 3*pow(p,2)*(1-pow(beta, 2))
    term_6a = 0.4*h*pow(beta, 2)
    term_6b = (sigma*(pow(h,2))/pow(r,2))
    term_7  = pow(sigma,2)*(pow(h,4))/(pow(r,4))
    
    return term_5 + 2*(term_6a - term_6b)*p + term_7/3

# Null pressure array for storing calculated values after each iterations
pressure = []    # pound per square feet
pressure_0 =[]   # pound per square inches(psi)

# Null array for storing the number of iterations in future
num = []

for i in h:

    # Assuming initial pressure value
    p_guess = 100

    # Tolerance  
    tolerance = 1e-6

    # Loop counter for incrementing the number of iterations
    counter = 1

    while(abs(f(p_guess, i, r, beta, sigma)) > tolerance):
        p_guess = p_guess - alpha*(f(p_guess, i, r, beta, sigma)/fprime(p_guess, i, r, beta, sigma))
        counter = counter+1
        
    pressure.append(p_guess)
    pressure_0.append(p_guess/144)
    num.append(counter)
    
print('*Thickness of the Ice wedge :n',h)
print('*Cushion Pressure (pound per square feet) :n',pressure)
print('*Cushion Pressure (psi) :n',pressure_0)
print('*Number of iterations :n',num)

# Plotting graph
plt.figure(1)
plt.plot(h, pressure, marker = '.',color='black',markerfacecolor='red',markersize=10)
plt.title('Cushion pressure versus Ice field thickness')
plt.xlabel('Thickness of ice field (feet)')
plt.ylabel('Cushion Pressure (pound per square feet)')
plt.grid()
plt.show()

plt.figure(2)
plt.plot(h, pressure_0, marker = '.',color='black',markerfacecolor='green',markersize=10)
plt.title('Cushion pressure versus Ice field thickness')
plt.xlabel('Thickness of ice field (feet)')
plt.ylabel('Cushion Pressure (pound per square inches)')
plt.grid()
plt.show()



''' Part 2 of the program to calculate optimal relaxation factor for given thickness of ice field'''


# h_1 the thickness of the ice field
h_1 = 0.6
# Relaxation factor array assumed considering underrelaxing and overrelaxing values
alpha_array = np.linspace(0.1,1.9,50) 
# Null array for storing the number of iterations in future
num_1 = []

for j in alpha_array:

    # Assuming initial pressure value
    p_guess = 100

    # Tolerance
    tolerance = 1e-6

    # Loop counter for incrementing the number of iterations
    counter = 1

    while (abs(f(p_guess,h_1,r,beta,sigma)) > tolerance):
            
        p_guess = p_guess - (j*f(p_guess,h_1,r,beta,sigma))/fprime(p_guess,h_1,r,beta,sigma)
            
        counter = counter + 1

    num_1.append(counter)
    print(j)

# Plotting graph
plt.figure(3)
plt.plot(alpha_array,num_1,marker = '.',color='black',markerfacecolor='red',markersize=7.5)
plt.xlabel('Relaxation factor (alpha)')
plt.ylabel('No. of iterations')
plt.title('Iterations versus Relaxation factor')
plt.grid()
plt.show()


 

 OUTPUTS :

• The above program when executed yields pressure values for different ice thicknesses. From here it was observed that when ice thickness was 0.6, the pressure corresponding to it was 4.24 pounds per square feet or 0.0294 psi. 

• The variation of pressure with ice thickness is as follows-

                                                                        Figure(1)

                                                                        Figure(2)

• For finding optimal relaxation factor for given ice field thickness, the following plots were observed-

                                                                           Figure(3)

• From the figure(3), it can be concluded that curve rapidly decreases as the number of iterations decreases and relaxation factor increases and become constant around $\alpha$ =1 and again increases with the number of iterations and hence optimal relaxation factor lies in the vicinity of 1 as $\alpha$ could be somewhere in between of -

ERRORS: Python throws an error if there is indentation mismatching. It is important to take care of this. The error encountered was-

The reason for this error was that the pressure values which were iterated inside while loop was stored in the same loop without coming out of it. It was rectified by just removing spaces before line 61,62,63 so that they were out of while loop and now values calculated will be stored in that.

CONCLUSION: In this way, we can observe that the Newton Raphson method is a powerful tool for calculating solutions for non-linear functions as it is a more fast technique and gives approximately close solutions in fewer iterations compared to other methods.