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

Objective:

To solve the given equation and find the minimum pressure required to break a given thickness of the ice by using \'Newton Raphson\' iterative technique in python.

The formula used to determine the minimum pressure required to break a given thickness of the ice is given by

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.
β = 0.5, r = 40 feet and σ = 150 pounds per square inch (psi) and

h = [0.6, 1.2, 1.8, 2.4, 3.0,3.6,4.2]

To solve the above equation newton raphson method is used. Newton Raphson method is an iterative method which is used to find the roots. Multivariate newton raphson method used to solve stiff systems of coupled equations to find the roots. Here we are using the explicit differential because it is not a stiff system. The formula for the raphson method is given by

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

The first partial derivative concerning p for the above equation is given by

`f^\' (x)=3p^2 (1-β^2 )+2p*(0.4hβ^2-(σh^2)/r^2 )+ (σ^2 h^4)/(3r^4 )-((σh^2)/(3r^2 ))^3`.

Code:

• Import the matplotlib module for plotting.
• Define the f1 function (equation for calculation of minimum pressure as shown above) and fprime function(a derivative of the main function(f1)).
• In the range of h solve the equation using the newton raphson method.
• Solve the newton raphson method by using the while condition and since it is an iterative process there is always error associated with it which is called tolerance and the absolute value of the function should always greater than the tolerance value.
• Create a table for height and corresponding pressure. Define the for loop of the range of 7 and print the values of the height and pressure as shown.

import matplotlib.pyplot as plt 

def f1(p,beta,sigma,h,r):
	return pow(p,3)*(1-pow(beta,2))+(0.4*h*pow(beta,2)-(sigma*pow(h,2)/pow(r,2)))*pow(p,2)+ p*((pow(sigma,2)*pow(h,4))/(3*pow(r,4))) - pow(((sigma*pow(h,2))/(3*pow(r,2))),3)
def fprime(p,beta,sigma,h,r):
	return 3*pow(p,2)*(1-pow(beta,2))+ (0.4*h*pow(beta,2)-sigma*(pow(h,2)/pow(r,2)))*2*p +((pow(sigma,2)*pow(h,4))/(3*pow(r,4))) - pow(((sigma*pow(h,2))/(3*pow(r,2))),3)

beta = 0.5 # width of ice wedge
r = 40 # size of air cuhsion
sigma = 150 # tensile strength of the ice
h = [0.6,1.2,1.8,2.4,3.0,3.6,4.2]
p_guess=100
alpha = 1
tol=1e-8
#print(tol)
iter =0
p=[]
for H in h:
	#print(\'Thickness of the ice:\',H)
		
	while (abs(f1(p_guess,beta,sigma,H,r))>tol):
		p_guess=p_guess-alpha*(f1(p_guess,beta,sigma,H,r)/(fprime(p_guess,beta,sigma,H,r)))
		iter=iter+1
		 
	p.append(p_guess)
\"\"\"
Creating the table
\"\"\"
for i in range(0,7):
	if i==0:
		print(\"s.no\\tHeight\\t\\t\\tPressure\")
	print(i+1,\'\\t\\t\',h[i],\'\\t\',p[i])
#print(\'thickness of ice:\',h,\'Value of the pressure :\',p)


plt.plot(h,p,\'*\')
plt.xlabel(\'Thickness of the ice\')
plt.ylabel(\'Pressure\')
plt.show()

Results:

The pressure required to break the thickness of the height is given below.

Plot:

From the above results, It is clear that as the thickness of the ice (h) increases the pressure required to break the ice also increased proportionally. Hence, the lowest pressure required to break the ice

Optimum7 relaxation factor:

Code:

• Import the matplotlib module for plotting and numpy module for the array operations.
• Define the f1 function (equation for calculation of minimum pressure as shown above) and fprime function(a derivative of the main function(f1)).
• In the range of alpha solve the equation using the newton raphson method.
• Solve the newton raphson method by using the while condition and since it is an iterative process there is always error associated with it which is called tolerance and the absolute value of the function should always greater than the tolerance value.
• Define the empty array called iterations and append it through the c.
• Create a plot between the alpha and iterations.

import matplotlib.pyplot as plt 
import numpy as np 
def f1(p,beta,sigma,h,r):
	return pow(p,3)*(1-pow(beta,2))+(0.4*h*pow(beta,2)-(sigma*pow(h,2)/pow(r,2)))*pow(p,2)+ p*((pow(sigma,2)*pow(h,4))/(3*pow(r,4))) - pow(((sigma*pow(h,2))/(3*pow(r,2))),3)
def fprime(p,beta,sigma,h,r):
	return 3*pow(p,2)*(1-pow(beta,2))+ (0.4*h*pow(beta,2)-sigma*(pow(h,2)/pow(r,2)))*2*p +((pow(sigma,2)*pow(h,4))/(3*pow(r,4))) - pow(((sigma*pow(h,2))/(3*pow(r,2))),3)

beta = 0.5 # width of ice wedge
r = 40 # size of air cuhsion
sigma = 150 # tensile strength of the ice
alpha = np.linspace(0.1,2,50)
#print(tol)
p=[]
iterations=[]
for a in alpha:
	#print(a)
	
	p_guess=100

	c=0
	h = 0.6
	tol=1e-6

	while (abs(f1(p_guess,beta,sigma,h,r))>tol):
		p_guess=p_guess-a*(f1(p_guess,beta,sigma,h,r)/(fprime(p_guess,beta,sigma,h,r)))
		c=c+1
	iterations.append(c)
			 
	p.append(p_guess)

plt.figure(1)
plt.plot(alpha,iterations)
plt.show()
 

Result:

In general, the optimization of the relaxation factor is chosen where the number of iterations is the lowes. Here, the range of alpha is from 0.1 to 2 and the optimum relaxation factor is around 1.7. As shown in the plot, after 1.7 the newton raphson method takes more iterations to converge the solution.