Optimization and Finding Pressure values in Air Cushion Vehicle using Newton Rapshon Method

Objective:-

To determine the minimum cushion pressure to break the ice using an air cushion ice breaker with the help of the Newton Raphson Method.

In the given code I have found the results of the following points:-

1. The minimum pressure for different thickness.
2. Finding a suitable relaxation factor at h=0.6
3. Using the optimized relaxation factor determining the minimum pressure and comparing it with the results obtained for the relaxation factor of 1.

Theory:-

In determining the minimum cushion pressure needed to break a given thickness of 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 = cushion pressure,

h = thickness of the ice field,

R = size of the air cushion,

σ = tensile strength of the ice,

β = width of the ice wedge.

Basically, in this method, it is a repetitive nature to solve the problem. We have to follow certain no. of iterations to get the desired value.

Newton Raphson Method(NRM):-

Procedure:-

1. First I have imported the required import files from python library.
2. I have defined two function definitions that will return the values of two muller functions (f) and (df/dp). I have taken a scaling factor of 1e-26 as the functional values are of very high orders and to nullify it through scaling factor.
3. I have taken the variable called \'time\' to calculate the number of iterations occurred to obtain the minimum pressure for all values of thickness.
4. an array t of size 7*2 is taken to store the values of cushion pressure and thickness of the ice obtained using the Newton Raphson Method.
5. initially, the approximation is taken as 150 psi and relaxation factor of 1 is taken to find the minimum pressure with a tolerance of 1e-04.
6. the results are plotted and displayed in table format in the command window.
7. then in the second case to find the optimum relaxation factor, I have fixed the thickness to 0.6 feet and NRM is carried out on different relaxation factors with the same tolerance.
8. suitable plots have been plotted i.e Relaxation factor vs No.of iteration, Tolerance achieved vs Relaxation factor.
9. the upper steps are again done to find the minimum pressure against the thicknesses under the optimized relaxation factor and a plot is find out.
10. at least the difference between the time in the 2 cases are found out and displayed in the command window.
11. Most relaxation factor is the one that can produce the muller function with the least deviation from zero at the same time arriving at this solution in the least number of steps.

Given below is the required Code:-

\'\'\'
MINIMUM PRESSURE GOR AIR CUSHION ICE BREAKER
BY SHRABAN
\'\'\'

#Library Import
import numpy as np
import matplotlib.pyplot as plt 


#Functions
def f(p,h,r,beta,sigma):
 		a=pow(p,3)*(1-pow(beta,2))
 		b=((0.4*h*pow(beta,2))-(sigma*pow(h,2)/pow(r,2)))*pow(p,2)
 		c=pow(sigma,2)*pow(h,4)*p/(3*pow(r,4))
 		d=sigma*h*h/3*r*r

 		f_res=a+b+c-pow(d,3)
 		return ((f_res)*(1e-26))

#first derivative of the function
def df(p,h,r,beta,sigma):
		a=3*p*p*(1-pow(beta,2))
		b=2*p*((0.4*h*pow(beta,2))-(sigma*pow(h,2)/pow(r,2)))
		c=pow(sigma,2)*pow(h,4)/(pow(r,4)*3)

		f_res=a+b+c
		return ((f_res)*(1e-26))

#Input
beta=0.5
r=40		#in feets
r=r*12		#in inches
sigma=150	#pound/inch^2(psi)

#pressure minimisation for various ice thickness using Newton Rapson Method
time=0
t=np.zeros((7,2,1))		#initialising the table t

for i in range(1,8):
			#thickness array generation
			h=0.6*i 			#feet
			t[i-1,0]=h			#column 1 of table with thickness in feet
			h=h*12			#converting to inches
			#corresponding minimum pressure array generation
			p=150				#psi
			j=1
			tolerance=1e-04

			while (abs(f(p,h,r,beta,sigma))>tolerance):
				func= f(p,h,r,beta,sigma)
				dfunc=df(p,h,r,beta,sigma)
				p=p-(func/dfunc)		#NRM
				j=j+1
				time=time+1
			t[i-1,1]=p 				#column 2 of table with pressure in psi

