So far in this chapter, we have applied the finite difference method to steady heat transfer problems. In this section we extend the method to solve transient problems.

We applied the finite difference method to steady problems by discretizing the problem in the space variables and solving for temperatures at discrete points called the nodes. The solution obtained is valid for any time since under steady conditions the temperature do not change with time. In transient problems, however, the temperatures change with time as well as position, and thus the finite difference solution of transient problems requires discretization in time in addition to discretization in space as shown in Figure. 5.16. This is done by selecting a suitable time step and solving for the unknown nodal temperature repeatedly for each until the solution at the desired time is obtained. For example, consider a hot metal object that is taken out of the oven at an initial temperature of T_i at time t=0 and is allowed to cool in ambient air. If a time step of = 5 min is chosen, the determination of the temperature distribution in the metal piece after 3 hours requires the determination of the temperatures times, or in 36 time steps. Therefore, the computation time of this problem will be 36 times that of a steady problem. Choosing a smaller will increase the accuracy of the solution, but it will also increase the computation time.

In transient problems, the superscirpt i is used as the index of counter of times steps, with i=0 corresponding to the specified initial condition. In the case of the hot metal piece discussed above, i=1 corresponds to min, min, and a general time step i corresponds to t_i=i. The notation is used to represent the temperature at the node m at time step i .

Figure. 5.16 Finite Difference Formulation of Time Dependent Problems involves discrete points in time as well as space

The formulation of transient heat conduction problems differs from that of steady ones tin that the transient problems involve an additional term representing the change in the energy content of the medium with time. This additional term appears as a first derivative of temperature with respect to time in the differential equation, and as a change in the internal energy content during in the energy balance formulation. The nodes and the volume elements intransient problems are selected as they are in the steady case, and , again assuming all heat transfer is into the element for convenience, the energy balance on a volume element during a time interval can be expressed as

or

where the heat transfer rate normally consists of conduction terms for interior nodes, but may involve convection, heat flux, and radiation for boundary nodes.

Noting that , where r is the density and C is the specific heat of the element, dividing the relation above by gives

or, for any node m in the medium and its volume element,

where and are the temperatures of node m at times t_i= i and t_i+1=(i+1), respectively, and represents the temperature change of the node during the time interval between the time steps i and i+1.

Note that the ratio is simply the finite difference approximation of the partial derivative that appears in the differential equations of transient problems. Therefore, we would obtain the same result for the finite difference formulation if we followed a strict mathematical approach instead of the energy balance approach used above. Also note that the finite difference formulations of steady and transient problems differ by the

single term on the right side of the equal sign, and the format of that term remains the samein all coordinate systems regardless of whether heat transfer is one, two or three dimensional. For the special case of (ie., no change in temperature with time), the formulation reduces to that of steady case, as expected.

The nodal temperatures in transient problems normally change during each time step, and you may be wondering whether to use temperatures at the previous time step i or the new time step i+1 for the terms on the left side of Equation 5.38. Well, both are reasonable approaches and both are used in practice. The finite difference approach is called the explicit method in the first case and the implict method in the second case, and they are expressed in the general form as

It appears that the time derivative is expressed in forward difference form in the explicit case and backward difference form in the implicit case. Note that both formulations are simply expressions between the nodal temperatures before and after a time interval and are based on determining the new temperatures using the previous temperatures . The explicit and implicit formulations given above are quite general and can be used in any coordinate system regardless of the dimension of heat transfer. The volume elements in multidimensional cases simply have more surfaces and thus involve more terms in the summation.

The explicit and implicit methods have their advantages and disadvantages, and one method is not necessarily better than the other one. Below you will see that the explicit method is easy to implement but imposes a limit on the allowable time step to avoid instabilities in the solution, and the implicit method requires the nodal temperatures to be solved simultaneously for each time step but imposes no limit on the magnitude of the time step. We will limit the discussion to one and two dimensional cases to keep the complexities at a manageable level, but the analysis can readily be extended to three dimensional cases and other coordinate systems.

