This section provides some preliminary information about heat transfer equation, also known as the diffusion equation.

Return to computing page for the first course APMA0330
Return to computing page for the second course APMA0340
Return to Mathematica tutorial for the first course APMA0330
Return to Mathematica tutorial for the second course APMA0340
Return to the main page for the first course APMA0330
Return to the main page for the second course APMA0340
Return to Part VI of the course APMA0340
Introduction to Linear Algebra with Mathematica

Heat Conduction Equations

Animation in Mathematica:

The laws of thermodynamics tell us that no matter what the temperature distribution of the object is initially, the system must undergo a process that brings the object to thermal equilibrium. This process must obey the heat equation. This chapter deals with heat transfer processes that occur in solif matters without bulk motion of the matter. A solid (a block of metal, say) has one surface at a high temperature and one at a lower temperature. This type of heat conduction can occur, for example,through a turbine blade in a jet engine. The outside surface, which is exposed to gases from the combustor, is at a higher temperature than the inside surface, which has cooling air next to it. The level of the wall temperature is critical for a turbine blade.

The concept of heat flux was a key contribution of Joseph Fourier, in the analysis of heat transfer phenomena. According to the transport definition, flux may be a single vector, or it may be a vector field / function of position. In the latter case flux can readily be integrated over a surface. So flux is defined as the rate of flow of a property per unit area, which has the dimensions [quantity]·[time]−1·[area]−1. The area is of the surface the property is flowing "through" or "across". For example, the magnitude of a river's current, i.e. the amount of water that flows through a cross-section of the river each second, or the amount of sunlight energy that lands on a patch of ground each second, are kinds of flux. We denote the heat flux as

\[ {\bf q} = -\kappa \,\frac{\partial T}{\partial {\bf n}} \quad \left[ \mbox{with units: } \frac{W}{m^2} \right] , \]
where κ is the thermal conductivity, n is the outward normal to the surface. We also need the heat rate
\[ Q = {\bf q}\,A \quad \left[ \mbox{with units: } W \right] , \]
where A is the cross-section area. The heat equation in rectangular coordinates:
\[ \rho c\,\frac{\partial T}{\partial t} = \frac{\partial}{\partial x} \left( \kappa \,\frac{\partial T}{\partial x}\right) + \frac{\partial}{\partial y} \left( \kappa \,\frac{\partial T}{\partial y}\right) + \frac{\partial}{\partial z} \left( \kappa \,\frac{\partial T}{\partial z}\right) + f(x,y,z,t) . \]
For constant coefficients, we get the diffusion (or heat transfer) constant coefficient equation)
\[ \frac{\partial T}{\partial t} = \frac{\kappa}{\rho\,c} \,\nabla^2 T = \frac{\kappa}{\rho\,c} \left( \frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} + \frac{\partial^2 T}{\partial z^2} \right) . \]
The differential operator
\[ \Delta = \nabla^2 = \frac{\partial^2}{\partial x^2_1} +\frac{\partial^2}{\partial x^2_2} + \cdots + \frac{\partial^2}{\partial x^2_n} \]
is called the Laplacian in n-dimensional Euclidean space n. When a spherical body is homogeneous, isotropic and opaque, the thermal conductivity κ is temperature-dependent, while the density ρ and specific heat c are assumed to be constant, the heat equation becomes
\[ \rho c\,\frac{\partial T}{\partial t} = \frac{1}{r^2} \,\frac{\partial}{\partial r} \left[ r^2 \kappa (T)\, \frac{\partial T}{\partial r} \right] , \qquad \mbox{for}\quad 0 \le r < R, \quad t> 0 , \]
The diffusion or heat transfer equation in cylindrical coordinates is
\[ \frac{\partial T}{\partial t} = \frac{1}{r} \,\frac{\partial}{\partial r} \left( r\alpha\, \frac{\partial T}{\partial r} \right) . \]
Consider transient convective process on the boundary (sphere in our case):
\[ - \kappa (T)\, \frac{\partial T}{\partial r} = h \left( T - T_{\infty} \right) \qquad \mbox{at} \quad r = R. \]
If a radiation is taken into account, then the boundary condition becomes
\[ - \kappa (T)\, \frac{\partial T}{\partial r} = h \left( T - T_{\infty} \right) \epsilon\sigma \left( T^4 - T^4_{s} \right) \qquad \mbox{at} \quad r = R. \]
Here T is an environment fluid (or air) temperature (usually constant), Ts is the sink temperature, T denotes the temperature, t the time, r the spatial coordinate, h the convective heat transfer coefficient, ϵ the surface emissivity (0 ≤ ϵ ≤ 1), and σ = 5.67×10-8 (W/(m²K4)) the Stefan–Boltzmann constant, respectively.


