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/A₀) (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 Stress	Elastic 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.

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

Meshing

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

Creating Section

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

Assigning Section

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

Creating Partition

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

Setting Material Properties

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

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.

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.

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.

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

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

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.

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.

Showing relationships between moment and arc length, and rotation and arc length 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.

From the following comparison, we see that the buckling is higher when there is no pressure than the buckling with internal 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 = x³ + x² – 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 = x³ + x² – 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)

Isoparametric Full Integration Element Model Simulation by Mathematica and Abaqus

In this tutorial, a quadrilateral 4 node isoparametric full integration element with coordinates as node-1: (0, 0), node-2: (5, 0), node-3: (5, 5), and node-4: (0, 5) is considered. Let's assume, Young's modulus, E = 20 GPa, and Poisson's ratio, ν = 0.2. We are interested here to compare the stiffness matrices for the plane strain condition (having unit thickness), which are determined by Mathematica as well as a finite element software, ABAQUS. Let's consider that the node-1 and -2 are fixed, and we are required to calculate the displacement of the other two nodes (3 and 4) under gravitational load. The material density is 2400 kg/m³. We will also find out the stresses and strains at the integration points of the element. 

The following Mathematica program is developed to determine the stiffness matrix, stress, strain and displacement manually (with Mathematica) of the nodal points of a quadrilateral 4 node isoparametric full integration element.

Manual Approach (Mathematica):

(*Mathematica Codes to the Solve Problem*)

(*Definition of Nodal Coordinates from Problem 1*)

coordinates = {{0, 0}, {5, 0}, {5, 5}, {0, 5}};

(*Properties in Problem 1*)

t = 1;

Ee = 20*10^9;

Nu = 2/10;

Cc = Ee/((1 + Nu)*(1 - 2 Nu))*{{1 - Nu, Nu, 0}, {Nu, 1 - Nu, 0}, {0,

     0, 1/2 (1 - 2 Nu)}};

(*Defining the Shape Functions for Each Node*)

Shapefun = Table[0, {i, 1, 4}];

Shapefun[[1]] = (1 - xi) (1 - eta)/4;

Shapefun[[2]] = (1 + xi) (1 - eta)/4;

Shapefun[[3]] = (1 + xi) (1 + eta)/4;

Shapefun[[4]] = (1 - xi) (1 + eta)/4;

Nn = Table[0, {i, 1, 2}, {j, 1, 8}];

Do[Nn[[1, 2 i - 1]] = Nn[[2, 2 i]] = Shapefun[[i]], {i, 1, 4}];

x = FullSimplify[Sum[Shapefun[[i]]*coordinates[[i, 1]], {i, 1, 4}]];

y = FullSimplify[Sum[Shapefun[[i]]*coordinates[[i, 2]], {i, 1, 4}]];

J = Table[0, {i, 1, 2}, {j, 1, 2}];

J[[1, 1]] = D[x, xi];

J[[1, 2]] = D[x, eta];

J[[2, 1]] = D[y, xi];

J[[2, 2]] = D[y, eta];

JinvT = Transpose[Inverse[J]];

mat1 = {{1, 0, 0, 0}, {0, 0, 0, 1}, {0, 1, 1, 0}};

mat2 = {{JinvT[[1, 1]], JinvT[[1, 2]], 0, 0}, {JinvT[[2, 1]],

    JinvT[[2, 2]], 0, 0}, {0, 0, JinvT[[1, 1]], JinvT[[1, 2]]}, {0, 0,

     JinvT[[2, 1]], JinvT[[2, 2]]}};

mat3 = Table[0, {i, 1, 4}, {j, 1, 8}];

Do[mat3[[1, 2 i - 1]] = mat3[[3, 2 i]] = D[Shapefun[[i]], xi];

  mat3[[2, 2 i - 1]] = mat3[[4, 2 i]] = D[Shapefun[[i]], eta], {i, 1,

   4}];

