**MEEG 342: Computer Project I**

** **

**Numerical Solution of Two-Dimensional, Steady-State
Temperature
Distribution in a Rectangular Domain**

__ __

__Purpose__

The purpose of this assignment is to apply the finite difference numerical technique to solve the steady-state heat conduction equation, given by

_{}
(1)

where T is the temperature, x and y are the
physical
coordinates of the domain, g is the volumetric rate of internal heat
generation
(in W/m^{3} ), and k is the thermal conductivity. Note
that such
a heat generation would, for example, be caused by passing electric
current
through the material. Also note that a common mistake is not to
include
the heat generation rate in the proper units (recall that it represents
energy
per unit volume ).

__ __

__Problem Definition__

A two-dimensional solid (the third dimension
is so
long that the temperature does not vary in that direction) generates
heat
internally. One surfaces of the solid is insulated, two of the other
surfaces
are in contact with a body at constant temperatures To and T_{1}
respectively and one surface is cooled with a coolant. We are
interested in the
steady state temperature distribution.

As you know, the solution of Equation (1) gives the distribution of temperature in a Cartesian domain. This expression is a second-order, linear partial differential equation. As in analytical solutions, we need to specify boundary conditions in order to solve the problem. In this assignment, we will use the mixed type of boundary conditions.

Let us consider the domain shown below in
Figure
1. The block is assumed to have a height H and a width W. A
two-dimensional analysis assumes that the domain is either infinitely
deep
(such as a long bar, with temperature independent of the z direction),
or that
it is a finite block where the two surfaces parallel to the x-y plane
are
perfectly insulated.

__The boundary conditions are:__

_{} at y = 0

_{} at y = H

_{}
at x = 0

_{}
at x = W

**Figure 1. Problem
definition
**

__Numerical Solution Approach__

Divide the domain into (m-1) parts in the x-direction and (n-1 ) parts in the y-direction. Discretization of the domain requires m nodes by n nodes for the mesh size, requiring a dimension statement for the temperature of real T(m,n). The individual nodes are numbered (i,j), where i=1,2,3,.....m and j=1,2,3,....n, related as shown below in Figure 2.

**Figure 2. Grid cell
notation**

In addition, we will use the following for the width, height, and boundary temperatures

W = 60 cm
H = 30 cm To=
160C T_{1}
= 100 C Tamb = 20 C

**SOLUTIONS**

*Temperature solution for all cases**Plot the temperature solution for all cases.*

__ __

__Case 1a: m=21, n=11, k=2 W/m C
, g=0,
h=0 __

__Case 1b: m=41, n=21, k=2 W/m C , g=0, h =
0__

__Case 2a: m=21, n=11, k=2 W/m C
, g=0, h
= 500 W/m ^{2} C__

__Case 2b: m=41, n=21, k=2 W/m C , g=0, h =
500 W/m ^{2}
C__

__Case 3a: m=21,
n=11, k= 50
W/m C , g=0.9 MW/m ^{3} , h=0 __

__Case 3b: m=41,
n=21, k= 50
W/m C , g=0.9MW/m ^{3} , h=0 __

__ __

__ __

__ __

__ __

__Case 4a: m=21,
n=11, k=50
W/m C , g=0.9 MW/m ^{3}, h = 500 W/m^{2} C __

__Case 4b: m=41, n=21, k=50 W/m
C ,
g=0.9 MW/m ^{3}, h = 500 W/m^{2} C __

** **

*For Case 1 as h=0 it corresponds to insulated boundary condition at the top. As the bottom is also insulated the temperature should vary only with x and should be the same along the y direction.*

_{}

Since T is not a function of x and y, the partial differential equation turns to an ordinary differential equation. And with g = 0;

_{}

where *a, b* are
constants. From BC’s, _{}
at x = 0 and _{}
at x = W, the analytical solution is

_{}

With numerical values

_{}

*Calculate q that needs
to be supplied along the left wall to maintain the wall temperature at
T0
analytically and compare it with the calculations of q for case 1a and
1b. Repeat
the procedure at the right wall for this case. Remember for Case 1, the
q
supplied through the left wall should be leaving through the right wall
under
steady state conditions with no heat generation. *

And the heat flux from left to right

_{}

_{}

With numerical values for case 1a, at the left the heat flux for unit depth is

_{}

_{}

For case 1b with the same formulas, _{} and _{}

For case b, the results are closer to the analytical solution.

The error for the left wall decreases from %31 to %26 and for the right wall, from %5 to %3.

* *

* *

* *

*Compare the
temperature solution for case 1a and 1b with analytic solution for
the
temperature in the x-direction for steady state and no heat
generation.*

*Case 1a*

*Case 1b*

As it seen from the figures, although the numerical solution is very close to the analytical one for both cases, in the second one the accuracy is better.

*Compare the temperature solution of case 2a with 2b. Which one is more accurate? To prove your solution is correct conduct a global energy balance for case 2.*

The Energy balance,

_{}

_{} and *g
= 0*

For unit depth and Fourier law at the left and the right wall,

_{}

For the numerical solution, one can get the heat flux by using the sum,

_{}
_{}

* *

*Calculate the error in
the energy balance for case 2a and 2b.*

*Case 2a:*

_{}

* *

*Case 2b:*

_{}

In the second case, better accuracy is obtained. This is simply a result of finer mesh.

*Case 3 is similar to case 1 except in this case there is internal heat generation.*

*Find
the analytical solution for the temperature along the x direction and
compare
it with the numerical solutions obtained for cases 3a and 3b. Also
calculate
the q along the left and right wall and compare it with the analytic
solution.*

