1.

A square plate domain of side 8 cm and thickness 1 cm has a concentrated point heat
source of 200 W at the center. All the sides are exposed to air medium at a far field
temperature of 200C. Using symmetry of the problem and formulating only one finite
element, determine the temperature at the center and at the corners. Take k = 50
W/m.K and h = 500 W/m2.K.

Solution
Considering symmetrical condition given for the problem, we can consider a single triangular
element as shown in Figure.

Coordinates of each node and corresponding coefficients of the linear interpolation functions are
given as

x1= 0 y1 = 0 x2= 4 y2 = 0 x3 = 4 y3 = 4
1
∴ 𝐴= × 4 × 4 = 8 𝑐𝑚2 and ∆= 2𝐴 = 16 𝑐𝑚2
2
𝑥2 𝑦3 − 𝑥3 𝑦2 16 − 0 𝑥3 𝑦1 − 𝑥1 𝑦3 0 − 0
𝑎1 = = =1 𝑎2 = = =0
∆ 16 ∆ 16
𝑥1 𝑦2 − 𝑥2 𝑦1 0 − 0
𝑎3 = = =0
∆ 16
𝑦2 − 𝑦3 0 − 4 1 𝑦3 − 𝑦1 4 − 0 1
𝑏1 = = =− 𝑏2 = = =
∆ 16 4 ∆ 16 4
𝑦1 − 𝑦2 0 − 0
𝑏3 = = =0
∆ 16
𝑥3 − 𝑥2 4 − 4 𝑥1 − 𝑥3 0 − 4 1
𝑐1 = = =0 𝑐2 = = =−
∆ 16 ∆ 16 4
𝑥2 − 𝑥1 4 − 0 1
𝑐3 = = =
∆ 16 4

Verify that 𝑎𝑖 = 1, 𝑏𝑖 = 0 and 𝑐𝑖 = 0.

1
 Then stiffness matrix for the element due to the effect of conductivity is given by
1 2 3
2 2
𝑏1 + 𝑐1 𝑏1 𝑏2 + 𝑐1 𝑐2 𝑏1 𝑏3 + 𝑐1 𝑐3 1
𝑘
𝐾 = 𝑘𝐴𝑡 𝑏1 𝑏2 + 𝑐1 𝑐2 𝑏2 2 + 𝑐2 2 𝑏2 𝑏3 + 𝑐2 𝑐3 2
𝑏1 𝑏3 + 𝑐1 𝑐3 𝑏2 𝑏3 + 𝑐2 𝑐3 𝑏3 2 + 𝑐3 2 3
1 2 3
1 1 1
− 0
16 16
1 2 1 2
= 50 × 10−2 × 8 × 1 − −
16 16 16
1 1
0 −
16 16 3
1 2 3
1 1 1
− 0
4 4
1 2 1 2
= − −
4 4 4
1 1
0 −
4 4 3

Where k = 50  10-2 W/cm.K, A= 8 cm2 and t= 1cm.

Then stiffness matrix for the element due to the effect of convective boundary condition on edge
2-3 is given by
1 2 3
0 0 0 1
1 1
𝐾 ℎ 0
= ℎ𝑇∞ 𝑡 3 6 2
1 1
0
6 3 3
1 2 3
0 0 0 1
1 1
0
= 500 × 10−4 × 20 × 1 3 6 2
1 1
0
6 3 3
1 2 3
0 0 0 1
100 100
0
= 3 6 2
100 100
0
6 3 3

2
 Where h = 500  10-4 W/cm2.K, T = 20 and t= 1cm.

Then global stiffness matrix for the element is given by

𝐾 𝑘 ℎ
= 𝐾 + 𝐾
1 2 3
1 1
− 0 1
4 4
1 203 197
= − 2
4 6 12
197 403
0 3
12 12

Force vector due to heat generation is given by

200 1 Heat generated by the point source is distributed to total eight
𝑓 𝑄 = 8 elements of the complete domain. This contributes to node 1
0 2 only.
0 3

Similarly, force vector due to the effect of convective boundary condition on edge 2-3 is given
by
0 1 0 1
1 1 0 1
ℎ −4
𝑓 = ℎ𝑇∞ 𝑡𝑠23 2 2 = 500 × 10 × 20 × 1 × 4 2 2 = 200 2
1 1 200 3
2 3 2 3

Where s23 is the length of edge 2-3.

There is not any contribution on force vector due to heat flux on any edges [𝑞12 = 𝑞31 = 0].
Therefore, global force vector is then given by
𝑓 𝑄 ℎ
= 𝑓 + 𝑓
25 1
= 200 2
200 3

Hence the equation for the unknown temperatures is then given by

𝐾 𝑇 = 𝑓

3
 or, 1 1
− 0 𝑇1 25
4 4
1 203 197
− 𝑇2 = 200
4 6 12
197 403
0 𝑇3 200
12 12

Solving, we get
𝑇1 = 104.980 𝐶, 𝑇2 = 4.980 𝐶 and 𝑇3 = 3.520 𝐶.

Hence the temperature at the center of the plate is 104.980 𝐶 and each corner is 3.520 𝐶.