One dimensional Heat Transfer Equation in infinite strip

The one dimensional heat conduction equation

\[ u_t = \alpha\, u_{xx} \qquad\mbox{or} \qquad \frac{\partial u}{\partial t} = \alpha\,\frac{\partial^2 u}{\partial x^2} , \]
where \( \alpha = \kappa/(\rho c_p) \) is a constant known as the thermal diffusivity, κ is the thermal conductivity, ρ is the density, and cp is the specific heat of the material in the bar. You can also change the border style to thick, dotted You can also change the border style to thick, dotted
Metals   Ag     Cu     Al     Fe   Steel
κ [W/m-K] 420 390 200 70 50

Non-metals   H2O     Air   Engine oil   H2     Brick     Wood     Cork  
κ [W/m-K] 0.6 0.026 0.15 0.18 0.6 0.2--0.5 0.04

To make the solution more meaningful and simpler, we group as many physical constants together as possible. Therefore, it is convenient to introduce dimensionless variables.

Convective heat transfer occurs whenever an object is either hotter or colder than a surrounding gas or fluid. Convection is a process of heat transfer between a surface and a fluid flowing in contact with it, never in a solid. It occurs always within a thin stagnant fluid film layer on the heat transfer surface. Heat transfer by convection is more difficult to analyze than heat transfer by conduction because no single property of the heat transfer medium, such as thermal conductivity, can be defined to describe the mechanism. Convective heat transfer is complicated by the fact that it involves fluid motion as well as heat conduction. Heat transfer by convection varies from situation to situation (upon the fluid flow conditions), and it is frequently coupled with the mode of fluid flow.

The primary equation for the rate of convection heat transfer is known as Newton's Law of Cooling:

\[ Q = hA \left( T_{\infty} - T \right) , \]
where q is the convective heat transfer rate (units: W), h is the convective heat transfer coefficient (in units W/(m²K), A (units: m²) is the surface area of the object being cooled or heated, T is the bulk temperature of the surrounding gas or fluid, and T is the surface temperature (units: K) of the object. The ratio q/A is the heat flux. The algebraic sign of Newton's Law of Cooling is positive for T∞ > T (heat transfer into the object) and negative when T∞ < T (heat transfer out of the object). The convective heat transfer coefficient h is usually a positive, experimentally determined value. It depends upon the surface geometry, the nature of the fluid motion, the properties of the fluid, and the bulk fluid velocity.

There are two types of convection: natural convection and forced convection. Natural convection is produced by density differences in a fluid due to temperature differences (e.g., as in “hot air rises”). Global atmospheric circulation and local weather phenomena (including wind) are due to convective heat transfer. In forced convection, the fluid motion is generated by a source like a pump or a fan. It is one of the main forms of heat transfer used by engineers, because large amounts of thermal energy can be transported efficiently. It is used in heating and air conditioning systems, electronics cooling, and in numerous other technologies. Forced convection is used in designing heat exchangers, in which one fluid stream is used to heat or cool another fluid stream.

One of the most common examples of convection is natural convection in a non-mechanically ventilation/air conditioned building. As people enter a building, the lights get turned on and the sun heats the building, the air in the building begins to get warmer. The warm air is less dense than the air around it and begins to rise up and out of the building. The empty space left by the warm air is then replaced by cooler outside air and the cycle continues. This convective heat transfer through the movement of air is called natural convection. It is referred to as natural because it does not rely on a mechanical source like a fan to move the air.

Let u(x, t) denote the solution of the heat equation subject to the initial condition:

\[ \frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2} , \qquad u(x,0) = f(x) , \]
where f(x) is a given function on the interval \( (-\infty , \infty ) . \) Its solution is
\[ u(x, t) = \left( 4\pi t \right)^{-1/2} \int_{-\infty}^{\infty} f(\xi )\, e^{-(x-\xi )^2 /(4t)} \,{\text d}\xi , \]
or, writing
\[ \eta = (x-\xi ) /\sqrt{4t} \qquad \Leftrightarrow \qquad \xi = x - \eta \sqrt{4t} , \]
\[ g(x, \eta ) = f(\xi ) = f \left( x - \eta \sqrt{4t} \right) , \]
we obtain
\[ u(x, t) = \pi^{-1/2} \int_{-\infty}^{\infty} g(x, \eta )\, e^{-\eta^2} \,{\text d}\eta = E \left[ g(x, \eta ) \right] , \]
where η is a Gaussian distributed random variable, i.e., the probability that η lie between y and \( y + {\text d}y \) is
\[ \Pr [ y < \eta < y + {\text d}y ] = \pi^{-1/2} \, e^{-y^2} \,{\text d}y , \]
and \( E \left[ g(x, \eta ) \right] \) denotes the expected value of \( g(x, \eta ) . \) Assume that \( g(x, \eta )\, e^{-\eta^2 /2} \) is a square integrable on the line, then we can expand it in Hermite series as a function of η:
\[ g(x, \eta ) \, e^{-\eta^2 /2} = \sum_{n\ge 0} a_n \, H_n (\eta )\, e^{-\eta^2 /2} , \]
\[ a_n = \pi^{-1/2} \, \int g(x, \eta ) \, e^{-\eta^2 /2}\, H_n (\eta )\,{\text d}\eta = E \left[ H_n (\eta )\,g(x, \eta ) \right] . \]
Note that \( a_0 = u(x,t) . \)

