Welcome to the World of Modelling and Simulation

What is Modelling?

This blog is all about system dynamics modelling and simulation applied in the engineering field, especially mechanical, electrical, and ...

Implementation of Newton Raphson and Modified Newton Raphson Method by Mathematica

This article discusses the Mathematica approaches to code Newton Raphson (NR) and Modified Newton Raphson (MNR) and the differences if any in computation. Let's assume that stress and strain of a rigid bar is defined by the following relationship:

𝜎 = 100000 (√𝜀 + (1/80) sin (𝜋𝜀/0.002))

Now, taking the first derivative of the above relationship, we may write,

𝜎' = 100000 ((1/2√𝜀) + 19.653 cos (1570.80𝜀))

The NR method implements an iterative approach by considering a tangent to a curve, for which the function under consideration becomes zero. If f(xⱼ) is any function and f′(xⱼ) is the derivative of that function at point xⱼ, then according to NR approach, we may write as,

xⱼ₊₁ = xⱼ - f(xⱼ) / f′(xⱼ)

The MNR works in a similar way, but is helpful for converging a solution efficiently for roots with multiplicity and may be represented as,

xⱼ₊₁ = xⱼ - ⟮ f(xⱼ)  f′(xⱼ)⟯ /⟮ f′(xⱼ)² -  f(xⱼ) f′′(xⱼ)⟯



The following Mathematica program is developed to plot the function and its derivatives.

(*Function Defined in the Problem*)
sigma = 100000*(Sqrt[epsilon] + 1/80*Sin[Pi*epsilon/0.002]);
y1 = 4000;
y2 = 6000;
sigmap = D[sigma, epsilon];
sigmap2 = D[sigma, {epsilon, 2}];

(*Plotting the Function, its Derivatives, Upper and Lower Bounds*)
A1 = Plot[{sigma, y1, y2}, {epsilon, 0, 0.005}, ImageSize -> Large, 
  AxesLabel -> {"Epsilon", "Sigma"}, PlotStyle -> {Red, Green, Black, Thick}]
A2 = Plot[sigmap, {epsilon, 0, 0.005}, ImageSize -> Large, 
  AxesLabel -> {"Epsilon", "Sigma'"}, PlotStyle -> {Red, Thick}]
A3 = Plot[sigmap2, {epsilon, 0, 0.005}, ImageSize -> Large, 
  AxesLabel -> {"Epsilon", "Sigma''"}, PlotStyle -> {Red, Thick}]

(*Initial Guesses*)
sigma /. epsilon -> 7.928*10^-4
sigma /. epsilon -> 7.924*10^-4
sigma /. epsilon -> 3.882*10^-3
sigma /. epsilon -> 3.879*10^-3


plotting the function

plotting the function

plotting the function

function results














Mathematica Program for the NR Method:


(*Mathematica Codes for Newton-Raphson Approach*)
Sigma = 100000*(Sqrt[epsilon] + 1/80*Sin[Pi*epsilon/0.002]) - 4000;
SigmaDerivative = D[Sigma, epsilon];
epsilontable = {0.0002};
Error = {1};
SetTolerance = 0.0005;
MaximumIteration = 100;
i = 1;

(*Newton-Raphson Algorithm Implementation*)
While[And[i <= MaximumIteration, Abs[Error[[i]]] > SetTolerance],
  epsilonnew = 
   epsilontable[[i]] - (Sigma /. 
       epsilon -> epsilontable[[i]])/(SigmaDerivative /. 
       epsilon -> epsilontable[[i]]); 
  epsilontable = Append[epsilontable, epsilonnew]; 
  Errornew = (epsilonnew - epsilontable[[i]])/epsilonnew; 
  Error = Append[Error, Errornew];
  i++];

L = Length[epsilontable];
SolutionTable = 
  Table[{i - 1, epsilontable[[i]], Error[[i]]}, {i, 1, L}];
