Ice breaking

\"\"\"
Finding minimum pressure in Equation for Breaking Ice with Air cushion Vehicle using Newton-Raphson method
\"\"\"

import numpy as np
import math
import matplotlib.pyplot as plt

# Equation

def eqn(beta, sigma, r, h, p):
	e = np.arange(0.0,2.0)
	t1 = 1-np.power(beta,2)
	t2 = 0.4*h*np.power(beta,2) - sigma*np.power(h/r,2)
	t3 = np.power(sigma,2)*np.power(h/r,4)/3
	t4 = np.power(sigma*np.power(h/r,2)/3,3)
	e[0] = np.power(p,3)*t1+np.power(p,2)*t2+p*t3-t4    # function
	e[1] = 3*np.power(p,2)*t1+2*p*t2+t3                 # Derivative
	return e

# Iterative Newton Raphson method solving	

def sol(ini, iter, Tol, alpha, beta, sigma, r, h):
	p = ini
	for i in np.arange(1,iter+1):
		e = eqn(beta, sigma, r, h, p)
		if (abs(e[0])>Tol): p = p - alpha*(e[0]/e[1])
		else: break	
	return p,e,i

# Inputs
r 		= 40.0            								  # in feet
h 		= np.array([0.6, 1.2, 1.8, 2.4, 3.0, 3.6, 4.2])   # in feet	
ini 	= 23040.0   					# in Pounds per square foot
sigma 	= 150.0*144.0					# in Pounds per square foot
p 		= np.arange(0.0,len(h))			# in Pounds per square foot
iter 	= 100     	
Tol 	= 1e-5    	  
alpha 	= 1         	
beta 	= 0.5

print(\'Table\')
print(\'h  		p 				Value of function 	Number of iterations\')
for i in np.arange(0,len(h)):
	p[i],e,ni = sol(ini, iter, Tol, alpha, beta, sigma, r, h[i]) 
	# print(\'nNumber of iterations = {}, Equation value = {}\'.format(ni,e[0]))	
	print(\'{}		{:.6f}		{:.6f}			{}\'.format(h[i],p[i],e[0],ni))

# Plot
pressure = np.arange(-200.0,200.1,0.1)
e 		 = np.zeros((len(h), len(pressure), 2))

for j in np.arange(0,len(h)):
	for i in np.arange(0,len(pressure)):
		x = eqn(beta, sigma, r, h[j], pressure[i])
		e[j,i,0] = x[0]
		e[j,i,1] = x[1] 

plt.plot(h,p)
plt.xlabel(\'Thickness of ice(h) in feet\')
plt.ylabel(\'Cushion pressure(p) in psf\')
plt.grid()
plt.show()

plt.figure()
for i in np.arange(0,len(h)):
	plt.plot(pressure,e[i,:,0])
	plt.plot([-100,100],[p[i],p[i]])
	plt.title(\'h = {} , p = {}\'.format(h[i],p[i]))
	plt.xlabel(\'Pressure(p) in feet\')
	plt.ylabel(\'Value of function\')
	plt.grid()
	plt.show()

MODULE AND FUNCTION DEFINITION

The modules numpy, math, matplotlib.pyplot are imported

The function is defined in terms of beta, sigma, r , h, p where beta is related to the width of the ice wedge, sigma is the tensile strength of the ice in N/m^2, r is the size of the air cushion in m, h is the thickness of the ice field in m, p is cushion presssure in N/m^2. For returng evenly spaced values within a given interval using np.arrange between 0.0 and 2.0 and is read to a variable e. The whole equation is broke down into four terms t1, t2, t3, t4

t1 = 1-np.power(beta,2)

t2 = 0.4*h*np.power(beta,2) - sigma*np.power(h/r,2)

t3 = np.power(sigma,2)*np.power(h/r,4)/3

t4 = np.power(sigma*np.power(h/r,2)/3,3)

The function is read towards e[0] and its derivative is read towards to e[1] and e is returned.

NEWTON RAPHSON SOLVING

Primarily the counter of the iterations are created . The value of root is calculated using Newton Raphson method by defining an initial value p = ini. And for each operation of the loop the tolerance range is verified and loop continues till the condition breaks and the value of p, e, I are returned.

Defining the inputs r, h, ini, sigma, p, iter, Tol, alpha ( relaxation factor), beta.

COMPUTING THE OUTPUTS

The pressure p[i] for each value of h is calculated by the execution of a for loop which is governed up to a range by using np.arrange till len(h). The value of h, pressure p[i], Value of function e[0] and Iterations ni are printed as a table.

For returng evenly spaced values within a given interval of 0.1 using np.arrange between -200.0 and 200.1 and is read to a variable pressure. Two arrays with initial elements zero are created which is between length of matrix h and the length of pressure. Using for loop governed by the elements 0 and length of h using np.arrange. One more for loop is created which is governed by the elements 0 and length of pressure using np.arrange. The root x is calculated using the function and is read in e[j,i,0] and e[j,i,1] as x[0] and x[1].

PLOTTING THE RESULTS

Plot is created between h and p. And respective labels are given.

For plotting for different h , np.arrange is used. And plots pressure and e , p[i] are done. And respective labels are given.