National Polytechnique School -Mechanical department

Academic Year 2024-2025

Option: Mechanical Engineering
Semester: 3
Tutorial
Fundamental Teaching Unit: UEF 211
Subject: Finite Element Method CST Element
EXAMPLE Constant Strain Triangle
For a thin plate subjected to the surface traction shown in Figure 01, determine the nodal
displacements and the element stresses. The plate thickness t = 20 mm, E = 210 GPa,
and
0.30.
V =

400mm
T = 7MPa

■ Figure 01 Thin plate subjected to tensile stress

SOLUTION:
Discretization
To illustrate the finite element method solution for the plate, we first discretize the plate
into two elements, as shown in Figure 02. It should be understood that the coarseness of
the mesh will not yield as true a predicted behavior of the plate as would a finer mesh,
particularly near the fixed edge. However, since we are performing a longhand solution, we
will use a coarse discretization for simplicity (but without loss of generality of the
method).
In Figure 02, the original tensile surface traction in Figure 01 bas been converted to
nodal forces as follows:
1
F = -TA
2
F = _!_ (7 X 106 )(20 X 10-3)(200 X 10-3)
2
F = 14,000 N

■ Figure 02 Discr etized plate

01
 The governing global matrix equation is

{F} = [K]{d} (1)

Expanding matrices in Eq. (1), we obtain

Fix Rix u1 0
Fix R1y V1 0
F2, R2x U2 0
F2y v2 0
R2y
= [K] = [K] (2)
14,000 U3 U3
F3x

0
F3y V3 V3

F4x 14,000 U4 U4

F4 y 0 V4 V4

where [K] is an 8 X 8 matrix (two degrees of freedom per node with four nodes) before delet­
ing rows and columns to account for the :fixed boundary support conditions at nodes 1 and 2.

Assemblage of the Stiffness Matrix

We assemble the global stiffness matrix by superposition of the individual element
stiffness matrices. The stiffness matrix for an element is

[ k] = tA[Bf [D][B] (3)

In Figure 3 for element 1, we have coordinates Xi = 0, Yi = 0, xj = 400 mm, yj = 10, Xm
= 0, and Ym = 200 mm, since the global coordinate axes are set up at node 1, and

A= .!_ bh

A= ( ½ }o.4)(0.2) = 0.04 m 2

m=2 j=3

i = I

■ Figure 3 Element 1 of the discretized plate

We will now evaluate [B], where [B] is given by

:. l
0 /3j 0 f3m

[B] = _l 0 0 0 (4)

