2. For the problem shown in Fig.

3, use any type of elements and determine for the reaction

at the support, the displacements at each nodes, the principal stresses, and principal
strains distribution on a vertical plane 0.125 m from the fixed support. Data:  = 𝑦𝑜𝑢𝑟
𝑎𝑔𝑒 𝑖𝑛 𝑦𝑟𝑠/ 150.000, = 𝑦𝑜𝑢𝑟 𝑎𝑔𝑒 𝑖𝑛 𝑦𝑟𝑠/ 0.1500 𝑘𝑁/𝑚2, t=0.030 m, E=209 Gpa, P =
𝑦𝑜𝑢𝑟 𝑎𝑔𝑒 𝑖𝑛 𝑦𝑟𝑠/ 2.500 𝑘𝑁 and your calculation should be in three decimal places.

Given data

 w= 27/0.2 = - 135.000 kN /m2 ,

 P=27/2.5= 10.8 kN ,
 E = 208 Gpa ,
 v = 27/100 = 0.270,
 t=0.02 m ,
Required

 Displacement at each node

 Reactions forces at support
 Major principal stress on vertical plane 0.125 from support
Solution

Solution for this problem is done using FEM by following necessary steps as plain stress
problem analysis.

ASTU Page 1
Step 1- Discretizing the Domain:

Consider the problem as plane stress problem. And it is used linear triangular element to
discretize and a having a total of six elements with eight nodes.

Figure discretization of elements using linear triangle

The total force due to distributed load is divided in to node 5 and 8. With resultant force of
0.3375 kN downward and 10.4625 kN upward respectively.

Element Connectivity for the system

Element number Node (i) Node (j) Node (m)

1 1 5 6

2 1 2 5

3 2 4 5

4 2 3 4

5 5 8 7

6 5 7 6

ASTU Page 2
Step 2 -Writing the Element Stiffness Matrices

Six element stiffness matrices are obtained by making call to MATLAB Linear triangular
element stiffness. Each matrix has size of 6x6.

>> E=208e9;

>> t=0.02;

>> age=27;

>> NU=age/100;

>> % Element Stiffness Matrix

>> k1 = LinearTriangleElementStiffness1(E,NU,t,0,0,0.25,0.25,0,0.25,1);

>> k2 = LinearTriangleElementStiffness1(E,NU,t,0,0,0.25,0,0.25,0.25,1);

>> k3 = LinearTriangleElementStiffness1(E,NU,t,0.25,0,0.5,0.25,0.25,0.25,1);

>> k4 = LinearTriangleElementStiffness1(E,NU,t,0.25,0,0.5,0,0.5,0.25,1);

>> k5 = LinearTriangleElementStiffness1(E,NU,t,0.25,0.25,0.25,0.5,0,0.5,1);

>> k6 = LinearTriangleElementStiffness1(E,NU,t,0.25,0.25,0,0.5,0,0.25,1);

Step 3 –Assembling the Global Stiffness Matrix:

The system has 8 nodes thus, global stiffness matrix K become 16x16. So to found K first set up
zero matrix of size 16x16 and call for MATLAB function of K=LinearTriangleAssemble . As the
matrix is large it is suppressed in the analysis.

>> % assemble of ESM to get global stiffness matrix

>> K=zeros(16);

>> K = LinearTriangleAssemble(K,k1,1,5,6);

>> K = LinearTriangleAssemble(K,k2,1,2,5);

>> K = LinearTriangleAssemble(K,k3,2,4,5);

ASTU Page 3
>> K = LinearTriangleAssemble(K,k4,2,3,4);

>> K = LinearTriangleAssemble(K,k5,5,8,7);

>> K = LinearTriangleAssemble(K,k6,5,7,6);

Step 4 – Applying the Boundary Conditions:

The boundary conditions are:

 U1x= U1y=U6x=U6y=U7x=U7y= 0
 F2x=F2y=F3x=F3y=F4x=F5x=F8x=F8y=0
 F5y= -0.3375*103 N, F4y=10.4625*103 N
