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

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

 TrueStrain 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-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 = 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)

### 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/m3. 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.