B = Chop[FullSimplify[mat1.mat2.mat3]];

Kbeforeintegration = FullSimplify[Transpose[B].Cc.B];

xiset = {-1/Sqrt[3], 1/Sqrt[3]};

etaset = {-1/Sqrt[3], 1/Sqrt[3]};

Kk = t*FullSimplify[

    Sum[(Kbeforeintegration*Det[J] /. {xi -> xiset[[i]],

        eta -> etaset[[j]]}), {i, 1, 2}, {j, 1, 2}]];

Kk = Round[Kk, 0.0001];

Print["The Stiffness Matrix is,"]

ScientificForm[Kk // MatrixForm, 4]

Print[""]

(*Determining the Force Vector*)

rb = {0, -9.81*2400};

fbeforeintegration = t*Transpose[Nn].rb*Det[J];

Ff = FullSimplify[

   Sum[(fbeforeintegration /. {xi -> xiset[[i]],

       eta -> etaset[[j]]}), {i, 1, 2}, {j, 1, 2}]];

Print["The Force Vector is,"]

FF = Transpose[

   Join[{{"Fx_1 =", "Fy_1 =", "Fx_2 =", "Fy_2 =", "Fx_3 =", "Fy_3 =",

      "Fx_4 =", "Fy_4 ="}}, {Ff}]];

ScientificForm[FF // MatrixForm, 4]

Print[""] 

(*Determining the Displacement*)

ID = {5, 6, 7, 8};

Kkr = Kk[[ID, ID]];

Ffr = Ff[[ID]];

Uur = Inverse[Kkr].Ffr;

Uu = Table[0, {i, 1, 8}];

Uu[[ID]] = Uur;

Print["The Displacement Vector is,"]

UU = Transpose[

   Join[{{"u_1 =", "v_1 =", "u_2 =", "v_2 =", "u_3 =", "v_3 =",

      "u_4 =", "v_4 ="}}, {Uu}]];

ScientificForm[UU // MatrixForm, 4]

Print[""] 

(*Determining the Stress and Strain*)

Eps = B.Uu;

Eps[[3]] = Eps[[3]]/2;

Sig = Cc.Eps;

EPS = ConstantArray[0, {4, 8}];

Do[EPS[[i, 1]] = i, {i, 1, 4}];

ii = 1;

Do[EPS[[ii, 1]] = xiset[[i]];

 EPS[[ii, 2]] = etaset[[j]];

 EPS[[ii, {3, 4, 5}]] = Eps /. {xi -> xiset[[i]], eta -> etaset[[j]]};

 EPS[[ii, {6, 7, 8}]] = Sig /. {xi -> xiset[[i]], eta -> etaset[[j]]};

 ii++, {i, 1, 2}, {j, 1, 2}]

Print["Stress and Strain Components are,"]

EPS = Join[{{"x_i", "eta", "e_xx", "e_yy", "e_xy", "s_xx", "s_yy",

     "s_xy"}}, EPS];

EPS1 = EPS;

ScientificForm[EPS // MatrixForm, 4]

Print[""]

 

Program Outputs:

Now, a model of a quadrilateral 4 node isoparametric full integration element is developed in Abaqus to determine the same stiffness matrix, stress, strain and displacement of the nodal points, which were obtained previously manually.

Using Abaqus:
In Abaqus, the geometry is defined by using the given nodal points in assignment. The material is elastic, where Young’s modulus is 20000000000 pa, and Poisson’s ratio is 0.2. In meshing, only a single element is considered with full integration option. The following two figures show the variations of displacements in magnitude.

The following figure shows the variations of stress in magnitude.

The Output File for the Stresses:

The Output File for the Strains:

Comparison between Stiffness Matrices:

Mathematica Calculation

From Abaqus

By comparing all the results (stress, strain, and displacement) between manual calculation (Mathematica) and Abaqus model, we see that except the stress values, all other parameters are close. However, the stresses differ in the Abaqus model.