_{}

Since T is not a function of x and y, the partial differential equation turns to an ordinary differential equation. And with nonzero g term;

_{}

where *b, c*
are constants. From BC’s, _{}
at x = 0 and _{}
at x = W, the constants can be found. _{} and _{} then,

_{}

With numerical values

_{}

And the heat flux for unit depth

_{}

At the left and the right wall the analytical solution of heat flux is,

_{}

_{}

With numerical values for case 3a, at the left the heat flux is,

_{}

_{}

And for case 3b,

_{}

_{}

The errors are found as follows, *Case
a*

* Case a @right, Error =
%3.1*

* Case b*

* Case b @right, Error =
%3.0.*

For case b, the error decrease from %5 to %3.

For both cases the analytical and the numerical solution are same and the temperature distribution curves fit.

*For Case 4, compare the temperature solution of cases 4a with 4b. Conduct a global energy balance for this case. Calculate the error in the energy balance and confirm that it goes down, as the mesh is refined.*

By using the same method in question 4 for case 2 and this time adding the heat generation term to the energy balance,

The Energy balance,

_{}

_{} * *

For unit depth and Fourier law at the left and the right wall,

_{}

For the numerical solution, one can get the heat flux by using the partial sum which is actually the discretization of the integral,

_{} _{} _{}

* *

*Calculate the error in
the energy balance for case 4a and 4b.*

* *

*Case 4a:*

_{}

* *

*Case 4b:*

_{}

The effect of finer mesh can be seen here from these results. In the second case, the result is closer to zero.

*Discussions comparing the results of cases a with cases b, should be included. The discussion should focus on the convergence of the numerical solution and the effect of changing properties on the temperature distribution. Will any of the materials start melting?*

The accuracy of a solution in finite difference method is related with the number of meshes used in the discretization. Therefore, in cases b the accuracy is better than those of in cases a.

Convergence of a numerical solution depends on the convergence accuracy defined in the sample code. Decreasing the convergence accuracy leads to more accurate results.

The effect of thermal conductivity coefficient can be seen by comparing the first two and the last two cases. When k is increased from 2 to 50, heat transfer increases. The heat generation leads to higher temperatures and higher thermal conductivity is required in order to transfer heat from the system and avoid high temperatures that might cause melting of the material.

**Extra Notes…**

Derivation of the equations by using finite difference method for the top boundary:

The Energy balance,

_{}

_{}

With *dx*=*dy*=_{}
and for unit depth,

_{}

Since we have half of the cell,

_{}

_{} in the figure above.
Divide by k all terms,

_{}

Rearrange the equation by multiplying all terms
with _{}

_{}

_{}

**Sample Code**

'MEEG302 Heat Transfer

'example program for computer assignment 1

'this program calculates the temperature distribution in a simple

'rectangular domain with constant heat flux at the top and constant

'temperature at the other three sides

Sub sample()

'define constants

Const M = 21 'grid size

Const N = 11 '

Const X = 0.6 'rectangular dimensions (m)

Const Y = 0.3 '

Const g = 900000 'heat generation

Const h = 0 'convection coefficient

Const T0 = 160 'constant temperature at side : left

Const T1 = 100 'constant temperature at side : right

Const Ts = 20 'surrounding temp

Const EPS = 0.01 'convergence accuracy

Const k = 50 'conduction coeff

'define arrays and variables

Dim Tnew(M, N), T(M, N) As Single 'temperatures

Dim isConverged As Boolean 'flag

Dim dx, dy As Single 'delta x and delta y

'initialize variables

dx = X / (M - 1)

dy = Y / (N - 1)

'initialize the temperature array

For i = 1 To M

For j = 1 To N

T(i, j) = 125

Tnew(i, j) = 125

Next j

Next i

'**** main loop

'this loop keeps recalculating the temperature field until

'it has converged

Do

'Tnew() holds the newly calculated temperatures

'T() holds the values from the last iteration

'top boundary (convection)

For i = 2 To M - 1

j = N

Tnew(i, j) = (1 / (2 * (2 + h * dx / k))) * (2 * T(i, j - 1) + T(i - 1, j) + T(i + 1, j) + 2 * h * dx * Ts / k + g * dx * dy / k)

Next i

'left boundary (constant temp)

For j = 1 To N

i = 1

Tnew(i, j) = T0

Next j

'right boundary (constant temp)

For j = 1 To N

i = M

Tnew(i, j) = T1

Next j

'bottom boundary (constant temp)

For i = 2 To M - 1

j = 1

Tnew(i, j) = (1 / 4) * (2 * T(i, j + 1) + T(i - 1, j) + T(i + 1, j) + g * dx * dy / k)

Next i

'interior nodes

For i = 2 To M - 1

For j = 2 To N - 1

Tnew(i, j) = (1 / 4) * (T(i - 1, j) + T(i + 1, j) + T(i, j - 1) + T(i, j + 1) + g * dx * dy / k)

Next j

Next i

'check to see if the temperatures have converged

isConverged = True

For i = 1 To M

For j = 1 To N

If Abs(T(i, j) - Tnew(i, j)) > EPS Then isConverged = False

Next j

Next i

'copy new temperatures to old temperatures and display on the worksheet

For i = 1 To M

For j = 1 To N

'this command displays the elements of T on the worksheet

'the command is Cells(row,column)

'because Excel counts rows from the top the grid is

'actually printed upside down

Cells(j + 1, i + 1) = Tnew(i, j)

'copy new temperatures to old array T

T(i, j) = Tnew(i, j)

Next j

Next i

'stop the loop if the values have converged

Loop While (isConverged = False)

End Sub