Step 5 – Solving the Equation:

Solve the system of equation in step 3 will performed by partitioning manually and use
MATLAB. First call submatrix

>> % applying boundary conditions

>> Kp=[K(3:10,3:10),K(3:10,15:16);K(15:16,3:10),K(15:16,15:16)]; % partioned stiffness

matrix

>> Fp=[0,0,0,0,0,10462.5,0,-337.5,0,0]';

>> Up=Kp\Fp;

>> % global nodal displacement vector

>> U=[0;0;Up(1:8);0;0;0;0;Up(9:10)];

Nodal displacement vector

ASTU Page 4
Step 6–post – processing:

Now it is possible to obtain nodal forces

Nodal Force Vector

ASTU Page 5
As ∑ Fx= 0 and ∑ Fy= 0 force equilibrium is satisfied.

Then set up the element nodal displacement vector and solve for element stress sigma.

>> U=[0;0;Up(1:8);0;0;0;0;Up(9:10)];

>> u1=[U(1:2);U(9:12)];

>> u2=[U(1:4);U(9:10)];

>> u3=[U(3:4);U(7:10)];

>> u4=[U(3:8)];

>> u5=[U(9:10);U(13:16)];

>> u6=[U(9:10);U(11:14)];

>> str1 = LinearTriangleElementStresses1(E,NU,t,0,0,0.25,0.25,0,0.25,1,u1);

>> str2 = LinearTriangleElementStresses1(E,NU,t,0,0,0.25,0,0.25,0.25,1,u2);

>> str3 = LinearTriangleElementStresses1(E,NU,t,0.25,0,0.5,0.25,0.25,0.25,1,u3);

>> str4 = LinearTriangleElementStresses1(E,NU,t,0.25,0,0.5,0,0.5,0.25,1,u4);

>> str5 = LinearTriangleElementStresses1(E,NU,t,0.25,0.25,0.25,0.5,0,0.5,1,u5);

>> str6 = LinearTriangleElementStresses1(E,NU,t,0.25,0.25,0,0.5,0,0.25,1,u5);

Now it is possible to calculate the principal stresses and principal angle for each element by
making

call to the MATLAB function .

>> pstr1 = LinearTriangleElementPStresses(str1)

pstr1 =

1.0e+06 *

0.7357

-2.4370

-0.0000

ASTU Page 6
pstr2 = LinearTriangleElementPStresses(str2)

pstr2 =

1.0e+06 *

3.9639

0.4101

-0.0000

>> pstr3 = LinearTriangleElementPStresses(str3);

>> pstr4 = LinearTriangleElementPStresses(str4);

>> pstr5 = LinearTriangleElementPStresses(str5);

>> pstr5 = LinearTriangleElementPStresses(str5)

>> pstr6 = LinearTriangleElementPStresses(str6)

The major principal stress distribution on a vertical face 0.125m from the fixed support is similar
to the principal stresses calculated at the cancroids for element no 1 ,2 , 5 and 6 as the problem is
solved by CST method

Element Coordinates Major principal stress

No.
x(m) Y(m) (MPa)

1 0.125 0.167 2.43

2 0.125 0.083 0.41

5 0.125 0.42 3.97

6 0.125 0.33 2.43

ASTU Page 7
 2. The rectangular beam in Fig. 3 has a square opening of side length S. Use the load
magnitudes, type and number of elements, length of side of the square given in Table- 1
to compute the displacement, forces at A,B, C and at the load application points. Point B
is at center and use the material data given in Question number 1.

Given data

The material is considered to have similar properties as given on question 1

E = 209 Gpa , v = 27/100 = 0.27 , t=0.02m, S= 0.25m

W= 125kN /m, F1=18 kN , F2=20 kN

- Analysis Type – Plane stress.

- Finite element type – Constant Strain Triangle (CST)

Required

- Displacement at A, B, C and at load application points

- Forces at A ,B ,C and at load application points

Solution

ASTU Page 8
Solution for this problem is done using FEM by following necessary steps as plane stress
problem analysis.

