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

Python for Mechanical Engineers

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

Objective:

To determine 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”) using the following 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).

Given that:

1. To 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.6and 4.2 feet. (relaxation factor = 1)

2. To find the optimum relaxation factor by plotting a graph. (h = 0.6)

3. Explaining the results with the help of appropriate relaxation factor.

Air Cushion Vehicle:

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”)

The equations:

Iteration: Numerical methods are often, of a repetitive nature. These consist in repeated execution of .The same process where at each step the results of the preceding step is used. This is known as

iteration process and is repeated till the result is obtained to a desired degree of accuracy.

Relaxation factor: - it represents the value of reduction of the current load step for the restart providing the solution fails to converge. A new value of the assumed load step is found from the express.

Newton Raphson method:

The Newton-Raphson method uses an iterative process to approach one root of a function. The specific root that the process locates depends on the initial, arbitrarily chosen x-value.

- Depending on the conditions under which we are attempting to solve this equation, several of the variables may be changing.

So, it may be necessary to use partial derivatives.

- For the purposes of this example, we are assuming that height is the only thing changing, and that these values are all functions of time. This avoids the use of a partial derivative, we simply differentiate

all variables with respect to pressure, as shown above.

Procedure:

1. Find f'(x) and define NR equation.
2. Guess an initial value of x for the first iteration.
3. substitute x in NR equation and calculate x*.
4. If |x-x*|>tolerance, stop iteration and outputthe root:x*
5. If number of iterations reaches an assumed maximum, stop.
6. Else let x=x* and repeat the steps from 3 until a condition in 4 or 5 is satisfied.

Coding:

# Program to calculate pressure for given thickness of ice using Newton Raphson Method
import math
import numpy as np
import matplotlib.pyplot as plt
# Taking inputs
# Beta is related to the width of the ice wedge
beta = 0.5
# r is the size of the air cushion (feet)
r = 40
# sigma the tensile strength of the ice (pounds per square inch (psi))
sigma = 150
# sigma the tensile strength of the ice (pounds per square foot), 1foot = 12inches, 1 square feet =
#12*12 square inches
sigma = 21600
# h the thickness of the ice field
h = [0.6, 1.2, 1.8, 2.4, 3.0, 3.6, 4.2]
# Relaxation factor assumed
alpha = 1.1
# Defining the governing function
def f(p,h,r,beta,sigma):
	term_1 = pow(p, 3)*(1-pow(beta, 2))
	term_2a = 0.4*h*pow(beta, 2)
	term_2b = (sigma*(pow(h,2))/pow(r,2))
	term_3 = pow(sigma,2)*(pow(h,4))/(pow(r,4))
	term_4 = (sigma/3)*(pow(h,2))/(pow(r,2))
	return term_1 + (term_2a - term_2b)*pow(p,2) + term_3*(p/3) - pow(term_4, 3)
# Defining the derivative of the governing function
def fprime(p,h,r,beta,sigma):
	term_5 = 3*pow(p,2)*(1-pow(beta, 2))
	term_6a = 0.4*h*pow(beta, 2)
	term_6b = (sigma*(pow(h,2))/pow(r,2))
	term_7 = pow(sigma,2)*(pow(h,4))/(pow(r,4))
	return term_5 + 2*(term_6a - term_6b)*p + term_7/3
# Null pressure array for storing calculated values after each iterations
pressure = [] # pound per square feet
pressure_0 = [] # pound per square inches(psi)
# Null array for storing the number of iterations in future
num = []
# Creating a for loop for interation and performing the newtion-raphson technique
for i in h:
	# Assuming initial pressure value
	p_guess = 100
	# Defining Tolerance
	tolerance = 1e-6
	# Loop counter for incrementing the number of iterations
	counter = 1
	while(abs(f(p_guess, i, r, beta, sigma)) > tolerance):
		p_guess = p_guess - alpha*(f(p_guess, i, r, beta, sigma)/fprime(p_guess, i, r, beta, sigma))
		counter = counter+1
	pressure.append(p_guess)
	pressure_0.append(p_guess/144)
	num.append(counter)

print('Thickness of the Ice wedge :',h)

