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

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:

Stiffness Matrix

Force Vector
Displacement Vector
Stress and Strain matrix
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.

Variations of Displacements

Variations of Displacements










The following figure shows the variations of stress in magnitude.


Variations of Stress

The Output File for the Stresses:


The Output File for the Stresses

The Output File for the Strains:


The Output File for the Strains

Comparison between Stiffness Matrices:

Mathematica Calculation


Stiffness Matrix

From Abaqus


Stiffness Matrix
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.

A Short Introduction to MATLAB for Beginners

In this post, you will find series of basic tutorial videos on MATLAB. As you may know that MATLAB is a very powerful computational tool, which has been widely used in academics as well as industries. You may find these videos useful if you haven't yet used this tool. 


   

  
  


2D Laplace Equation Solution by 5 Point Finite Difference Approximation

The temperature distribution of a rectangular plate is described by the following two dimensional (2D) Laplace equation:

Txx + Tyy = 0

The width (w), height (h), and thickness (t) of the plate are 10, 15, 1 cm, respectively. There is a heat source at the top edge, which is described as, T = 100 sin (πx / w) Celsius, and all other three edges are kept at 00 C. There is no heat flow along the plate thickness since the faces are insulated. Also, assume that there is no internal heat generation so that the problem becomes 2D temperature distribution.

Solve the Laplace equation with a 5 point finite difference approximation by MATLAB where Δx = Δy = 2.5 cm. The temperature distribution is defined by T(x,y)

The following MATLAB program develops a 5-point approximation stencil to solve the 2D Laplace equation.

MATLAB Program:

close all;
clc;
% Width, height & thickness of the thin plate (dimension, cm)
W=10.0;
h=15.0;
t=1.0;
% Gird size definition along x and y axes
xx=2.5;
yy=2.5;

imax=(W/xx)-1; % Array dimension along the width
jmax=(h/yy)-1; % Array dimension along the height
N_array=imax*jmax; % Total number of grid points

% Generating the Matrix, A
for j=1:N_array
    for i=1:N_array
        if i==j
            a(i,j)=-4.0;
        end 
        if j==(i+1)
            a(i,j)=1.0;
        end
        if j==(i-1)
            a(i,j)=1.0;
        end
        if j==(i+5)
           a(i,j)=1.0;
        end
        if j==(i-5)
           a(i,j)=1.0;
        end
        if i==(imax+2)
            a(i,i+1)=0.0;
        end
        if i==(imax+3)
            a(i,i-1)=0.0;
        end
        if i==(imax+2+imax+2)
            a(i,i+1)=0.0;
        end
        if i==(imax+2+imax+2+1)
            a(i,i-1)=0.0;
        end
    end
end

for i=1:imax
    x(i)=(i)*xx;
    Tb(i)=100*sin(pi*x(i)/W); % Top edge temperature (in Celcius)
end

% Generating the Matrix, b
k=imax+2;
for i=1:imax
    b(k)=-Tb(i);
    k=k+imax+2;
end

%Displaying A and b matrices
B=b';
T=GaussElimination(a,B); % Calling the function defined elsewhere
TT(:,1) = T(1:5);
TT(:,2) = T(6:10);
TT(:,3) = T(11:15);
TT = [TT;Tb];
TT = padarray(TT,[1 0],'pre');
TT = padarray(TT, [0 1]);

imagesc([0:xx:xx*imax],[0:yy:yy*jmax],TT)
h=colorbar;
colormap jet;
ylabel(h,'Temperature, C')
set(gca,'YDir','normal')
xlabel('Width, cm')
ylabel('Height, cm')


The following function is created in MATLAB to solve the system of algebraic equations using Gauss elimination method:


function x = GaussElimination(A,b)
% Gauss Elimination without Pivoting
% x = GaussElimination(A,b)
% input:
% A = Coefficient Matrix
% b = Right Hand Side Vector
% output:
% x = Solution Vector

[m,n] = size(A);
if m~=n
    error('Matrix A must be square');
end
nb = n+1;
ATM = [A b];

% Forward Elimination
for k = 1:n-1
for i = k+1:n
factor = ATM(i,k)/ATM(k,k);
ATM(i,k:nb) = ATM(i,k:nb)-factor*ATM(k,k:nb);
end
end

% Backward Substitution
x = zeros(n,1);
x(n) = ATM(n,nb)/ATM(n,n);
for i = n-1:-1:1
x(i) = (ATM(i,nb)-ATM(i,i+1:n)*x(i+1:n))/ATM(i,i);
end

disp('Solution Vector')
disp(x)
disp('Coefficient Matrix')
a1 = ATM(:,1:n); 
disp(a1)
disp('The b Vector')
b1 = ATM(:,nb);
disp(b1)
disp('Verification of the Solution')
C = a1*x-b1;  
disp(C)


Program Output:

Solution Vector
    1.3046
    3.3734
    7.4183
   15.8087
   33.4596
    1.8450
    4.7707
   10.4910
   22.3568
   47.3190
    1.3046
    3.3734
    7.4183
   15.8087
   33.4596

Coefficient Matrix
Columns 1 through 11

-4.0000    1.0000         0         0         0    1.0000         0         0         0         0         0
      0   -3.7500    1.0000         0         0    0.2500    1.0000         0         0         0         0
      0         0   -3.7333    1.0000         0    0.0667    0.2667    1.0000         0         0         0
      0         0         0   -3.7321    1.0000    0.0179    0.0714    0.2679    1.0000         0         0
      0         0         0         0   -3.7321    0.0048    0.0191    0.0718    0.2679    1.0000         0
      0         0         0         0         0   -3.7321    1.0718    0.0192    0.0051    0.0013    1.0000
      0         0         0         0         0         0   -3.4050    1.0824    0.0220    0.0055    0.2872
      0         0         0         0         0         0         0   -3.3673    1.0839    0.0210    0.0964
      0         0         0         0         0         0         0         0   -3.3638    1.0786    0.0343
      0         0         0         0         0         0         0         0         0   -3.3861    0.0124
      0         0         0         0         0         0         0         0         0         0   -3.7047
      0         0         0         0         0         0         0         0         0         0         0
      0         0         0         0         0         0         0         0         0         0         0
      0         0         0         0         0         0         0         0         0         0         0
      0         0         0         0         0         0         0         0         0         0         0

  Columns 12 through 15

         0         0         0         0
         0         0         0         0
         0         0         0         0
         0         0         0         0
         0         0         0         0
         0         0         0         0
    1.0000         0         0         0
    0.3179    1.0000         0         0
    0.1088    0.3219    1.0000         0
    0.0385    0.1094    0.3206    1.0000
    1.0947    0.0323    0.0114    0.0037
   -3.3489    1.1156    0.0393    0.0124
         0   -3.2968    1.1193    0.0365
         0         0   -3.2919    1.1072
         0         0         0   -3.3318

The b Vector
         0
         0
         0
         0
  -70.7107
   -0.0907
   -0.3887
   -1.4838
   -5.5569
 -120.7386
   -0.5983
   -1.9828
   -5.5408
  -14.9919
 -111.4802

Verification of the Solution
   1.0e-13 *

   -0.0022
         0
         0
   -0.0711
    0.1421
    0.0053
   -0.0078
   -0.0644
   -0.1421
         0
   -0.0011
   -0.0044
    0.0089
    0.0355
         0

The following plot shows the desired temperature distribution within the plate.
Temperature distribution in the plate