Step 1- Discretizing the Domain:

The problem is needed to analyze as plane stress problem. And it is used Constant Strain
Triangle (CST) to discretize and a having a total of 12 elements with 12 nodes.

Element Connectivity for the system

Element Nodes

Node i Node j Node m

1 1 4 12

2 1 2 4

3 2 5 4

ASTU Page 9
 4 2 3 5

5 3 6 5

6 5 6 7

7 6 8 7

8 7 8 9

9 7 9 10

10 10 9 11

11 10 11 12

12 12 4 10

Step 2 -Writing the Element Stiffness Matrices

A total of 12 element stiffness matrices are obtained by making call to MATLAB

LinearTriangleElementStiffness1. Each matrix has size of 6x6.

>> E=209e9;

>> NU=0.27;

>> t=0.02;

>> k1 = LinearTriangleElementStiffness1(E,NU,t,0,0,0.375,0.25,0,0.375,1);

>> k2 = LinearTriangleElementStiffness1(E,NU,t,0,0,0.5,0,0.375,0.25,1);

>> k3 = LinearTriangleElementStiffness1(E,NU,t,0.5,0,0.625,0.25,0.375,0.25,1);

>> k4 = LinearTriangleElementStiffness1(E,NU,t,0.5,0,1,0,0.625,0.25,1);

>> k5 = LinearTriangleElementStiffness1(E,NU,t,1,0,1,0.375,0.625,0.25,1);

>> k6 = LinearTriangleElementStiffness1(E,NU,t,0.625,0.25,1,0.375,0.625,0.5,1);

>> k7 = LinearTriangleElementStiffness1(E,NU,t,1,0.375,1,0.75,0.625,0.5,1);

>> k8 = LinearTriangleElementStiffness1(E,NU,t,0.625,0.5,1,0.75,0.5,0.75,1);

>> k9 = LinearTriangleElementStiffness1(E,NU,t,0.625,0.5,0.5,0.75,0.375,0.5,1);

ASTU Page 10
 >> k10 = LinearTriangleElementStiffness1(E,NU,t,0.375,0.5,0.5,0.75,0,0.75,1);

>> k11 = LinearTriangleElementStiffness1(E,NU,t,0.375,0.5,0,0.75,0,0.375,1);

>> k12 = LinearTriangleElementStiffness1(E,NU,t,0,0.375,0.375,0.25,0.375,0.5,1);

Step 3 –Assembling the Global Stiffness Matrix:

The system has 12 nodes thus, global stiffness matrix K become 24x24. So to found K first
set up zero matrix of size K and call for MATLAB function of K=LinearTriangleAssemble .

>> % Assembling the ESM to get Global stiffness matrix

>> K=zeros(24);

>> K = LinearTriangleAssemble(K,k1,1,4,12);

>> K = LinearTriangleAssemble(K,k2,1,2,4);

>> K = LinearTriangleAssemble(K,k3,2,5,4);

>> K = LinearTriangleAssemble(K,k4,2,3,5);

>> K = LinearTriangleAssemble(K,k5,3,6,5);

>> K = LinearTriangleAssemble(K,k6,5,6,7);

>> K = LinearTriangleAssemble(K,k7,6,8,7);

>> K = LinearTriangleAssemble(K,k8,7,8,9);

>> K = LinearTriangleAssemble(K,k9,7,9,10);

>> K = LinearTriangleAssemble(K,k10,10,9,11);

>> K = LinearTriangleAssemble(K,k11,10,11,12);

>> K = LinearTriangleAssemble(K,k12,12,4,10);

Step 3 – Applying the Boundary Conditions:

The boundary conditions are:

Su=Essential Boundary condition

U1=U2=U21=U22=U23=U24=0

ASTU Page 11
 St= Natural Boundary Condition

F3=F4=F5=F6=F7=F8=F9=F10=F11=F12=F13=F14=F17=F19=F20=0

F5=20, F15=18, F16=-31.25, F18=-62.5, F22=-31.25

