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 Output File for the Stresses:
The Output File for the Strains:
Comparison between Stiffness Matrices:
Mathematica Calculation
From Abaqus