Transient Heat Conduction in a Plane Wall

Consider transient one dimensional heat conduction in a plane wall of thickness L with heat generation that may vary with time and position and constant conductivity k with a mesh size of D_x= L/M and nodes 0,1,2,… M in the x -direction, as shown in Figure 5.17.

Figure 5.17 The nodal points and volume elements for the transient finite difference formulation of one dimensional conduction in a plane wall

Noting that the volume element of a general interior node m involves heat conduction from two sides and the volume of the element is V_element= A, the transient finite difference formulation for an interior node can be expressed on the basis of Equation 5.38 as

(5.41)

Canceling the surface area A and multiplying by /k , it simplifies to

(5.42)

where is the thermal diffusivity of the wall material. We now define a dimensionless mesh Fourier Number as

(5.43)

Then Equation 5.42 reduces to

(5.44)

Note that the left side of this equation is simply the finite difference formulation of the problem for the steady case. This is not surprising since the formulation must reduce to the steady case for . Also, we are still not committed to explicit or implicit formulation since we did not indicate the time step on the left side of the equation. We now obtain the explicit finite difference formulation by expressing the left side at time step i as

(5.45)

This equation can be solved explicitly for the new temperature (and thus the name explicit method) to give

(5.46)

For all the interior nodes m = 1,2,3,…., M- 1 in a plane wall. Expressing the left side of Equation 5.44 at time step i+1 instead of i would give the implicit finite difference formulation as

(Implicit)

(5.47)

which can be rearranged as

(5.48)

The application of either explicit or implicit formulation above to each of the M-1 interior nodes gives M-1 equations. The remaining two equations are obtained by applying the same method to the two boundary nodes unless, of course, the boundary temperatures are specified as constants (invariant with time). For example, the formulation of the convection boundary condition at the left boundary (node 0) for the explicit case can be expressed as (Figure. 5.18)

Figure. 5.18 Schematic of the explicit finite difference formulation of the convection condition at the left boundary of a plane wall

(5.49)

which simplifies to

(5.50)

Note that in the case of no heat generation and t=0.5, the explicit finite difference formulation for a general interior node reduces to which has an interesting interpretation that the temperature of an interior node at the new time step is simply the average of the temperatures of its neighboring nodes at the previous time step.

Once the formulation (explicit or implicit) is complete and the initial condition is specified, the solution of transient problem is obtained by marching in time using a step size of and determine the nodal temperatures from the initial condition. Taking the initial temperatures as the previous solution at t=0 , obtain the new solution at t= 2using the same relations. Repeat the process until the solution at the desired time is obtained.

Stability Criterion for Explicit Method: Limitation on

The explicit method is easy to use, but is suffers from an undesirable feature that severely restricts utility: the explicit method is not unconditionally stable, and the largest permissible value of the time step is limited by the stability criterion. If the time step is not sufficiently small, the solutions obtained by the explicit method may oscillate wildly and diverge from the actual solution. To avoid such divergent oscillations in nodal temperatures, the value of must be maintained below a certain upper limit established by stability criterion. It can be shown mathematically or by a physical argument based on the second law of thermodynamics that the stability criterion is satisfied if the coefficients of all in the expressions (called the primary coefficients) are greater than or equal to zero for all nodes m. Of course, all the terms involving for a particular node must be grouped together before this criterion is applied.

Different equations for different nodes may result in different restrictions on the size of the time step , and the criterion that is most restrictive should be used in the solution of the problem. A practical approach is to identify the equation with the smallest primary coefficient since it is the most restrictive and to determine the allowable values of by applying the stability criterion to that equation only. A value obtained this way will satisfy the stability criterion for all other equations in the system.

For example, in the case of transient one dimensional heat conduction in a plane wall with specified wall temperatures, the explicit finite difference equations for all the nodes (which are interior nodes ) are obtained from Equation 5.44. The coefficient of in the expression is 1 - 2t , which is independent of the node number m, and thus the stability criterion for all the nodes in this case is 0 or

