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


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.


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

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:
		x_guess = x_guess-alpha*(f(x_guess,i)/fprime(x_guess,i))
		iter = iter+1

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

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


tol = 1e-4

for i in alpha:
	x_guess = 10				
		x_guess = x_guess-i*(f(x_guess,h[0])/fprime(x_guess,h[0]))
		iter = iter+1

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

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


The Table we get in the result is
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 :

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

                                        Relaxation factor vs Iteration


  • 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.

Projects by Jitesh Sahjwani

OBJECTIVE In this challenge, we check and create the properties for a Neon side crash-BIW model and run the simulation & observe results. Question: Neon side crash -BIW Check the unit system and either follow [Mg mm s] or [Kg mm ms]. Create an appropriate int Read more

OBJECTIVE In this challenge, we check and create the properties for a Neon frontal crash-BIW model and run the simulation & observe results. Question: Frontal crash-BIW Check the unit system and either follow[Mg mm s] or [Kg mm ms]. Create an appropriate inte Read more

OBJECTIVE In this report, we mesh a Bumper assembly.Perform simulation on Crash tube and compare the Type-7 simulation results, Type-11 along with Type-7 and different kinematic condition simulation results. Create six different cases  Run the crash tube model Read more

OBJECTIVE The main objective is to analyze and compare the failure behavior of a plate using different material models LAW-1, LAW-2, LAW-27 & LAW-36 along with Johnson-cook Failure card. For LAW-2, Once use Johnson\'s Failure card along with some certain p Read more

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

Google drive link for Hypermesh files:   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

Google drive link for the \'Engine_data.out\' file:   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

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

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


The End