2. A triangular domain with unit thickness is shown in Figure P.6. The domain is under
steady state heat conduction with side 2-3 subjected to a fluid at ambient temperature
250C and the nodes 1 and 2 are maintained at 1000C and 2000C respectively. Determine
the temperature at the centroid of the triangle. Take k = 50 W/m.K and h = 500
W/m2.K.

Solution
Coordinates of each node and corresponding coefficients of the linear interpolation functions are
given as

x1= 0 y1 = 0 x2= 4 y2 = 0 x3 = 2 y3 = 2
1
∴ 𝐴= ×4×2=4 and ∆= 2𝐴 = 8
2
𝑥2 𝑦3 − 𝑥3 𝑦2 8 − 0 𝑥3 𝑦1 − 𝑥1 𝑦3 0 − 0
𝑎1 = = =1 𝑎2 = = =0
∆ 8 ∆ 8
𝑥1 𝑦2 − 𝑥2 𝑦1 0 − 0
𝑎3 = = =0
∆ 8
𝑦2 − 𝑦3 0 − 2 1 𝑦3 − 𝑦1 2 − 0 1
𝑏1 = = =− 𝑏2 = = =
∆ 8 4 ∆ 8 4
𝑦1 − 𝑦2 0 − 0
𝑏3 = = =0
∆ 8

4
 𝑥3 − 𝑥2 2 − 4 1 𝑥1 − 𝑥3 0 − 2 1
𝑐1 = = =− 𝑐2 = = =−
∆ 8 4 ∆ 8 4
𝑥2 − 𝑥1 4 − 0 1
𝑐3 = = =
∆ 8 2

Then stiffness matrix for the element due to the effect of conductivity is given by
1 2 3
2 2
𝑏1 + 𝑐1 𝑏1 𝑏2 + 𝑐1 𝑐2 𝑏1 𝑏3 + 𝑐1 𝑐3 1
𝑘
𝐾 = 𝑘𝐴 𝑏1 𝑏2 + 𝑐1 𝑐2 𝑏2 2 + 𝑐2 2 𝑏2 𝑏3 + 𝑐2 𝑐3 2
𝑏1 𝑏3 + 𝑐1 𝑐3 𝑏2 𝑏3 + 𝑐2 𝑐3 𝑏3 2 + 𝑐3 2 3
1 2 3
1 1 1
0 −
8 8
1 1 2
= 50 × 4 0 −
8 8
1 1 1
− −
8 8 4 3
1 2 3
25 0 −25 1
= 0 25 −25 2
−25 −25 50 3

Where k = 50 W/m.K, A = 4 m2.

Then stiffness matrix for the element due to the effect of convective boundary condition on edge
2-3 is given by

1 2 3
0 0 0 1
1 1
𝐾 ℎ 0
= ℎ𝑇∞ 3 6 2
1 1
0
6 3 3
1 2 3
0 0 0 1
1 1
0
= 500 × 25 3 6 2
1 1
0
6 3 3

5
 1 2 3
0 0 0 1
12500 6250
0
= 3 3 2
6250 12500
0
3 3 3

Where h = 500 W/m2.K, T= 25.

Then global stiffness matrix for the element is given by

𝐾 𝑘 ℎ
= 𝐾 + 𝐾
1 2 3
25 0 −25 1
12575 6175
0
= 3 3 2
6175 12650
−25
3 3 3

Force vector due to the effect of convective boundary condition on edge 2-3 is given by
0 1 0 1
1 1 0 1
ℎ = ℎ𝑇∞ 𝑠23 = 500 × 25 × 1 × 2 2 2 = 17677.67 2
𝑓 2 2 2
1 1 17677.67 3
2 3 2 3

Where s23 is the length of edge 2-3.

There is not any contribution on force vector due to heat flux on any edges [𝑞12 = 𝑞31 = 0] as
well as heat generation. Therefore, global force vector is then given by
𝑓 ℎ
= 𝑓
0 1
= 17677.67 2
17677.67 3

Hence the equation for the unknown temperatures is then given by

𝐾 𝑇 = 𝑓
or, 25 0 −25 𝑇1 𝑄
12575 6175
0
3 3 𝑇2 = 17677.67
6175 12650
−25
3 3 𝑇3 17677.67

Substituting 𝑇1 = 1000 𝐶, 𝑇2 = 2000 𝐶 and solving for 𝑇3 , we get

6
 𝑇3 = −92. 840 𝐶.

Then temperature variation within the element is given by

𝑇 = 𝑁1 𝑇1 + 𝑁2 𝑇2 + 𝑁3 𝑇3
= (𝑎1 + 𝑏1 𝑥 + 𝑐1 𝑦) 𝑇1 + (𝑎2 + 𝑏2 𝑥 + 𝑐2 𝑦) 𝑇2 + (𝑎3 + 𝑏3 𝑥 + 𝑐3 𝑦) 𝑇3
𝑥 𝑦 𝑥 𝑦 𝑦
= 1 − − 𝑇1 + − 𝑇2 + 𝑇
4 4 4 4 2 3
𝑥 𝑦 𝑥 𝑦 𝑦
= 1− − × 100 + − × 200 + × −92.84
4 4 4 4 2
= 25𝑥 − 121.42𝑦 + 100

