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

Objective:

In determining the minimum cushion pressure needed to break a given thickness of the ice
using an air cushion vehicle, Muller (“Ice Breaking with an Air Cushion Vehicle”) in
Mathematical Modeling: Classroom Notes in Applied Mathematics, SIAM 1987) derived
the equation

where p denotes the cushion pressure, h the thickness of the ice field, R the size of the air
cushion, σ the tensile strength of the ice, and β is related to the width of the ice wedge.
Take β = 0.5, r = 40 feet and σ = 150 pounds per square inch (psi). Construct a table
with h in the first column, and p in the second column for h = 0.6, 1.2, 1.8, 2.4, 3.0, 3.6
and 4.2 feet.

Newton Raphson:

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) = 0f(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)f(x), and you know the root you are looking for is near the point x = x_0x=x0​. Then Newton\'s method tells us that a better approximation for the root is

`x_(guss) = x_(guss)-alpha*[f(x)/f^\'(x)]`

x_guss is the intial guss value and finding the roots of the function. It is used for the non-linear equations.

Python Code:

I have written python code for the pressure required to break the Ice Brick and using matplotlib,numpy and tabulate modules. The function has derivated by the function of P. That has been split into terms.

Inputs

Case: 1

beta = 0.5
r = 40
sigma = 21600
p_guss = 30
height = [0.6,1.2,1.8,2.4,3.0,3.6,4.2]
tol = 1e-4
alpha = 1

In this case The thickness of the Ice Brick

Case: 2

beta = 0.5
r = 40
sigma = 21600
p_guss = 30

h = 0.6

tol = 1e-4
alpha = 1

In this case we have taken the relaxiation factor is 1

Case: 3

In this case, The relaxation Factor has been taken in the range of 0.1 to 2.

\'\'\'
	Breaking Ice with air cushion Vehicle
\'\'\'
import matplotlib.pyplot as plt
from tabulate import tabulate
import numpy as np


def f(p_guss,beta,r,sigma,h):
	term1 = pow(p_guss,3)*(1-pow(beta,2))
	term2 = (0.4*h*pow(beta,2)-(sigma*pow(h,2))/pow(r,2))*pow(p_guss,2)
	term3 = ((pow(sigma,2)*pow(h,4))/(3*pow(r,4)))*p_guss
	term4 = pow((sigma*pow(h,2))/(3*pow(r,2)),3)
	return (term1+term2+term3-term4)

def fprime(p_guss,beta,r,sigma,h):
	term1 = 3*pow(p_guss,2)*(1-pow(beta,2))
	term2 = 2*(0.4*h*pow(beta,2)-(sigma*pow(h,2))/pow(r,2))*p_guss
	term3 = ((pow(sigma,2)*pow(h,4))/(3*pow(r,4)))
	return (term1+term2+term3)
# Inputs

beta = 0.5
r = 40
sigma = 21600
p_guss = 30
height = [0.6,1.2,1.8,2.4,3.0,3.6,4.2]
tol = 1e-4
alpha = 1
itr = 1
p = []
it = []

# Case 1

for h in height:
	while (abs(f(p_guss,beta,r,sigma,h))>tol):
		p_guss = p_guss - alpha*(f(p_guss,beta,r,sigma,h)/fprime(p_guss,beta,r,sigma,h))
		itr= itr+1
	it.append(itr)
	p.append(p_guss)

headers = [\'thickness\',\'iterations\',\'pressure\']
mydata = tuple(zip(height,it,p))
print(tabulate(mydata,headers=headers,tablefmt=\'pipe\'))

# Case 2
\'\'\'
	Relaxiation Factor= 1 and thickness = 0.6
\'\'\'
p_guss = 30
h = height[0]
itr = 1

while (abs(f(p_guss,beta,r,sigma,h))>tol):
		p_guss = p_guss - alpha*(f(p_guss,beta,r,sigma,h)/fprime(p_guss,beta,r,sigma,h))
		itr = itr+1

print(p_guss)
print(itr)

# Case 3
\'\'\'
	The range of Realxiation factor
\'\'\'
alpha = np.arange(0.1,2,0.1)
h = height[0]
it_1 = []

for al in alpha:
	p_guss = 30
	itr = 1
	while (abs(f(p_guss,beta,r,sigma,h))>tol):
		p_guss = p_guss - al*(f(p_guss,beta,r,sigma,h)/fprime(p_guss,beta,r,sigma,h))
		itr = itr+1
	it_1.append(itr)
header = [\'alpha\',\'iterations\']
mydata = tuple(zip(alpha,it_1))
print(tabulate(mydata,headers=header,tablefmt=\'pipe\'))

plt.figure()
plt.plot(height,p)
plt.xlabel(\'Thickness\')
plt.ylabel(\'Pressure\')
plt.title(\'Thickness vs Pressure\')
plt.legend([\'Pressure\'])

plt.figure()
plt.plot(alpha,it_1)
plt.xlabel(\'Relaxiation Factor\')
plt.ylabel(\'iterations\')
plt.title(\'Alpha vs iterations\')
plt.legend([\'No of iteration for minimization\'])
plt.show()

Result:

Case: 1

In this case, I have written the code for finding the pressure for a different kind of Thickness. As seen above Data, gives how much pressure and Iteration took for Convergence the solution. The first one has less thickness so that the pressure also required very less compared to the last one and the number of iteration also required vey less it depends on the Relaxation Factor. If the relaxation factor changes above one there number of iteration also increase for convergence. As the below graph tells us thickness vs Pressure. Where the thickness increases pressure. Where the thickness increases as well as pressure increases.

Case: 2

In this case, I have written the code for finding the pressure for thickness 0.6. As below Data, gives how much pressure and Iteration took for Convergence the solution. for the thickness has taken 0.6. The optimal relaxation factor has taken 1 because it needs less iteration for converging the solution. When the relaxation factor increases more than 1 or less than 1 as well as number iteration increases.

Case:3

In this case, I have written the code for finding the pressure for a different kind of Relaxation Factor. As seen above Data, gives how much pressure and Iteration took for Convergence the solution. Since the relaxation factor 0.1 to 1 number of iterations decreases for converging the solution. The Relaxation factor 1.1 to 1.9 the number of iteration has been increased. So that for relaxation factor 1 gives an optimal solution less number of iterations. The below figure shows how the number iterations are required for converging depending on the Relaxation factor.