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