#plotting
plt.plot(t[:,0],t[:,1])
plt.plot(t[:,0],t[:,1],\'*\')
plt.xlabel(\'THICKNESS(FEETS)\')
plt.ylabel(\'PRESSURE(PSI)\')
plt.title(\'MINIMUM PRESSURE VS THICKNESS\')
plt.show()

#printing table
print(\'|THICKNESS|		|PRESSURE|		\')
for i in range(0,7):
			print (\'|	\',t[i,0],\'	|\',t[i,1],\'	|\')


\'\'\'
Relaxation factor optimisation at h=0.6
\'\'\'

h=0.6*12
tolerance=1e-04
rf=0.2
RelaxationFactor=[]
counter=[]
func_val=[]
option=[]
Result=[100,rf]

while (rf<1.9):
		RelaxationFactor.append(rf)
		p=150
		j=0
		while (abs(f(p,h,r,beta,sigma))>tolerance):
			func=f(p,h,r,beta,sigma)
			dfunc=df(p,h,r,beta,sigma)
			p=p-rf*(func/dfunc)
			j=j+1

		counter.append(j)
		func_val.append(func)
		Op=i*abs(func)
		option.append(Op)
		if(Result[0]>Op):
			Result=[Op,rf]
		rf=rf+0.05

#plotting
plt.plot(RelaxationFactor,counter)
plt.plot(RelaxationFactor,counter,\'o\')
plt.xlabel(\'RelaxationFactor\')
plt.ylabel(\'Number of Iteration\')
plt.show()

plt.plot(RelaxationFactor,func_val,\'x\')
plt.plot([0,2],[0,0])
plt.xlabel(\'RelaxationFactor\')
plt.ylabel(\'Tolerance Achieved\')
plt.show()

print(\'Optimised Relaxation Factor\')
print(Result[1])

\'\'\'
PRESSURE MINIMISATION FOR VARIOUS ICE THICKNESS USING NEWTON RAPHSON METHOD(NRM) WITH OPTIMAL RELAXATION FACTOR
\'\'\'

ORF=Result[1]
time2=0
t2=np.zeros((7,2,1))		#initialising the table t

for i in range(1,8):
			#thickness array generation
			h=0.6*i 			#feet
			t2[i-1,0]=h			#column 1 of table with thickness in feet
			h=h*12			#converting to inches
			#corresponding minimum pressure array generation
			p=150				#psi
			j=1			
			tolerance=1e-04

			while (abs(f(p,h,r,beta,sigma))>tolerance):
				func= f(p,h,r,beta,sigma)
				dfunc=df(p,h,r,beta,sigma)
				p=p - ORF*(func/dfunc)	#NRM
				j=j+1
				time2 = time2 + 1
			t2[i-1,1]=p 				#column 2 of table with pressure in psi


#plotting
plt.plot(t2[:,0],t2[:,1],\'b-\')
plt.plot(t[:,0],t[:,1],\'r--\')
plt.xlabel(\'THICKNESS(FEETS)\')
plt.ylabel(\'PRESSURE(PSI)\')
plt.title(\'MINIMUM PRESSURE VS THICKNESS WITH OPTIMISED RELAXATION FACTOR\')
plt.legend([\'Optimised Relaxation\',\'Relaxation=1\'])
plt.show()

#printing
print(\'Time saved using RelaxationFactor in terms of no. of Iteration\')
print(time - time2)

Outputs:-

PLOT 1:-

PLOT 2:-

PLOT 3:-

PLOT 4:-

COMMAND WINDOW:-

The above plots and outputs give us a clear indication that the optimal relaxation factor(RF) can be obtained and results can be found in a much shorter value of time. Actually, RF is not dependent on the no. of iteration rather it depends on the function that has to be optimized. We see the results of pressure variation on the RF values at OPTIMAL and \'1\' value doesn\'t change much.

References:-

1. Google

2. Engineering Mathematics Books of authors B.S.GREWAL, KREYSIG, etc.

3. Youtube