SolutionTable1 = {"Iteration Number", "Epsilon", "Error"};
L = Prepend[SolutionTable, SolutionTable1];
Print["Showing Results for the Newton-Raphson Approach when Sigma is \
4,000 and Initial Guess is 0.0002"]
ScientificForm[L // MatrixForm, 4]



Mathematica Program for the MNR Method:

(*Mathematica Codes for the Modified Newton-Raphson Approach*)
Sigma = 100000*(Sqrt[epsilon] + 1/80*Sin[Pi*epsilon/0.002]) - 4000;
epsilontable = {0.0002};
SigmaDerivative = D[Sigma, epsilon] /. epsilon -> epsilontable[[1]];
Error = {1};
SetTolerance = 0.0005;
MaximumIteration = 100;
i = 1;

(*Implementation of the Modified Newton-Raphson Algorithm*)
While[And[i <= MaximumIteration, Abs[Error[[i]]] > SetTolerance],
  epsilonnew = 
   epsilontable[[i]] - (Sigma /. epsilon -> epsilontable[[i]])/
     SigmaDerivative; epsilontable = Append[epsilontable, epsilonnew];
   Errornew = (epsilonnew - epsilontable[[i]])/epsilonnew; 
  Error = Append[Error, Errornew];
  i++];

T = Length[epsilontable];
SolutionTable = 
  Table[{i - 1, epsilontable[[i]], Error[[i]]}, {i, 1, T}];
SolutionTable1 = {"Iteration Number", "Epsilon", "Error"};
T = Prepend[SolutionTable, SolutionTable1];
Print["Showing Results for the Modified Newton-Raphson Approach when \
Sigma is 4,000 and Initial Guess is 0.0002"]
ScientificForm[T // MatrixForm, 4]


Results:

The following table shows the comparison of results from NR and MNR methods. We see that at iteration number 17, the Modified Newton-Raphson approach converges, whereas the Newton-Raphson approach is still iterating. However, it is noticeable that the step is faster than the Modified approach. Nevertheless, Newton-Raphson approach fails to converge early while the Modified approach, although iterates slowly, but converges eventually.

Comparison between Newton-Raphson and Modified Newton-Raphson Approach























#NewtonRaphson #Mathematica #MatrixForm #Blog #Blogger

Discussion on MATLAB Function and Function Handle

In MATLAB, the functions and function handles (@) are used frequently, especially when writing codes in m-files. We can create many functions, but when we call a particular function, and if that function has another function inside the input parameters - the function handle (@) is the way to address and call the second function. For example, we have created a function to find the roots of a polynomial using the Newton-Raphson method and named it as 'newtonraphson'. The function input parameters are: polynomial, initial guess, error tolerance and maximum number of iterations. And the output is the roots of the polynomial.

Now, to define the polynomial, we can create another function and name it as 'polyfunction'. This function only has the polynomial, and we are interested to find its roots. Let's see how we can call the function 'newtonraphson' in MATLAB command prompt.

>> newtonraphson(@ polyfunction, 1.0, 1e-5, 100)
 
Here,

newtonraphson = a function for finding roots of a polynomical
@ = function handle to call another function
1.0 = Initial guess
1e-5 = Error tolerance
100 = Maximum number of iterations 

If you do not put the function handle (@) in front of the 'polyfunction', there will be an error you will see instantly in the command prompt. As you can see, a function handle is used to pass information to another function. There are other uses of the function handle in MATLAB, but the above is most frequently and widely used expression.

The aforementioned example shows the application of the function handle (@) to a known or defined function by the user (e.g. 'polyfunction'). Now, if there is an anonymous equation that we need to define in our coding, then we can also use the function handle to do the job. For example, in our program, we need to define an equation: 3x - 4y + 2z. The x, y, and z are the unknowns. To define it with the function handle, we can write as,

>> L = @ (x, y, z) (3*x - 4*y + 2*z)

Here,

L = is any parameter referring the equation
@ = function handle defining x, y, z with respect to the equation for evaluation



#MatlabFunction #FunctionHandle #PolyFunction #NewtonRaphson #Matlab #Blog #Blogger

Simulation of Buckling Behavior of a Pipe by Abaqus

In this tutorial, I am going to show how we can simulate the buckling characteristics of a pipe in Abaqus. The length of the pipe is 4.5 m, outside diameter (OD) is 0.9 m, and thickness is 9 mm. First, let us consider the stress-strain relationship, which is defined by the following table.


Table 1: Stress-Strain Relation for a Pipe

True Strain     True Stress (N/mm^2)
0                              0
0.001                     200
0.002                     400
0.003                     600
0.004                     800
0.0055                     900
0.008                     1000
0.0125                     1100
0.026                     1200


The relationship between the true and engineering stress-strain is governed by:

σtrue = F/A = (F/A0) (1+εeng) = σeng(1+εeng)


Now, to find the elastic strain, the following relation is used.

εe = σ/E


The following table is based on the above definitions.


Table 2: Elastic and Plastic Strain for a Pipe

True
Strain
True Stress
(N/mm^2)
Tr StressElastic
Strain (εe)
Plastic Strain (εp)
 0 0 0 0 0
 0.001 200 200.2 0.001001 -1e-06
 0.002 400 400.8 0.002004 -4e-06
 0.003 600 601.8 0.003009 -9e-06
 0.004 800 803.2 0.004016 -1.6e-05
 0.0055 900 904.95 0.004525 0.000975
 0.008 1000 1008 0.00504 0.00296
 0.0125 1100 1113.75 0.005569 0.006931
 0.026 1200 1231.2 0.006156 0.019844


The internal pressure inside the pipe is given by following relation:

P = [0.8(2t)σy]/OD

P = 12.8 MPa


This is the ultimate internal pressure of the pipe. Now, we are ready to create the model in Abaqus.

Steps in Abaqus Modelling:

Creating the Geometry
First, a part is created in Abaqus CAE, where a 3-point arc is chosen to draw a half circle with radius 450 mm. Since, the part is selected as 3D shell with extrusion, the length is 2250 mm. We will use symmetry and that’s why, these dimensions are selected to use symmetry boundary conditions about X and Z-axes.


3D shell with extrusion


The S4R shell element is chosen with reduced integration for analysis.

3D shell with extrusion


Meshing

Next, we create the mesh. The mesh size is set to 18 mm, which is twice the size of the shell thickness.


Meshing

Meshing


Creating Section
After that, the sections are created. The shell thickness is set 9 mm as per the problem.

Creating Section






















Creating Section






















Assigning Section
After creating the section, we assign it to the part. Then, we use “Render Shell Thickness” to show the shell thickness.

Assigning Section

Assigning Section


Creating Partition
A partition is created with 200 mm offset from one of the pipe edges.

Creating Partition

Setting Material Properties

For the elastoplastic behaviour, the following material properties are given:


Setting Material Properties


Creating Steps with Loads and Boundary Conditions

In this step, we define the boundary condition (BC) and load. A reference point (RP) is created to fix one end of the pipe. For the RP, a rigid body constraint is defined. Three BCs are generated as follow:

BC-1: for symmetry about X-axis

BC-2: for symmetry about Z-axis

BC-3: for the reference point (RP1), which is fixed along the axes.

Then, a pressure load of magnitude 12.8 MPa is applied inside the pipe.

Creating Steps with Loads and Boundary Conditions

Results and Discussions: 

Since the model is created with necessary boundary conditions and loads, in this section, we are interested to consider two circumstances.

Buckling when both pressure and bending are considered

Buckling when bending is considered without pressure

We apply the static Rik’s approach as a force control. The following figure shows the Von Mises stress when only pressure is applied inside the wall of the pipe.

Von Mises stress when pressure is applied

Buckling when both pressure and bending are considered

At first, we apply a moment of 5x10^10 Nmm at the reference point by creating a new step where Rik’s static algorithm is selected. In this scenario, the pressure is still present inside the pipe. The following figure demonstrates the three steps for the simulation.

Showing 3 steps: no pressure (left), with pressure (center), and with pressure and rotation (right)

The following figure shows two relationship when pressure is present: moment vs arc length and rotation vs arc length. The intersection point between two graphs in the figure is the convergence point for the simulation.

Showing relationships between moment and arc length, and rotation and arc length with pressure

Showing relationships between moment and arc length, and rotation and arc length with pressure


Buckling when both pressure and bending are considered

The following plot shows the relation between moment and angular displacement when pressure is present in the pipe.

relation between moment and angular displacement





















Buckling when bending is considered without pressure

In this step, we just suppress the load leaving all other conditions same. We are interested to see the effect of the bending only on the buckling rate of the pipe.

Buckling when bending is considered without pressure
Showing 2 steps: no pressure (left), and with angular displacement (right)



The following figure shows two relationship when pressure is absent: moment vs arc length and rotation vs arc length. The intersection point between two graphs in the figure is the convergence point for the simulation.



relationships between moment and arc length
Showing relationships between moment and arc length, and rotation and arc length without pressure
























Buckling when bending is considered without pressure
Moment without pressure (left) and rotation without pressure (right) when arc length is 0.8016 mm













The following plot shows the relation between moment and angular displacement when pressure is absent in the pipe.


relation between moment and angular displacement


From the following comparison, we see that the buckling is higher when there is no pressure than the buckling with internal pressure.


buckling is higher when there is no pressure


To learn in detail of FEA, you may visit here.



#VonMisesstress #Buckling #Bending #FiniteElementAnalysis #FEA #Abaqus #ModellingSimulation #Blog #Blogger

A MATLAB Program to Find the Roots of a Function by Bisection Method

In this article, a MATLAB program is developed to find the roots of any function by Bisection approach. There are two functions created for this purpose. The first one implements the Bisection method. And the second one is for the function itself, for which, we are interested to find the roots. Let us assume a function,

y = x3 + x2 – 5x – 3 

First Function: Bisection Method

function [root,f_x,e_r,iter] = bisectionmethod(func_name,x_lower,x_upper,e_tol,max_iter,varargin)

% INPUTS:

% func_name: input function to find its roots

% x_lower: lower limit of the bracket

% x_upper: upper limit of the bracket

% e_tol: defining error tolerance for this method

% max_iter: Total iteration number

% varargin: to take any number of inputs for the function

 

% OUTPUTS:

% root: root of the given function

% f_x: function value at the root

% e_r: relative error

% iter: number of iteration taken


close all;

clc;

if nargin < 3  % Defines the number of input function arguments

    error('Minimum 3 Input Arguments are Required to Run the Function')

end

% Product of the function at upper & lower interval

product_of_functions = func_name(x_lower,varargin{:}) * func_name(x_upper,varargin{:});

 

if product_of_functions > 0

    error('No Roots Found within the Given Range'),

end


% Logical operator || 'or'

if nargin <4 || isempty(e_tol)

    e_tol=0.0001;

end 

if nargin <5 || isempty(max_iter)

    max_iter=50;

end


% Assuming some parameters to initiate 'while' loop

iter = 0; xr = x_lower; e_r = 1;


while(1)

    xrold=xr;

    xr=(x_lower + x_upper)/2; % Bi-section Approach

    iter = iter + 1;


    if xr~=0  % Logical operation ~= 'not equal'

       e_r=abs((xr-xrold)/xr);

    end

    product_of_functions = func_name(x_lower,varargin{:}) * func_name(xr,varargin{:});

    if product_of_functions < 0

        x_upper = xr;

    elseif product_of_functions > 0

        x_lower = xr;      

    else

        e_r=0;

    end

    if e_r <= e_tol || iter >= max_iter

        break,

    end

end

 root = xr;

 f_x = func_name(xr, varargin{:});

end


Second Function: The Polynomial (y = x3 + x2 – 5x – 3) 

function y = func_example(x)
% Defining a function to find its roots
y = x.^3 + x.^2 - 5.*x - 3; 
end 


Now, when you are going to call the Bisection Method from the MATLAB command prompt, use the following syntax where the function handle is used to call another function, which is the polynomial in our case.

>> bisectionmethod(@func_example,-100,100)