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

OBJECTIVE

The key objective is to determine the minimum cushion pressure needed to break a given thickness of the ice using an air-cushion vehicle. Solve the equation by using the Newton Raphson method and Write a Python code.

Methodology that we use;

Whole idea of newton Rapson method is to Linearize the problem.

(“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) convert this unit into the pound per sq. feet.

Construct a table with pressure (p) in the first column, no. of iteration in the second column and thickness(h) in the third column for h = 0.6, 1.2, 1.8, 2.4, 3.0, 3.6 and 4.2 feet.

Use relaxation factor of 1 and explain your results and choose the optimum relaxation factor for this case. Show with the help of a graph, assume h = 0.6 and at last, we Explain results with respect to the appropriate relaxation factor.

PROCEDURE

1. Split the equation into four parts to make the calculation easy.

2. Import all the necessary modules;

• import matplotlib.pyplot as plt
This module helps us to plot the graph of the given result.
• import math
This module helps us to use mathematical functions like pie, sine, cos, etc.
• from prettytable import PrettyTable
This is very helpful in making tables for our given result.
• import numpy
This module helps us to do mathematical operations on an element of an array.
Whenever we want to define more than one value in a single variable we have to use an array and numpy.linspace is also the same thing it stores value in a range like this; numpy.linspace (lower limit, upper limit, no. of division)
This works acc. to a formula:
Size of division = (upper limit-lower limit)/(no. of division-1)
• From scipy.integrate import odeint
Odeint function which integrates the system of ODE.

3. Use an empty array for storing a value that we get from FOR loop.

4. Append is a command that will continuously append the value of a variable. It keeps on adding the values into an existing array.

5. Define the function: We can call the function whenever and wherever we required.
def fprime(p,h): is the derivative function of def f(p,h):.

6. x_guess is the first guess value of pressure which we use assume for the iteration process in the N-R method and we assume tolerance tol= 1e-4 and relaxation factor is alpha = 1.

7. Performing N-R iteration

we make an empty array it= [] that stores no. of iterations value.

we use FOR loop for using the value of h in that.
we use While loop for doing iteration up till pressure value we get is less than the tolerance.