Coordinates of the centroid of the triangle is

0+4+2 0+0+2 2
𝑥= =2 𝑦= =
3 3 3

Then temperature at the centroid of the triangle is then given by

𝑇(𝑥, 𝑦) = 25𝑥 − 121.42𝑦 + 100
2
= 25 × 2 − 121.42 × + 100
3
= 69.050 𝐶

3. The nodal values of the element, as shown in Figure P.6, in a heat conduction problem
are: T4 = 1200C, T9 = 2100C, T10 = 1500C.

Determine where the 1600C isotherm cuts the boundaries of the element. Also
determine the [K] matrix for the element. Take k =1 W/m.K.

Solution
Coordinates of each node and corresponding coefficients of the linear interpolation functions are
given as

7
 x4= 2 y4 = 4 x9= 8 y9 = 10 x10= 2 y10 = 10
1
∴ 𝐴= × 6 × 6 = 18 and ∆= 2𝐴 = 36
2
𝑥9 𝑦10 − 𝑥10 𝑦9 80 − 20 5 𝑥10 𝑦4 − 𝑥4 𝑦10 8 − 20 1
𝑎4 = = = 𝑎9 = = =−
∆ 36 3 ∆ 36 3
𝑥4 𝑦9 − 𝑥9 𝑦4 20 − 32 1
𝑎10 = = =−
∆ 36 3
𝑦9 − 𝑦10 10 − 10 𝑦10 − 𝑦4 10 − 4 1
𝑏4 = = =0 𝑏9 = = =
∆ 36 ∆ 36 6
𝑦4 − 𝑦9 4 − 10 1
𝑏10 = = =−
∆ 36 6
𝑥10 − 𝑥9 2 − 8 1 𝑥4 − 𝑥10 2 − 2
𝑐4 = = =− 𝑐9 = = =0
∆ 36 6 ∆ 36
𝑥9 − 𝑥4 8 − 2 1
𝑐10 = = =
∆ 36 6

Then the shape functions for the triangular element are given as
5 𝑦
𝑁4 = 𝑎4 + 𝑏4 𝑥 + 𝑐4 𝑦 = −
3 6
1 𝑥
𝑁9 = 𝑎9 + 𝑏9 𝑥 + 𝑐9 𝑦 =− +
3 6
1 𝑥 𝑦
𝑁10 = 𝑎10 + 𝑏10 𝑥 + 𝑐10 𝑦 =− − +
3 6 6

Then temperature variation within the element is given by

𝑇 = 𝑁4 𝑇4 + 𝑁9 𝑇9 + 𝑁10 𝑇10
= (𝑎4 + 𝑏4 𝑥 + 𝑐4 𝑦) 𝑇4 + (𝑎9 + 𝑏9 𝑥 + 𝑐9 𝑦) 𝑇9 + (𝑎10 + 𝑏10 𝑥 + 𝑐10 𝑦) 𝑇10
5 𝑦 1 𝑥 1 𝑥 𝑦
= − 𝑇4 + − + 𝑇9 + − − + 𝑇10
3 6 3 6 3 6 6
5 𝑦 1 𝑥 1 𝑥 𝑦
= − × 120 + − + × 210 + − − + × 150
3 6 3 6 3 6 6
= 10𝑥 + 5𝑦 + 80

Hence the equation for the 1600C isotherm is given by

160 = 10𝑥 + 5𝑦 + 80
∴ 10𝑥 + 5𝑦 = 80 ……….(1)

8
 It can be understood from the given nodal temperatures that 1600C isotherm will intersect edges
4-9 and 4-10 only.

Equation of edge 4-9 is then given by

−𝑥 + 𝑦 = 2 ……….(2)

Point where the 1600C isotherm intersect the edge 4-9 is given by the solution of simultaneous
equations (1) and (2), i.e.,
14 20
𝑥= 𝑦=
3 3

Similalrly equation of edge 9-10 is then given by

𝑦 = 10 ……….(3)

Point where the 1600C isotherm intersect the edge 4-9 is given by the solution of simultaneous
equations (1) and (3), i.e.,
𝑥=3 𝑦 = 10

Then stiffness matrix for the element due to the effect of conductivity is given by
4 9 10
2 2
𝑏4 + 𝑐4 𝑏4 𝑏9 + 𝑐4 𝑐9 𝑏4 𝑏10 + 𝑐4 𝑐10 4
𝑘
𝐾 = 𝑘𝐴 𝑏4 𝑏9 + 𝑐4 𝑐9 𝑏9 2 + 𝑐9 2 𝑏9 𝑏10 + 𝑐9 𝑐10 9
𝑏4 𝑏10 + 𝑐4 𝑐10 𝑏9 𝑏10 + 𝑐9 𝑐10 𝑏10 2 + 𝑐10 2 10
4 9 10
1 1 4
0 −
36 36
1 1 9
= 1 × 18 0 −
36 36
1 1 1
− −
36 36 18 10
4 9 10
1 1 4
0 −
2 2
1 1 9
= 0 −
2 2
1 1
− − 1 10
2 2