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.