Finding the minimum value of pressure for breaking ice with Air-Cushion vehicle using Newton-Raphson method

AIM

To write a program to solve for the minimum value of pressure for breaking ice with Air-Cushion vehicle using Newton-Raphson method

OBJECTIVES

• To use the Newton-Raphson method to find out the value of pressure for h = 0.6 ft.
• To find the optimal relaxation factor for this problem with the help of a suitable plot.
• 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

Air-Cushion Vehicle

Ice sheet can be broken into small slabs with a heavy air cushion vehicle (ACV) traveling on float ice sheet at the critical speed. Ice-breaking with ACV is a new ice-breaking method in inland rivers. Ice sheet cracks when the stress and strain in ice sheet accumulated to the critical value of failure. The critical speed of ACV should consist with the velocity of the surface wave around ACV.The relation between the parameters for breaking the ice is given by-

`p^3(1-beta^2)+(0.4hbeta^2-frac{sigmah^2}{r^2})p^2+frac{sigma^2h^4}{3r^4}p-(frac{sigmah^2}{3r^2})^3=0`

where

p= cushion pressure

h= thickness of ice field

r=size of the air cushion

σ=Tensile strength of the ice

β is related to the width of the ice wedge

The given values are-

β=0.5, r=40 feet,

σ=150 psi(pounds per square inch)=150*144 psf(pounds per square foot)

Newton-Raphson Method

The Newton-Raphson method (also known as Newton's method) is a way to quickly find a good approximation for the root of a real-valued function f(x)=0.

It uses the idea that a continuous and differentiable function can be approximated by a straight line tangent to it.

Suppose you need to find the root of a continuous, differentiable function f(x), and you know the root you are looking for is near the point `x = x_0`. Then Newton's method tells us that a better approximation for the root is

`x_1 = x_0 - \frac{f(x_0)}{f'(x_0)}`

This process may be repeated as many times as necessary to get the desired accuracy. In general, for any x-value `x_n`​, the next value is given by

`x_{n+1} = x_n -\frac{f(x_n)}{f'(x_n)}`

Note: the term "near" is used loosely because it does not need a precise definition in this context. However, `x_0​`should be closer to the root you need than to any other root (if the function has multiple roots).

Geometric Representation

Here is a picture to demonstrate what Newton's method actually does:

We draw a tangent line to the graph of `f(x)` at the point `x = x_n`. This line has slope `f'(x_n)`and goes through the point `(x_n, f(x_n))`. Therefore it has the equation `y = f'(x_n)(x - x_n) + f(x_n)`. Now, we find the root of this tangent line by setting `y = 0`and `x=x_{n+1}` for our new approximation. Solving this equation gives us our new approximation, which is

`x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}`

So, for the given function the required values of f and f' are-

`f=p^3(1-beta^2)+(0.4hbeta^2-frac{sigmah^2}{r^2})p^2+frac{sigma^2h^4}{3r^4}p-(frac{sigmah^2}{3r^2})^3`

`(df)/(dp)=f'(p)=3p^2(1-beta^2)+(0.4hbeta^2-frac{sigmah^2}{r^2})2p+frac{sigma^2h^4}{3r^4}`

CODE

Programming Language used- PYTHON

#A program to find minimum pressure with Newton-Raphson method
#importing required modules
import math
import matplotlib.pyplot as plt
import numpy as np

#p denotes the cushion pressure
#h denotes the thickness of the ice field
h=[0.6,1.2,1.8,2.4,3.0,3.6,4.2]
#r is the size of the air cushion(in foot)
r=40
#sigma is the tensile strength of the ice(in pounds per square foot)
sigma=150*144
#beta is related to the width of the ice wedge
beta=0.5

#function
def f(p,h):
	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

#Derivative of the function
def fprime(p,h):
	t1=3*pow(p,2)*(1-pow(beta,2))
	t2=((0.4*h*pow(beta,2))-(sigma*pow(h,2)/pow(r,2)))*2*p
	t3=((pow(sigma,2)*pow(h,4))/(3*pow(r,4)))
	return t1+t2+t3
	
#Performing the iterations
p_guess=10

#Creating a empty array
pressure=[]

#Tolerance 
tol=1e-4

#Relaxation faction
alpha=0.6

for i in h:
	while(abs(f(p_guess,i))>tol):
		p_guess=p_guess-alpha*(f(p_guess,i)/fprime(p_guess,i))
	pressure.append(p_guess)

print ('The minimun value of pressure,p for h=0.6 at relaxation factor 0.6 is :',pressure[0])
print('Tabulation of the results of p for h=[0.6,1.2,1.8,2.4,3,3.6,4.2] for a relaxation factor of 0.6')
print('Height : Pressure')
for i in range(len(h)):
	print(h[i],'   :',pressure[i])

#To find optimal relaxation factor
h1=0.6
alpha_all=np.linspace(0.1,1.5,50)

#Iterations
iterations=[]

for alpha in alpha_all:
	iter=1
	p_guess=10
	while(abs(f(p_guess,h1))>tol):
		p_guess=p_guess-alpha*(f(p_guess,h1)/fprime(p_guess,h1))
		iter=iter+1
	iterations.append(iter)

#Plotting
plt.plot(alpha_all,iterations)
plt.xlabel('Relaxation factor')
plt.ylabel('Number of Iterations')
plt.show()

PROCEDURE

1. First the required modules(i.e math, matplotlib and nump) are imported.
2. Then the array for different thickness is defined.
3. Then rest of the input parameters are declared i.e. size of air cushion,’r’, tensile strength of ice,’σ’ and β related to the width of the ice wedge.
4. Then the function is defined by dividing it into 4 terms.
5. Similarly, its derivative with respect to pressure,p is defined.
6. Then p_guess is assumed as 10 and relaxation factor,α is assumed as 0.6, taking a tolerance value of 1e-4(or 0.0001).
7. Then an empty array for pressure is created.
8. Using Newton-Raphson method the p_guess is calculated and appended into the empty pressure array till the absolute value of f(p_guess) is greater tolerance(using for loop and while loop respectively).
9. Then the minimum pressure for different values of ‘h’ is tabulated.
10. To find a optimal value of relaxation factor an array of 50 values between 0.1 to 1.5 was declared using ‘linspace()’ command.
11. Then an empty array was created for storing number of iterations.
12. Then taking h=0.6, and a counter variable for number of iterations set to1, the same process is carried out appending the number of iterations into the empty array.
13. Finally the graph is plotted between No. of iterations and different Relaxation factor.

RESULTS

Minimum value of Pressure for h=0.6ft and Tabulation of results for h=[0.6,1.2,1.8,2.4,3,3.6,4.2] at `alpha=0.6`

Plot for finding optimal Relaxation Factor

DISCUSSION

• The ‘numpy’ module is used for ‘linspace()’ command. Similarly the ‘matplotlib’ module is used for plotting graphs.
• The minimum value of cushion pressure for thickness of ice 0.6ft was found out. And minimum cushion pressure for different values of thickness is also tabulated.
• From the graph it can be seen that the number of iterations is minimum for value of relaxation factor around 1.
• So, we can say the optimal value of relaxation factor is nearly equal to 1.

CONCLUSION

It can be said that using Newton Raphson method we can easily minimise a function.Here, the trend of change between Relaxation factor and Iterations is observed as expected.