print('Cushion Pressure (pound per square feet) :',pressure)
print('Cushion Pressure (psi) :',pressure_0)
print('Number of iterations :',num)
# Plotting graph
plt.figure(1)
plt.plot(h, pressure, marker = '.',color='blue',markerfacecolor='Yellow',markersize=10)
plt.title('Cushion pressure versus Ice field thickness')
plt.xlabel('Thickness of ice field (feet)')
plt.ylabel('Cushion Pressure (pound per square feet)')
plt.grid()
plt.show()

plt.figure(2)
plt.plot(h, pressure_0, marker = '.',color='blue',markerfacecolor='red',markersize=10)
plt.title('Cushion pressure versus Ice field thickness')
plt.xlabel('Thickness of ice field (feet)')
plt.ylabel('Cushion Pressure (pound per square inches)')
plt.grid()
plt.show()

# Part 2 of the program to calculate optimal relaxation factor for given thickness of ice field

# h_1 the thickness of the ice field
h_1 = 0.6
# Relaxation factor array assumed considering underrelaxing and overrelaxing values
alpha_array = np.linspace(0.1,1.9,50)
# Null array for storing the number of iterations in future
iterations = []
# Creating a for loop for interation and performing the newtion-raphson technique
for j in alpha_array:
	# Assuming initial pressure value
	p_guess = 100
	# Tolerance
	tolerance = 1e-6
	# Loop counter for incrementing the number of iterations
	counter = 1
	while (abs(f(p_guess,h_1,r,beta,sigma)) > tolerance):
		p_guess = p_guess - (j*f(p_guess,h_1,r,beta,sigma))/fprime(p_guess,h_1,r,beta,sigma)
		counter = counter + 1
	iterations.append(counter)

print(":"" alpha"" : ""No. of iterations")

for j in range(0,len(alpha_array)):
	print(alpha_array[j],iterations[j])
for k in range(len(iterations)):
	if iterations[k] == min(iterations):
		print('The optimum value of under Relaxation Factor: ', str(alpha_array[k]))
# Plotting graph
plt.figure(3)
plt.plot(alpha_array,iterations,marker = '.',color='black',markerfacecolor='orange',markersize=7.5)
plt.xlabel('Relaxation factor (alpha)')
plt.ylabel('No. of iterations')
plt.title('Iterations versus Relaxation factor')
plt.grid()
plt.show()

EXPLANATION OF THE PROGRAM:

1. At first necessary modules such as math, matplotlib, numpy are imported and required inputs are given. For the first part of code an array of and an assumed value of is considered according to the challenge statement.

2. Two different functions named and are defined according to the formula given, and an empty array of number of interations and pressure(in psi and psft) are created.

3. For calculating minimum pressure for different values of ice thickness a for loop is executed taking an initial guessed value

of pressure, tolerance, and a loop counter =1 which will increment after every iterations. Then a while loop is initiated inside

the same for loop for checking the whether all inputs in the formula if is greater than tolerance provided, a new value of pressure will be calculated according to Newton Raphson method.

4. The calculated values of pressure and number of iterations are stored in their respective arrays using 'append' command

and necessary statements are printed in command window.

5. Graphs for pressure versus ice thickness is plotted using calculated data set.

6. For calculating optimal relaxation factor( ) for given thickness of ice same steps as above is repeated.

7. Number of iterations calculated is stored as an array and a graph is plotted between the same and relaxation factors.

RESULTS

output:

1. In the output window we can see the cushion pressure for the respective heights of the ice wedge.

2. Variation of cushion pressure with change in thickness of ice field.

3. No of iterations required for different relaxation factors.

It can be concluded that curve rapidly decreases as number of iterations decreases and relaxation factor increases

optimal relaxation factor lies in the vicinity of 1. By looking at this graph we can conclude that optimal relaxation factor for this problem is 1.01836 .

Limitations of Newton Raphson method:

1. If f(x) = 0, this method will not work.
2. If f''(x) = 0, this causes the solution to diverge.
3. Sometimes get trapped with local minima and maxima.

Conclusion:

-The plots and results show that by using a suitable relaxation factor we can obtain results with as much accuracy, in a shorter processing time.

-Newton-Raphson method is simplest method of calculating roots for Non-linear problems. Optimum choice under Relaxation factor helps reducing calculation time by reducing no. of iterations required to calculate the root. Also to restrict minimum value root a tolerance value must be chosen for simplicity in resolving problems.