Vectorization in MATLAB: Use of Colon (:)

We are familiar with 'for loop' application in MATLAB. We may use many 'loops' to handle complex calculations. But, sometimes too many 'loops' cost computational time, which eventually affects the overall efficiency of the codings. In MATLAB, instead of looping, we may use the method of vectorization to skip some 'loops'. We may use the operator colon (:) in MATLAB to form vectors that alternately eliminates the 'loops'. Let's explain the method with a very simple example. For example, we would like to divide a square matrix (for our example, 5x5) with a scalar number. Below, there are three scenarios. At first, we can use 2 'for loops' to do the operation. Secondly, we can skip one 'for loop' by using a colon (:) for the columns of the matrix a. Finally, we can avoid all 'for loops', and use colons (:) for both rows and columns of the example matrix a.

2 'for loop':

i=1;

j=1;

n=5;

a=[1 2 3 4 5; 6 7 8 9 10; 1 2 3 4 5; 6 7 8 9 10; 1 2 3 4 5;];

 for i=1:n

     for j=1:n

     a(i,j) = a(i,j)./n;

     end

 end

 disp(a)

1 'for loop':

i=1;

j=1;

n=5;

a=[1 2 3 4 5; 6 7 8 9 10; 1 2 3 4 5; 6 7 8 9 10; 1 2 3 4 5;];

 for i=1:n

     

     a(i,j:n) = a(i,j:n)./n;

    

 end

 disp(a)

0 'for loop':

i=1;

j=1;

n=5;

a=[1 2 3 4 5; 6 7 8 9 10; 1 2 3 4 5; 6 7 8 9 10; 1 2 3 4 5;];

 

a(i:n,j:n) = a(i:n,j:n)./n;

disp(a)

In all the above cases, the results are the same. You may run the codes and check for yourself. The operation is pretty straightforward here. When colon (:) is first used between j and n,  we eliminate a 'for loop'. This colon (:) is considering all the elements in columns of the matrix a, which is equivalent to the 'loop' operation. The second colon (:), between i and n (i:n), takes all the row elements of the matrix a, and the operation is also same compared to using a 'for loop'. Below is the result for all three demonstrated cases.

    0.2000    0.4000    0.6000    0.8000    1.0000

    1.2000    1.4000    1.6000    1.8000    2.0000

    0.2000    0.4000    0.6000    0.8000    1.0000

    1.2000    1.4000    1.6000    1.8000    2.0000

    0.2000    0.4000    0.6000    0.8000    1.0000

This is a very simple example, and for this, you won't even realize the actual computational efficiency. But, imagine! you have very large matrices where numerous complex operations involved. In that case, you would be able to see the differences. 

#Vectorization #Matrix #ForLoop #UseofColon #Matlab #Blog #Blogger

Comparison between 4 and 8 Node Integration Elements in Finite Element Analysis

Let us assume a cantilever beam with one end fixed and the other end has 1 unit load acting vertically downward. The length and height of the beam are 5 and 1 unit respectively. The Young's modulus is 10,000 unit. At first, we will solve this problem analytically using the famous Euler-Bernoulli beam equation to find the deflection at the end of the beam. Next, we will use Abaqus with 4 node quadrilateral  reduced integration element. Then, we will solve the same problem using 8 node quadrilateral reduced integration element. Our goal is to see which element performs better in terms of convergence and accuracy for the finite element analysis (FEA) compared to the analytical solution.

The following Mathematica program solves the Euler-Bernoulli equation for the deflection at the tip.

(*Exact Solution of Euler-Bernoulli Cantilever Beam*)

Clear[y, P, x, X2, EI, L];

y[x];

(*Given Property Values*)

L = 5;

h = 1;

P = 1;

Ie = 1/12;

Ee = 10000;

EI = Ee*Ie;

theta = D[y[x], x];

V = EI*D[y[x], {x, 3}];

M = EI*D[y[x], {x, 2}];

y1 = y[x] /. x -> 0; y2 = y[x] /. x -> L; theta1 = 

 theta /. x -> 0; theta2 = theta /. x -> L; M1 = M /. x -> 0; M2 = 

 M /. x -> L;

V1 = V /. x -> 0; V2 = V /. x -> L;

(*Solving the Equation*)

s = DSolve[{EI*y''''[x] == 0, y1 == 0, theta1 == 0, V2 == P, M2 == 0},

    y, x];

Print["Analytical Solution of the Exact Displacement"]

y = Simplify[y[x] /. s[[1]]]

Plot[y, {x, 0, 5}, PlotStyle -> {Blue, Thick}]

Print["Displacement Variation at Point A in both X and Y Direction"]

DA = y /. x -> 5;

theta = D[y, x] /. x -> L;

a1 = N[L - (-h/2)*Sin[theta]];

a2 = N[DA + (-h/2)*Cos[theta]];

A1 = a1 - 5

A2 = a2 + h/2

Print["The Stress at Point B, S11"]

S11 = -Ee*X2*D[y, {x, 2}];

StressatB = S11 /. {x -> 0, X2 -> 0.5}

Now, we begin to simulate the same problem with two types of elements: 4 and 8 node quadrilateral reduced integration. As in any FEA simulation, we are required to define the meshing parameters and elements prior to the analysis, and these two elements are frequently used in Abaqus. 

The above results highlight that there is a trend of increasing convergence of the results starting from the fourth elements. We see that the stresses have minute changes in its absolute values, which are also visible for the displacements. With the increase of the element numbers, the faster the different plots to converge. It is evident that the elements are not converging to the same solutions always. Moreover, for this problem, the rate of convergence and accuracy is better with the 8 node quadrilateral reduced integration element; however, with the increase in number of elements, both 4 and 8 node quadrilateral reduced integration elements have similar performance.

If you would like to gather more knowledge on FEA, I would encourage you to visit this site as this tutorial is based on the discussions and examples of that.

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

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.

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