^{3}. 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:__

**In Abaqus, the geometry is defined by using the given nodal points in assignment. The material is elastic, where Young’s modulus is 20000000000**

__Using Abaqus:__

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

## No comments:

## Post a Comment