x_(nguess) = x_(guess)-alpha(f(x_(guess))/(f'(x_(guess))))

8. According to our second question requirement: we use numpy module for finding optimum relaxation factor. we take 100 values of alpha in between the 0.1 to 1.9.

We use again k = [] is an empty array which stores an appended values of iteration.

"""
NR iteration method
By Jitesh
"""

#Inputs
import matplotlib.pyplot as plt
import math
from prettytable import PrettyTable
import numpy as np

#assign the values
beta=0.5        #width of ice
r = 40           #feet
sigma= 21600     #sigma = 150 *144 for pound #psi to feet
h = [0.6, 1.2, 1.8, 2.4, 3.0, 3.6, 4.2] #feet
p = []

def f(p,h):
x1 = pow(p,3)*(1-pow(beta,2))
x2 = (0.4*h*pow(beta,2) -(sigma*pow(h,2)/pow(r,2)))*p**2
x3 = pow(sigma,2)*pow(h,4)*p/(3*pow(r,4))
x4 = pow((sigma*h*h/(3*r*r)),3)
f1= x1+x2+x3-x4
return f1

def fprime(p,h):
y1 = 3*pow(p,2)*(1-beta*beta)
y2 = 2*p*((0.4*h*beta*beta)-(sigma*h*h/(r*r)))
y3 = pow(sigma,2)*pow(h,4)/(3*pow(r,4))
fprime1 = y1+y2+y3
return fprime1

#Performing N-R iteration
alpha = 1   #relaxtion factor
tol = 1e-4
x_guess = 10
it=[]  #stores no. of iteration value

for i in h:
iter=1
while(abs(f(x_guess,i))>tol):
x_guess = x_guess-alpha*(f(x_guess,i)/fprime(x_guess,i))
iter = iter+1
print(x_guess,iter,h)
p.append(x_guess)
it.append(iter)

table= PrettyTable(['x_guess','iter','h'])
for j in range (0,7):
print (table)

alpha = np.linspace(0.1,1.9,100)

k=[]

tol = 1e-4

for i in alpha:
iter=1
x_guess = 10
while(abs(f(x_guess,h[0]))>tol):
x_guess = x_guess-i*(f(x_guess,h[0])/fprime(x_guess,h[0]))
iter = iter+1
print(x_guess,iter,h)
k.append(iter)

print(k)
print('the number of iterations('+str(min(k))+') is least at :')
for i in range(len(k)):

if k [i]==min(k):
print(alpha[i])

plt.plot(alpha,k)
plt.show()

The Table we get in the result is
pressure
in the first column, no. of iteration in the second column and thickness in the third column.

+--------------------+------+-----+
|      x_guess       | iter |  h  |
+--------------------+------+-----+
| 4.239063886378783  |  8   | 0.6 |
| 17.23669061548758  |  7   | 1.2 |
| 38.99008116755394  |  24  | 1.8 |
|  69.4992804690224  |  18  | 2.4 |
| 108.76429695856204 |  40  | 3.0 |
| 156.7851333779438  |  21  | 3.6 |
| 213.56179088560847 |  60  | 4.2 |
+--------------------+------+-----+

The minimum no. of Iteration is 8 between the alpha value of 0.999 to 1.154.

the number of iterations(8) is least at :
0.9909090909090909
1.009090909090909
1.0272727272727273
1.0454545454545454
1.0636363636363637
1.1545454545454545

We get the graph after running the code between the relaxation factor and Iteration.

Relaxation factor vs Iteration

CONCLUSION

• Here we learn how to solve the no. of equations using program.
• we get to know that optimum relaxation factor (alpha) is 1.
• We get to know that the minimum no. of iteration required is 8.
• We learned how to draw a table in our result and how to plot the graph
• Understand the use of FOR loop and While loop.

Comparison of Simulation results and Improve Shell Property of rail model using Hypermesh RADIOSS Jitesh Sahjwani · 2019-09-20 14:57:35

OBJECTIVE Using the crash beam file, Comparison of simulation results of the base setup and improved shell element properties. Questions: 1.Using the crash beam file and change the run time to 55 ms. 2. Change the number of animation steps during simulation to a min Read more

Crash Worthiness Analysis of Housing and Arm Bracket using Hypermesh Radioss Jitesh Sahjwani · 2019-09-20 04:21:10

Crash Worthiness analysis of housing and arm bracket using Radioss Hypermesh OBJECTIVE The key objective is to create the 2D and 3D meshing for the given models and solve the following questions PROCEDURE & Answers 1. Material property for steel we can ch Read more

Geometry Cleanup and Mesh on 2D SHELL Elements Using Hypermesh Radioss Jitesh Sahjwani · 2019-09-17 07:48:47

Google drive link for Hypermesh files: https://drive.google.com/open?id=1MDhngswNUUJFrXmWwMTR4jfR0R9T6PKd   OBJECTIVE The key objective is to carry out the geometry cleanup operation & ignore holes less than dia 5mm, extract mid surface and create a 2D mesh Read more

Data analysis File parsing script using Python Jitesh Sahjwani · 2019-09-14 04:40:04

Google drive link for the \'Engine_data.out\' file: https://drive.google.com/open?id=1lIh6RdUZRex9klgtDhpQ0tmtRTBY11YC   OBJECTIVE The primary aim is to create a data visualizer tool using python and perform certain things by this code is: The scri Read more

Curve Fitting for Temperature and specific heat values using Python Jitesh Sahjwani · 2019-09-13 18:38:15

OBJECTIVE The key objective of writing this program is to perform curve fitting with the help of a linear and cubic polynomial for the Temp. & Cp data. #1. What does popt and pcov mean? #2. What does np.array(temperature) do? #3. What does the * in *popt Read more

Constraint Minimization of Equation using lagrange's method Jitesh Sahjwani · 2019-09-11 06:59:11

OBJECTIVE The main objective is to Minimize the given equation with respect to the given constraint. Problem: This function f(x,y) = 5-(x-2)^2 -2(y-1)^2 is subject to the following constraint x + 4y = 3 Solution: Make a new function as g(x,y) = Read more

Transient behavior of a simple pendulum using python Jitesh Sahjwani · 2019-09-10 15:38:28

OBJECTIVE  The key objective is to write a program that solves the following ODE and simulate & create an animation of the transient behavior of a simple pendulum using python. 2nd Order Equation given (d^2theta)/dt^2 + b/m (d theta)/dt +g/l sin theta = 0` Read more

Simulation of the forward kinematics for 2R Robotic Arm Using Python Jitesh Sahjwani · 2019-09-10 04:21:22

OBJECTIVE The main aim of this program in python to simulate the forward kinematics of a 2R Robotic Arm and also make animations of the plotted graph. 2R robotic arm contains two links and two rotational joints. Workspace refers to the possible locations that th Read more

Simulation of Otto Cycle with P-V curve plotting using Python Jitesh Sahjwani · 2019-09-10 04:09:38

OBJECTIVE The key objective of this program is to make OTTO cycle simulator.Calculate the unknown pressure & temperature and plot the PV diagram using thermal equations and find the thermal efficiency with the help of Python scripting. Otto cycle: The air-standard Read more