One dimensional Heat Equation for a finite rod

3D heat equation ??? to be fixed
Monitor[solution = NDSolveValue[{D[T[x, y, z, t], t] == b*Laplacian[T[x, y, z, t], {x, y, z}], T[x, y, z, 0] == If[0 <= y <= 3 && 0 <= z <= 3, 200, 25], DirichletCondition[ T[x, y, z, t] == 25, {x == 0 || x == 15 || z == 0 || y == 0}]}, T, Element[{x, y, z}, region], {t, 0, 20}, EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])], monitor];
Visualize the solution at t=20. First get this right before you think about animations.
SliceContourPlot3D[solution[x, y, z, 20], Element[{x, y, z}, region]]
2D heat equation:
RT = 1 R = 1 FUNC = NDSolve[{D[T[x, y, t], t] == RT*(D[T[x, y, t], x, x] + D[T[x, y, t], y, y]), T[x, y, 0] == 0, T[0, y, t] == R*t, T[9, y, t] == R*t, T[x, 0, t] == R*t, T[x, 9, t] == R*t}, T, {x, 0, 9}, {y, 0, 9}, {t, 0, 6}]; a = Table[ Plot3D[T[x, y, t] /. FUNC, {x, 0, 9}, {y, 0, 9}, Mesh -> 15, PlotRange -> {{0, 9}, {0, 9}, {-1, 10}}, ColorFunction -> Function[{x, y, z}, Hue[.3 (1 - z)]]], {t, 0, 6}] Export["c:\anim.gif", a]
ListAnimate[ Table[Plot3D[T[x, y, t] /. FUNC, {x, 0, 9}, {y, 0, 9}, Mesh -> 9], {t, 0, 6}]]
  1. Y.M. Ali and L.C. Zhang, Relativistic heat conduction, International Journal of Heat and Mass transfer, 48, 2005, 2397--2406.
  2. Brovman, M.Ya., Application of isothermal surfaces for calculation of unsteady heat transfer processes, Journal of Engineering Physics and Thermophysics, Vol 71, No. 5, 1998.
  3. Brody, Jed and Brown, Max, Transient heat conduction in a heat fin, American Journal of Physics, 85, Issue 8, 582 (2017);
  4. Chester, M., Second sound in solids, Physical Review, 1963, 131, No. 15, pp. 2013--2015.
  5. J. H. Lienhard V and J. H. Lienhard, A Heat Transfer textbook
  6. Swenson, Robert J., Generalized heat conduction equation, American Journal of Physics, 46, Issue , 76 (1978);
  7. V. Vuorinen, Newton’s Cooling Law, Fourier’s Law and Heat Conduction, Aalto University, 2016.
  8. Vernotte, P., Les paradoxes de la theorie continue de l'equation de la chaleur, Comptes rendus de l'Académie des Sciences 246 (22) (1958) 3154--3155.


Return to Mathematica page
Return to the main page (APMA0340)
Return to the Part 1 Matrix Algebra
Return to the Part 2 Linear Systems of Ordinary Differential Equations
Return to the Part 3 Non-linear Systems of Ordinary Differential Equations
Return to the Part 4 Numerical Methods
Return to the Part 5 Fourier Series
Return to the Part 6 Partial Differential Equations
Return to the Part 7 Special Functions