(5.51)

When the material of the medium and thus its thermal diffusivity a is known and the value of the mesh size is specified, the largest allowable value of the time step can be determined from the above relation. For example, in the case of a brick wall ( a = 0.45 x 10^-6m²/s) with a mesh size of = 0.01m , the upper limit of the time step is

The boundary nodes involving convection and/or radiation are more restrictive than the interior nodes and thus require smaller time steps. Therefore, the most restrictive boundary node should be used in the determination of the maximum allowable time step when a transient problem is solved with the explicit method.

The implicit method is unconditionally stable, and thus we can use any time step we please with that method (of course, the smaller the time step, the better the accuracy of the solution). The disadvantage of the implicit method is that it results in a set of equations that must be solved simultaneously for each time step. Both methods are used in practice.

Two dimensional Transient Heat Conduction

Consider a rectangular region in which heat conduction is significant in the x-and y directions, and consider a unit depth of = 1 in the z direction. Heat may be generated in the medium at a rate of , which may vary with time and position, with the thermal conductivity k of the medium assumed to be constant. Now divide the x-y- plane of the region into a rectangular mesh of nodal points spaced and apart in the x- and y- directions, respectively, and consider a general interior node ( m,n ) whose coordinates are x = m and y = m , as shown in Figure. 5.19.

Noting that the volume element centered about the general interior node ( m,n ) involves heat conduction from four sides (right, left, top and bottom) and the volume of the element is , the transient finite difference formulation for a general interior node can be expressed on the basis of Equation 5.38 as

(5.52)

Taking a square mesh ( = = 1) and dividing each term by k gives after simplifying

(5.53)

where again is the thermal diffusivity of the material and t =a/l² is the dimensionless mesh Fourier number. It can also be expressed in terms of the temperatures at the neighbouring nodes in the following easy-to-remember form:

(5.54)

Figure 5.19 The volume element of a general interior node (m,n) for two dimensional transient conduction in rectangular coordinates

Again the left side of the equation is simply the finite difference formulation of the problem for the steady case as expected. Also, we are still not committed to explicit or implicit formulation since we did not indicate the time step on the left side of the equation. We now obtain the explicit finite difference formulation by expressing the left side at time step i as

(5.55)

Expressing the left side at the time i+1 instead of i would give the implicit formulation. The equation above can be solved explicitly for the new temperature to give

(5.56)

For all the interior nodes ( m,n ) where m = 1,2,3,….. M -1 and n = 1,2,3… N-1 in the medium. In the case of no heat generation and t = 0.25, the explicit finite difference formulation for a general interior node reduces to , which has the interpretation that the temperature of an interior node at the new time step is simply the average of the temperatures of its neighboring nodes at the previous time step.

The stability criterion that requires the coefficient of in the expression to be greater than or equal to zero for all nodes is equally valid for two or three dimensional cases and severely limits the size of the time step that can be used with the explicit method. In the case of transient two dimensional heat transfer in rectangular coordinates, the coefficient of in the expression is 1-4t, and thus the stability criterion for all interior nodes in this is 1-4t > 0, or

(5.57)

Where = = l . When the material of the medium and thus its thermal diffusivity a are known and the value of the mesh size l is specified, the largest allowable value of the time step can be determined from the relation above. Again the boundary nodes involving convection and/or radiation are more restrictive than the interior nodes and thus require smaller time steps. Therefore, the most restrictive boundary node should be used in the determination of the maximum allowable time step when a transient problem is solved with the explicit method.

The application of Equation 5.53 to each of the interior nodes gives equations. The remaining equations are obtained by applying the method to the boundary nodes unless, of course, the boundary temperatures are specified as being constant. The development of the transient finite difference formulation of boundary nodes in two (or three) dimensional problems is similar to the development in the one dimensional case discussed earlier. Again the region is partitioned between the nodes by forming volume elements around the nodes, and an energy balance is written for each boundary node on the basis of Equation 5.38.