TEM – Numerical Solutions
Introduction
In chapters 1, 2 and 3 the basic concepts of heat transfer (conduction, convection and radiation) were introduced and technics to obtain solutions in common situations were presented. After this chapters, students are able:
a) To relate heat fluxes with properties of materials and temperature gradients inside solid bodies,
b) To relate heat fluxes on surfaces of bodies with a body surface temperature, the temperature of the surrounding fluid and the flow over the surface
c) To compute radiative budgets over surfaces knowing the temperature of the surface, its radiation/absorption properties and the temperature of the surrounding surfaces
d) To relate these fluxes knowing that in stationary conditions there is continuity of those eat fluxes when they are sequential (e.g. conduction and convection) or parallel (e.g. convection plus radiation or conduction through different material mounted in a parallel matrix.
Students have also learned the concept of resistivity and how to use this concept to simplify calculations both in case of serial and/or parallel.
Using the energy conservation principle, the evolution equation was obtained for non-stationary conditions, relating the temperature rate of change with the divergence of the heat fluxes and with internal heat generation. This equation was presented into its differential form and was integrated to obtain spatial evolutions in simple geometries and to assess temporal evolutions when the divergence of the heat fluxes is not balanced by heat generated inside the body.
The microscopic concept of thermal conductivity was explained and it was shown that microscopically thermal is identical to viscosity inside a flux or to mass diffusivity. It was also shown that this is the reason why diffusivity units are m2/s that are also the units of the kinematic viscosity of a fluid or of the mass diffusivity and that is also the reason why good heat conductors are also good electric conductors. Based on this concept it was also shown that the evolution equation for velocity, temperature and concentration are identical.
The propagation of heat along a fin was used as an example of application of the energy conservation equation. The equation was integrated using 4 different boundary conditions and assuming that radiation heat transfer is negligible and solutions were obtained for the temperature evolution along the fin and for the heat fluxes. The solutions obtained are able to illustrate the processes responsible for eat transfer and to show the effect of different boundary conditions, but leave room for numerical solutions that are more versatile and if adequately designed can solve heat or mass transfer, i.e. are transversal to the whole course.
Numerical solutions for heat transfer in a fin
A fin is a one-dimensional geometry as represented in Figure 1. Heat flow through its bases, propagates along the fin due to diffusion and is lost to the surrounding fluid by convection and eventually by radiation if the fluid is air and the temperature is high enough for this process to be relevant.
Figure 1: Schematic representation of a fin (adapted from figure 3.17 from the text book [1]).In chapter 3 of the text book [1] the analytical solutions are presented for simple boundary conditions in absence of radiation. Those solutions are obtained starting from a heat budget calculation to the geometry displayed in Figure 2, which shows an elementary volume with infinitesimal length (dx) and crossing the whole fin.
Figure 2: Elementary volume in a generic one dimensional geometry. (adapted from figure 3.17 from the text book [1]).Applying the energy conservation principle for heat which states that the rate “the rate of accumulation of heat inside the control volume is equal to the heat flowing in, minus the heat flowing out, plus the heat generated inside the volume one will get the algebraic equation presented in [1], that making the length of the control volume to shrink to zero gave the differential equation.
Which is valid in case of stationary conditions, no heat sources inside the fin and no radiative heat transfer. If we maintain the algebraic equation, we have to solve it numerically. The number of calculations will require a computer, but we can maintain complexity of the geometry and/or of the boundary conditions.
Algebraic equation (explicit case)
To obtain the differential equation only one control volume was necessary, because the aim was to obtain the budget and then to shrink the control volume to zero and transform differences into differentials. In the numerical approach the aim is to obtain an algebraic equation to each control volume and thus the first step is to divide the geometry into elementary volumes and to compute a budget for each of them.
A 1D geometry with the volumes aligned along one single direction needs only one vector to store the information. The fin can thus be represented schematically as shown in Figure 3. Element (or cell) number 1 is the first on the left and element N is the last one on the right. Element “i” is located between “i-1” and “i+1”. The length of each element is Δx.
Figure 3: schematic representation of a fin using a 1D discretization with N elements (or cells). The first element got the number one and the last got the number N. The element i is located between i-1 and i+1.The geometry definition still need values for the cross section, the thickness and the total fin length:
The conduction heat flux along the fin depends on the thermal conductivity and the accumulation depends on specific the thermal capacity and on specific mass, i.e. on the thermal properties of the material. Convective heat transfer depends on the temperature of the surrounding fluid and on the fluid speed. Radiation budget is a function of the material emissivity and of the incident radiation i.e. on the sky temperature.
Computing heat budgets one gets:
a) Entering conduction eat flux:
In this equation means “evaluated at the interface between cell “i-1” and cell “i”. In case of a rectangular or cylindrical fin and homogeneous material conductivity and area do not change along the fin. Identically for the surface between cell “i” and cell “i+1”.
The convective heat transfer in element “i” is:
Where is the surface area of element “i” is the convection heat transfer coefficient and is the temperature of the fluid surrounding the fin. Identically for the radiative heat transfer:
Where ε is the emissivity and σ is the Stefan-Boltzmann constant and is the temperature of the sky. The rate of accumulation of heat inside the volume is:
Where is the average cross section area of element “i” and ρ is the specific volumetric mass and is the specific heat capacity. The product of the cross section by the length of the element is its volume.
In these equations the header “t” means the value of the property at time “t” and “t+Δt” means the value ne time step latter. The difference between the two values of a property represents the increase in one time-step and if divided by the length of a time step, means the rate of change.
Values of any property at time “t” are known. When the calculation is initiated, those values are called “initial conditions”.
Using the conservation principle one gets:
Dividing by the heat capacity of the element and multiplying by Δt, the time step” the explicit equation for the new value of the temperature is obtained:
This equation gives new values for temperature in aby point knowing their present values, and parameters describing the material, the geometry and the surrounding conditions. In the particular case of uniform cross section and material, one can define some integrated coefficients and use them to write the equation into a more synthetically from, easy to read:
The coefficient “D” is called the Diffusion Number. It is the ratio of the product “time step by the diffusivity” divided by the square of the time step. We will see later that this parameter is numerically very important and that it cannot be larger than 0.5 in any circumstance in case of explicit methods (as the one used in this equation).
Replacing these coefficients into the previous equation one gets:
Or:
Or on a condensed way:
This equation shows that the new temperature in element “i” is a function of the old value at that point, of the neighboring points in the fin (i-1 and i+1) and on convection and on radiation.
Boundary conditions
This equation need as many boundary conditions as the differential equation (two), since the boundary conditions are a consequence of the physics and they both are representing the same physics. One can identify the need for the boundary conditions on the terms due to diffusion, the terms multiplied by “D”.
At “i=1”, the first point, the equation requires the left value to be known (we can designate the left element as element “i=0”. We do not have an equation to compute that value and thus we have to prescribe it. I the analytical assessment we called it i.e. the temperature at the basis of the fin. Here we can do the same:
At the end of the fin we have the same problem. When “i=N” we would need the value of for which we have no equation. We have two possibilities: (a) we prescribe the value of the temperature or we prescribe the value of the heat flux and thus we do not have to compute it. In fact the term responsible for the existence of in the equation is the diffusive heat flux across the right hand side of the control volume. In that case the equation would become:
A simple approach would be to assume null gradient. Another possibility would be to assume that on the right side the eat flux is the convective plus the radiative heat flux computed as any other surface.in that case one would get:
Or,
Initial Conditions
The simulation of a non-stationary process aims at representing the evolution of properties starting at some instant of time and progressing onwards. This means that to compute the solution at the first time step, one needs to know the initial boundary conditions. They can get any physically possible value . As the time evolves the system evolves to a solution that is imposed by the boundary conditions and not by the initial conditions and thus in common engineering situation they can be any[1] (e.g. the environmental temperature).
Solving our algebraic equations
In the explicit case, assuming uniform properties for the material and for the geometry, it is possible to solve the equations using an excel worksheet. The workbook Fin_Numerical.xls has two worksheets already prepared, one for the case of the right side temperature imposed and another for the case of null flux.
These worksheets are illustrative, but aim to support the student on setting up the rationale to solve the problem and to create the ground for a more generic problem involving the resolution of an evolution concentration, whatever is the property transported (heat or mass).
The worksheet heading lists all the variables necessary to solve the equation: (geometry, material parameters, initial conditions and boundary conditions). Initial conditions are written in one line, for time=0 and one for each computation point (element). In this example we are using 10 computation points. Boundary condition specification can require one extra point on the left and another on the right (“0” and “N+1”). At column “0” Tb is imposed. At Column N+1 (i.e. 11 in this case” one can specify a value and use it or not according to the boundary condition selected.
It is recommended to use a worksheet for each boundary condition. This makes easy to store the work.
In any case, inner points are computed as a function of the left and right side values in previous time step and of the atmospheric and sky boundary conditions. They are all computed on the same way and after writing the equation for one of them, to write for other points is just to drag the cell content over the range of cells, both horizontally and vertically, having in mind that parameters are into “frozen cells” in this process.
The left boundary condition is temperature imposed. This is the most reasonable condition, since usually this rod is attached to a body with a higher heat capacity at a fixed temperature (e.g. a pan). With his type of boundary condition, the first cell is computed like inner cells. In case of right imposed temperature, also the last cell (N=10) is computed as others.
The second worksheet shows how to do in case of zero flux, i.e. null temperature gradient at the right face of the last element. The third boundary condition (convective + radiative fluxes at the top of the fin) is not imposed. Please do it yourself.
Computing other properties
Having solved the equation and knowing space and temporal the temperature evolutions of the temperature, one can calculate heat fluxes, e.g. (1) how much heat is being lost through the fin? (2) How is it split into convection and radiation?
One can also study scenarios, e.g. what is the difference between day and night heat fluxes? How does wind speed (h values) change the importance of radiation? How would temperature at the end of the fin change with material thickness, etc.
Plotting the results
Excel is well equipped with graphics tools. It is very easy to visualize the results. Two typical outputs are: (1) spatial distributions for different instants of time and (2) time series for different points. The worksheet shows examples of spatial distributions.
Limitations of this method and steps forward
This methodology has two limitations: the excel worksheet do not allow complex calculations and the explicit methods is limited by stability criteria (not yet addressed).
The excel worksheet is a powerful tool to get results into a simple geometry and under simple boundary conditions and constant parameters. If any of these variables becomes more complex or in case of in when tests of the type “if, then, else” are required, then the excel becomes very limitative and a programing language is necessary. The excel worksheet is impossible to manage implicit calculations.
Implicit calculations are required if stability limitations get important. In this case the algebraic equation transforms into a system of equations, one for each grid point and excel is not designed for that. In this case a program has to be designed using a programing language.
Actually there are many different programing languages. They all share a common programing logics and differ in terms of (a) versatility (the more basic, the most versatile, usually) (b) speed of calculation (c) availability and (d) integration with other languages. VBA is slow (this is not a problem here) but is quite versatile and is freely available. For those reasons and because is embedded into excel it is recommended for this course.
A VBA program will be provided to the students to be used in the resolution of more complex problems involving heat and mass transfer.
Algebraic equation (implicit case)
In the explicit calculation we have assumed that during a time step fluxes /per unit of time) remain constants and equal to their value at the beginning of the time step. Using this approach one could calculate them as a function of the values already known. In the first iteration initial conditions were used and in subsequent iterations the results from previous iterations were used. This approach generated simple algebraic equations where only the new variable at each point was unknown and thus could be explicitly calculated as a function of other variables, all already known. We will see that this approach can be instable, i.e. there are conditions where the system diverges from the correct solution. The other extreme option would be to assume that during a time step the fluxes would remain constant, but equal to their value at the end of the time step. We will see that this approach generates a stable (convergent) algorithm, but it requires more complex programing.
If an implicit algorithm was used, the diffusive flux would be written as:
And identically for the right element face. Convective heat fluxes can also be computed implicitly since they depend on the local temperature on a linear way:
Radiation can also be computed implicitly, linearizing e.g. as:
Doing so, one gets:
That can be rearranged to get:
Where:
The same coefficients that we have got in the explicit case.
Thus, explicit or implicit methods requires the calculation of the same coefficients. The difference is that in the explicit case they multiply by the values of the variables at time “t”, i.e. at the beginning of the time step while in the implicit case they multiply by the values of the variables at the end of the time step and thus the calculation of the new temperature involves the inversion of a tridiagonal matrix.