Step 4 – Solving the Equation:

Solve the system of equation in step 3 will performed by partitioning manually and use
MATLAB. First call Partitioned Stiffness matrix Kp

>> % partitioned stiffness matrix

>> Kp=K(3:20,3:20);

>> Fp=[0,0,20,0,0,0,0,0,0,0,0,0,18,-31.25,0,-62.5,0,0]';

>> up=inv(Kp)*Fp; % Unknown nodal displacements vector

Step 5 –post – processing:

Now it is possible to obtain nodal forces

>> U=[0,0,up',0,0,0,0]'; %Global displacement vector

>> F=K*U; % Global force vector

∑ F x =¿ ∑ F y =0 ¿ force equilibrium is satisfied.

Nodal forces and nodal displacements at the required points

Point(Node Coordinates Nodal force (kN) nodal displacement(m)

no) 1*e-6

X Y Fx Fy Ux Uy

A(1) 0 0 69.08 15.96 0 0

B(12) 0 0.375 -11.49 59.81 0 0

C(11) 0 0.75 -95.58 17.69 0 0

3 1 0 20 0 -0.0236 -0.1403

8 1 0.75 18 -31.25 0.0653 -0.1564

ASTU Page 12
 5. Consider the footing in Fig. 5.

a) Assume that the two layers in the figure are of the same type of soil (Sand Soil). Compare
the stress at a depth Z from the surface using both numerical method and Boussinesq
formula for point load.

Given data:

Boundary limits X1=4.5 m, X2= 4.5m and assume Z2=2m, B=2m

Distributed Load = 120KN/m

Solution

Boussinesq (1883) solved the problem for stresses inside a semi-infinite mass due to a point load
acting on the surface. In rectangular coordinates, the stresses may be expressed as follows:

ASTU Page 13
At center X=0, Y=0

Q (point laod) =q*B=120*2=240KN

Therefore, result from numerical method and Boussinesq formula for point load will be:

Depth (m) Boussineq’s stress (Kpa)

1 114.59

2 28.65

3 12.73

4 7.16

5 4.58

6 3.18

7 2.33

8 1.7

ASTU Page 14
Solution using FEM

The solution for question 5 and 6 is done by Linear Constant Strain (LST) method. Or using
MATLAB function QuadTriangleElementStiffness, but the final result of the analysis cannot
defined by the program, the program displays the solution as not a number(Nan).

b) Assume that you have two layers as given in the figure above. Use both theoretical estimation
and computational method (FEM) to estimate the stresses, settlement due to the distributed load
from the water tank and self-weight of the soil at depth Z below the surface. Use the following
data.

Solution

For the first layer the stress distribution below the center of the loaded area is calculated using
the ‘Boussinesq’ formula for a uniform vertical loading on a circular area:

Where, b = 2m and q = 120 KN/m2 γ=19 KN/m2

Depth-z (m) b/z Ic ∆𝛔z=q*Ic γh Total

0.5 2 0.911 109.32 9.5 118.82

1 1 0.646 77.52 19 96.52

1.5 0.667 0.424 50.88 28.5 79.38

2 0.5 0.284 34.08 38 72.05

2.5 0.4 0.2 24 47.5 71.5

3 0.33 0.146 17.52 57 74.52

3.5 0.286 0.11 13.2 66.5 79.7

ASTU Page 15
 4 0.25 0.087 10.44 76 86.44

For two-layer system using Poisson’s ratio method, Equivalent Poisson’s ratio is used in

Westergaard (1938) equation to evaluate the stresses as given below

H=6m, h1=4m , h2=2m, μ1=0.2 μ2= 0.25, and μeq = 0.231

Equivalent Poisson’s ratio, η = 0.35

Depth-z(m) b/z Ic ∆𝛔z=q*Ic γh Total

4.5 0.22 0.412 49.44 85.5 89.9

5 0.2 0.399 47.88 95 142.9

5.5 0.182 0.343 41.16 104 145.2

6 0.167 0.322 38.64 114 152.6

ASTU Page 16