[�
'Yi ,'j
2A
')'; /3i ,'j /3j 'Ym f3m

02
and,

/3; = Yj - Ym = 200 - 200 = 0

/3j = Ym - y; = 200 - 0 = 200
f3m = Y; - Yj = 0 - 200 = -200
(5)
'Yi =xm -xj = 0 - 400 = -400
o/j = X; =0-0=0
- Xm

'Ym =xj -x; = 400 - 0 = 400

Therefore, substituting Eqs. (5) into Eq. (4), we obtain

0 0 10 0 -10
[ ] 20 X lQ-3[ 1
B = 0 -20 0 0 0 2�]- (6)
2 X 4 X 10-2 -20 0 0 10 20 -10
m

[ ]
For plane stress, the D matrix is conveniently expressed here as

1 V
E [v 1
[D] = (7)
(1 - y 2 )
0 0

With v = 0.3 and E = 210 X 109 N/m, we obtain

[ D]
=
210(l09 )
[10.3
0.3 1 � ]N/m (8)
0.91
O O 0.35

0 0 -20

L]
0 -20 0
0.3
(20 X 10-3)(210 X 10 9 ) 10 0 0[�3 1
Then [B]T[D] = (9)
(2 X 4 X 10-2)(0.91) 0 0 10
0
-10 0 20
0 20 -10

Simplifying Eq. (9) yields

0 0 -7
-6 -20 0
(52.5)(109 ) 10 3 0
[Bf[D] = (10)
0.91 0 0 3.5
-10 -3 7
6 20 -3.5

03
 Using Eqs. (10) and (6) in Eq. (3), we have the stiffness matrix for element 1 as

0 0 -7
-6 -20 0
(52.5)(10 9)
3 10 3 0
[k<ll] = (0.04)(20 X 10- ) 3.5
0.91 0 0
-10 -3 7
6 20 -3.5

20 -10ol
0 10 0 -10
20 X 10-s [ O
X O -20 0 0 0 20 (11)
2 X 4 X 10- 2 -20 0 0 10

Finally, simplifying Eq. (11) y ields

UJ VJ U3 V3 Uz Vz

140 0 0 -70 -140 70

0 400 -60 0 60 -400 N
(12)
[k(ll] = 10.5 X 10
6
0 -60 100 0 -100 60
m
0.91 -70 0 0 35 70 -35
-140 60 -100 70 240 -130
70 -400 60 -35 -130 435
where the labels above the columns indicate the counterclockwise nodal order of the degrees
of freedom in the element 1 stiffness matrix.

m = 3

i = 1 j = 4

■ Figure 4 Element 2 of the discretized plate

In Figure 4 for element 2, we ha ve x; = 0, y; = 0, Xj = 400 mm, yj = 0,

Xm = 400 mm, and Ym = 200 mm. Then, we have

/3; = Yj - Ym = 0 - 200 = -200

/3j = Ym - y; = 200 - 0 = 200
f3m = y; - yj = 0 - 0 = 0
(13)
'Yi = Xm - Xj = 400 - 400 = 0
'Yj = X; - Xm = 0 - 400 = -400
'Ym =xj -x; = 400 - 0 = 400

04
Therefore, using Eqs. (6.5.13) in Eq. (6.5.4) yields

-IO 0 10 0 0
20x10-'[ 0 0 -20 0
-1 (14)
[B] = rn
(2)(0.04) � -10 -20 10 20 0i

�
The [D] rnatrix is again given by

[D] =
0.91
' [1
0
0.3
210(10 ) 0.3 1
0
� ]Nhn
0.35
(15)

Then, using Eqs. (6.5.14) and (6.5.15), we obtain

�l
-10 0 0
0 0 -10
0.3
(2 X 10-3)(210 X 109 ) 10 0 -20 1
[�3
(16)
[Bf[D] = 0 -20 10
(2)(0.04)(0.91) 0 0.35
0 0 20
0 20 0

Sirnplifying Eq. (16) yields

-10 -3 0
0 0 -3.5
(52.5)(109 ) 10 3 -7 (17)
[Bf[D] = -6 -20 3.5
0.91
0 0 7
6 20 0

Finally, substituting Eqs. (17) and (14) into Eq. (3), we obtain the stiffness rnatrix for elernent
2 as
-10 -3 0
0 0 -3.5
(52.5)(109 ) 10 3 -7
10 -3 )(0.04)
[k<2l] = (20 X
0.91 -6 -20 3.5
0 0 7
6 20 0

�]
-!0 0 10 0 0
20 X 10-'[ 0 0 0 -20 0 (18)
X
2(0.04) 0 -10 -20 10 20

05
Equation (18) simplifies to

U1 V1 U4 V4 U3 V3

100 0 -100 60 0 -60

0 35 70 -35 -70 0
= 10.5 X 10 -100 70 240 -130 -140 (19)
6
60 N
[k<2l]
0.91 60 -35 -130 435 70 -400 m
0 -70 -140 70 140 0
-60 0 60 -400 0 400

where the degrees of freedom in the element 2 stiffness matrix are shown above the col-
umns in Eq. (19). Rewriting the element stiffness matrices, Eqs. (12) and (19), expanded
to the order of, and rearranged according to, increasing nodal degrees of freedom of the
total [K] matrix (where we have factored out a constant 5), we obtain

Element 1

U1 V1 U2 V2 U3 V3 U4 V4

28 0 -28 14 0 -14 0 0
0 80 12 -80 -12 0 0 0
-28 12 48 -26 -20 14 0 0
52.5 X 10 6 14 -80 -26 87 12 -7 0 0 N (20)
[k(ll] =
0.91 0 -12 -20 12 20 0 0 0 m
-14 0 14 -7 0 7 0 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0

Element 2

U1 V1 U2 V2 U3 V3 U4 V4

20 0 0 0 0 -12 -20 12
0 7 0 0 -14 0 14 -7
0 0 0 0 0 0 0 0
52.5 X 10 6 0 0 0 0 0 0 0 0 N (21)
[k<2l] =
0.91 0 -14 0 0 28 0 -28 14 m
-12 0 0 0 0 80 12 -80
-20 14 0 0 -28 12 48 -26
12 -7 0 0 14 -80 -26 87

Using superposition of the element stiffness matrices, Eqs. (20) and (21), now that the
orders of the degrees of freedom are the same, we obtain the total global stiffness
matrix as

06
 U1 V1 U2 V2 U3 V3 U4 V4

480 -28 14 0 -26 -20 12

87
0 12 -80 -26 0 14 -7
-28 12 48 -26 -20 14 0 0
52.5 X 106 14 -80 -26 87 12 -7 0 0 N (22)
[K] =
0 .91 0 -26 -20 12 48 0 -28 14 m
-26 0 14 -7 0 87 12 -80
-20 14 0 0 -28 12 48 -26
12 -7 0 0 14 -80 -26 87

[Altematively, we could have applied the direct stiffness method to Eqs. (12) and
(19) to obtain Eq. (22).] Substituting [K] into {F} = [K]{d} ofEq. (2), we have

Rix
48 -28
0 14 0 -26 -20 12 0
R 1y
0 87 12 -80 -26 0 14 -7 0

R2x -28 12 48 -26 -20 14 0 0 0

R2y 52.5 X 10 9 14 -80 -26 87 12 -7 0 0 0

U3 (23)
14, 000 0 .91 0 -26 -20 12 48 0 -28 14
0
-26 0 14 -7 0 87 12 -80 V3

14, 000 -20 14 0 0 -28 12 48 -26 U4

12 -7 0 0 14 -80 -26 87 V4
0

Applying the support or boundary conditions by eliminating rows and columns correspond­
ing to displacement matrix rows and columns equal to zero [ namely, rows and columns
1-4 in Eq. (23)], we obtain

14, 000 48 0 -28

6
0 52.5 X 10 0 87
1 } - [ (24)
14, 000 0 .91 -28 12
0 14 -80 -26

Premultiplying both sides ofEq. (6.5.24) by [Kr 1, we have

48
V3 0 .91 0 87 [
(25)
( U4 1 - 52.5 X 10 6 -28 12
U3
V4 14 -80

Solving for the displacements in Eq. (25), we obtain

l
0 .05024
V3 0 .91 0 .000 34
= 1 } (26)
U4 3750 0 .05470
U3 }
V4 0 .00 878

07
 Simplifying Eq. (26), the final displacements are given by

U3 12.19
V3 0.083
= X 10_ 6 m (27)
(U4 1 {13.27}
V4 2.08

Comparing the finite element solution to an analytical solution, as a first approximation,

we have the axial displacement given by
PL (28, 000)0.4
S = = = 13_33 X 10 _6 m
AE (20 X 10-3)(200 X 10-3)(210 X 109 )

for a one-dimensional bar subjected to tensile force. Hence, the nodal x displacement com­
ponents ofEq. (6.5.27) for the two-dimensional plate appear to be reasonably correct, con­
sidering the coarseness of the mesh and the directional stiffness bias of the model. (For more
on this subject see Section 7.5.) The y displacement would be expected to be downward at
the top (node 3) and upward at the bottom (node 4) as a result of the Poisson effect. However,
the directional stiffness bias due to the coarse mesh accounts for this unexpected poor result.
We now determine the stresses in each element by using the equation :

{u} = [D][B]{d} (28)

In general, for element 1, we then have

U]

1 (2Jr�
�1
V V]
0 /33 0 /32
E V 1 U3
{u} = i 'YI 0 'Y3 0 (29)
v 2
(1 - ) r 1� y X V3
0 0 'YI /31 'Y3 /33 'Y2 /32
2 U2
V2

Substituting numerical values for [B], given by Eq. (6); for [D], given by Eq. (8); and
the appropriate part of {d}, given by Eq. (27), we obtain

L]
[I 0.3
210(109 )(10-')
{u} = 0.3 1
0.91(4)
0 0

[ 0 20o
0 (30)

l
0 10 0 -10
12.19
X 0 -20 0 0 0
0.083
-20 0 0 10 20 -10
0
0
Simplifying Eq. (30), we obtain

08
 t'}
ay = 2110 kPa
l (31)
32

r 16.75
r
In general, for element 2, we have
xy

�1
U1

1 0 0 0
V1

1 0 X 0 0 0 (32)
A /34 /33

{cr} (1 !
v2) l2�,J
U4

1 -v
=
0 0 --
o/1 o/4
V4

2
'}'4 /34 o/3
r o/ 1 /31
U3

V3

Substituting numerical values into Eq. (32), we obtain

[ 0.3 1
l 0.3

0.91(4)
ts]
210(109 )(10-')
0 0
{a}=

0
0 33)
0 0 0
13.27
X [ 0 0 0 -20 0 2�]
-10 10
2.08
0 -10 -20 10 20
12.19
0.083

Simplifying Eq. (33), we obtain

6964.5
{ ay } = 1 -7.5) kPa (34)
ax

_ 16.1
r xy

09