Heat Transfer
| Contents | Appendix |
| Introduction | Nomenclature |
| Heat Equation | References |
| Boundary Conditions in Heat Transfer |
Contents
Introduction
This tutorial gives an introduction to modeling heat transfer. Governing equations and boundary conditions that are relevant for performing heat transfer analysis are derived and explained.
Heat transfer is a discipline of thermal engineering that is concerned with the movement of energy. The driving force for heat transfer is temperature differences. The temperature differences come about though different phenomena in the interior or on the boundary of the simulation domain and can be categorized into thermal conduction, thermal convection and thermal radiation. Combining all effects, the changes in a temperature field in a given region over time are then modeled with a heat equation.
The modeling process results in a partial differential equation (PDE) that can be solved with NDSolve. Furthermore, in this tutorial, different types of heat sources are introduced along with an overview of how various real-world thermal interactions can be modeled with the available thermal boundary conditions.
The accuracy and the effectiveness of the heat transfer PDE model are validated in the separate notebook entitled Heat Transfer Verification Tests.
Extended examples of heat transfer modeling can be found in the "Model Collection".
Many of the animations of the simulation results shown in this notebook are generated with a call to Rasterize. This is to reduce the disk space this notebook requires. The downside is that the visual quality of the animations will not be as crisp as without it.
The symbols and corresponding units used throughout this tutorial are summarized in the Nomenclature section.
This notebook makes use of finite element method functionality for various aspects during the solution of the PDE models.
Heat Equation
Introduction to Heat Equation
The heat equation (1), which is derived from the law of energy conservation, is used for modeling time-dependent heat flow within a thermally conductive medium:
The dependent variable in the heat equation is the temperature 
, which varies with time 
and position 
. The partial differential equation (PDE) model describes how thermal energy is transported over time in a medium with density 
and specific heat capacity 
. The specific heat capacity is a material property that specifies the amount of heat energy that is needed to raise the temperature of a substance with unit mass by one degree Kelvin.
Besides the time derivative part, the PDE is made up of several components. First and foremost, there is a diffusive term: 
with a thermal conductivity 
. The thermal conductivity or other quantities may very well depend on the temperature 
. This will then result in a nonlinear heat equation.
The second part is a convective term: 
with a flow velocity 
for modeling internal heat convection. This term is only present if the medium allows for an internal flow. If the simulation medium is a solid, then this term is zero.
The term 
denotes a heat source within the domain and is explained in the Source Types section.
Heat Equation Derivation
To derive the heat equation, start with energy conservation. Consider balancing the energy generated within a unit volume domain 
with the energy flowing through the boundary of the domain.
In the above graphics, 
is the mass density in units of 
and 
is the internal energy per unit mass in units of 
. The total energy within the control volume is then equal to the product 
. The red circle in the middle represents a heat source 
in units of 
. Since power 
is equivalent to energy per unit time 
, 
denotes thermal energy generated per unit time inside the domain 
. The heat flux 
in units of 
represents the net rate of energy that exits through the boundaries, such that 
.
The energy balance within the domain 
can then be described by the following equation:
That is, the time rate of change in the total energy 
is equal to the energy 
generated per unit time inside the domain minus the net rate of energy that exits the domain 
.
Here the energy flux 
can be divided into two parts: a convection term 
and a diffusion term 
. The convection term denotes the energy transported by a possible flow inside the medium and is in proportion to a flow velocity 
:
If the heat transfer occurs in a solid medium, then, because a solid can not have an internal velocity field 
by definition, the convection term is set to 
.
The diffusion term represents the energy flux resulting from the energy gradient and is proportional to its energy diffusivity 
:
For heat transfer modeling, the diffusion term (2) is often expressed in the form of the temperature gradient 
instead, which is known as Fourier's law of heat conduction:
Here, the energy diffusivity 
is represented by the thermal conductivity 
in units of 
, and the minus sign indicates that the heat diffusion is in the direction of decreasing temperature.
Note that the diffusive component is always present, regardless of the type of the medium.
The principal property of the diffusion term is smoothing; this is explained in the section The Smoothing Characteristic of the Diffusion Equation.
Inserting (3) and (4) into the energy balance equation (5) yields:
Note that the above equations hold for both continuous and discontinuous density/velocity fields. This is explained in more detail in the appendix Conservation Laws with Discontinuous PDE Coefficients.
Since the domain 
with a unit volume 
, the total mass within the domain is then equal to 
. Therefore, the term 
on the left-hand side can be interpreted as a mass conservation equation. Specifically, the time derivative term can be understood as the accumulation (or loss) of mass in the domain 
per unit time, while the divergence term part denotes the difference in the rate of the mass outflow and inflow.
The rate of the mass inflow is equal to the rate of the mass outflow plus the accumulation/loss of mass within the domain 
. When there is no internal mass generation/elimination, the term 
sums to zero and can be removed from the equation (6), leading to:
Note carefully that this leads to a restriction on the usage of the heat transfer model: Equation (7) and the resulting heat transfer PDEs (8) and (9) cannot be used when the mass in the medium changes.
Since the internal energy 
depends on the temperature: 
, equation (10) can be rewritten into the following form using the chain rule:
Here, the term 
is also known as the specific heat capacity 
, which denotes the ratio of the energy added to/removed from the domain to the resulting temperature change. With this definition, equation (11) simplifies to the heat equation:
The general heat equation describes the energy conservation within the domain 
and can be used to solve for the temperature field 
in a heat transfer model. Since it involves both a convective term and a diffusive term, the equation (12) is also called the convection-diffusion equation.
Note that the preceding heat equation (13) is written in a nonconservative form. That is, the mass density 
, the specific heat capacity 
and the velocity 
terms of the 
term are outside of the gradient operator. One could imagine a convective term inside the divergence operator 
. However, since both 
and 
might be space dependent, the general heat equation (14) cannot simply be transformed into a conservative form.
In a solid medium, however, the internal velocity field 
is set to zero, and the governing PDE simplifies to a pure conductive heat equation:
The heat equation may also be expressed in cylindrical and spherical coordinates. Please refer to the appendix section Special Cases of the Heat Equation for a detailed explanation.
Heat Transfer Model Setup
Defining a function that represents the spatial terms of a heat transfer model in Cartesian coordinates will make the setup of a heat equation more convenient.
The inputs needed for a heat transfer model are:
 | the temperature variable | |
 | the spatial independent variables | |
 | a velocity field | |
 | the thermal conductivity | |
 | the density | |
 | the heat capacity | |
 | a heat source/sink |
Note that this model definition uses inactive PDE operators. "Numerical Solution of Partial Differential Equations" has several sections that explain the use of inactive operators.
Model Parameter Setup
The following model parameters are used for the examples in this notebook. These parameters define the simulation domain 
, the simulation end time 
and thermal properties of a medium.
In some examples, a smoothed step function 
is used to prescribe a time profile for a transient parameter, for example, the heat flux 
or the surface temperature 
. The smoothed step function is defined as follows:
Here the minimum value and the maximum value the function 
can reach are denoted by 
and 
. The location of the step is controlled by 
, and the smoothed steepness is controlled by 
.
Basic Heat Transfer Example
The following 2D stationary example [15] demonstrates a typical workflow of heat transfer modeling.
The model domain 
of width 
and height of 
is a ceramic strip that is embedded in a high-thermal-conductive material. The side boundaries of the strip are maintained at a constant temperature 
. The top surface of the strip is losing heat via both thermal convection and thermal radiation to the ambient environment at 
. The bottom boundary, however, is assumed to be thermally insulated.
The goal is to find the steady-state temperature distribution of the ceramic strip.
The thermal conductivity 
, heat transfer coefficient 
, density 
, heat capacity 
and emissivity 
of the ceramic strip are given by:
A default thermally insulated boundary condition is implicitly applied on the remaining bottom boundary.
Note that the fluxes are part of the equation. The reason for this is explained in the section "Partial Differential Equations and Boundary Conditions".
In the steady state, the minimum temperature is found on the top surface that is cooled by both radiation and convection modeled with the heat transfer value. The defined temperature at the side is the maximum temperature, since heat diffuses into the medium from the sides.
The setting of thermal boundary conditions will be explained in detail in a following section: Boundary Conditions in Heat Transfer.
The Shrink Fitting of Assembly is another time-independent application example. A time-dependent application example can be found in the Laser Beam Welding application example.
Source Types
The source term 
in the heat equation (16) is used to model internal heat generation (
) or absorption (
) within the domain. Based on their shape, heat sources are categorized as a Volumetric Heat Source, Layer Heat Source and Point Heat Source.
It is important that the mesh conforms to the geometrical bounds of the source term 
. The best way to do this is by explicitly generating the mesh for them. The finite element method, which is used as the solution method, benefits from a mesh where the individual elements do not cross a material boundary. An alternative is to make use of the MeshRefinementFunction. In this case, elements will cross the material boundary but they will be small, if the mesh is sufficiently well refined. Combining both approaches is also possible.
Heat sources can be functions of time 
, space 
and other dependent variables, such as the temperature 
itself.
Volumetric Heat Source
A volumetric heat source 
can be used to model an arbitrarily shaped heat source (
) or heat sink (
) within the domain. A volumetric heat source 
has units of 
. The corresponding source strength 
denotes the rate of an internal heating or cooling per unit volume where the heat source is active.
The volumetric heat source value is always specified in units of 
. The name comes from the 3D incarnation of the heat source but is used in other dimensions as well.
Consider the 1D case, where 
. When you specify the volumetric heat source 
in units of 
, that value will be multiplied by the cross-sectional area 
in units of 
and results in a unit of 
.
Similarly, in a 2D domain, where 
, the volumetric heat source 
in units of 
is multiplied by the thickness 
in units of 
and results in units of 
.
In the following 2D example, a rectangular heat source 
is introduced to heat up the domain. The source strength is fixed at 
.
For more involved shapes of a heat source, the function RegionMember can be used to specify the source region 
.
In simple cases like this, evaluating the RegionMemberFunction at the symbolic spatial coordinates 
will lead to a simple Boolean expression that can directly be used in the heat source specification and will lead to an efficient time integration.
In some geometrically more complicated cases, evaluating a RegionMemberFunction at the symbolic spatial coordinates will return unevaluated.
While using a RegionMemberFunction directly is possible, it will be less efficient to time integrate. The reason the Boolean expression is more efficient is that it can be automatically compiled, while an unevaluated RegionMemberFunction cannot be compiled. Also, see this note about the setup of efficient PDE coefficients.
The most accurate and efficient method to deal with heat sources is by element markers, as then the mesh will have a specific subregion for the heat source, which will result in an accurate solution. More information on markers and their generation in meshes can be found in the "Element Mesh Generation" tutorial.
An example that uses element markers for the heat source setup is presented in the appendix section Modeling Heat Source by Element Markers.
See this note about improving the visual quality of the animation.
The simulation begins with an undisturbed domain where 
. With a volumetric heat source 
placed in the domain, thermal energy is generated and gradually heats up the domain. The speed of the heat transfer depends on the heat conductivity 
and the heat capacity 
of the material.
Point Heat Source
A point heat source 
can be used to model an internal heat source (
) or heat sink (
) that is considered to have no spatial extension. A point heat source can be made use of in 3D and in 2D axisymmetric domains.
In a 3D model, a point source can be positioned at any place.
In 2D axisymmetric models, a point source can only exist if it is specified on the axis of symmetry. A point off the axis of symmetry would imply a line source once the rotational aspect of the axisymmetric model is taken into consideration.
Point sources are not used in 2D models, as that would imply an equivalent to an out-of-plane line source. See the section Layer Heat Source for details on this.
A point source, shown in red, specified on the axis of symmetry of the 2D axisymmetric domain, in gray. One can see that revolving the axisymmetric domain does not transform the point source into a line source.
The units of the point source 
are always in units of 
, while the volumetric heat source is in units of 
. Generally speaking,
where 
is the total heat generation in units of 
.
where the 
is a Dirac delta function. The 
have a unit of 
in each of the directions of the source location 
This leads to an expression for the volumetric source 
that can be used:
The Dirac delta function, however, poses a problem in numerical simulations, as it cannot be resolved in the discretized spatial domain. This is because the Dirac delta function is singular at the source location 
. A second problem is that in the finite element method, the evaluation of coefficients always happens within mesh elements, never at the edges. Hence, an approximation to the Dirac delta function is needed. The process of approximating the Dirac delta function is called regularization.
There are various regularized delta functions 
available [17,18] such as:
where 
is the regularization parameter that controls the support of the regularized delta functions 
. Typically, 
should have a size comparable to the mesh spacing 
. 
represents the difference 
.
The concept for axisymmetric models is the same as before, and the expression for the volumetric heat source is given as:
Here the Dirac delta functions provide the units of 
at the source location 
in the axis of symmetry, where 
is on the 
plane and can be any number from 0 to 2
. In an axisymmetric case, 
does not need to be specified.
Layer Heat Source
A layer heat source models a heat source (
) or heat sink (
) that is too thin to have a thickness in the model geometry. A layer heat source can be made use of in 3D, 2D and 2D axisymmetric domains. In 1D domains, a layer heat source is the same dimension as the simulation domain, which would make it a volumetric source.
The units of 
are always specified in units of 
.
In 3D domains, layer sources can be specified along edges of a geometry, including lines floating arbitrarily in the geometry.
In 2D axisymmetric domains, a layer source is specified as a point that gets its length through the rotation of the simulation domain around the axis of symmetry.
In 2D, a point also represents a layer source because of the thickness 
of the model in the out-of-plane direction:
On the left: a 3D layer source and a plane. To the right: The same layer source from a 2D perspective.
In 3D models, a line charge can be specified at any edge of a geometry:
and the variations 
and 
are possible.
Here 
is a Dirac delta function at the source location 
and provides the units of 
, and 
is the value of the layer source.
To model a layer source out-of-plane, it is necessary to specify a point as a source.
Since the point has no spatial extension in all directions, the Dirac delta function 
should be applied on each dimension (i.e. 
) of the modeling domain 
:
Here 
is a Dirac delta function at the source location 
in units of 
, 
is the value of the layer source, and 
is the thickness in units of 
.
In the following 2D example, a layer heat source 
is added at 
to heat up the domain. The layer source has a value of 
.
To utilize the regularized delta function 
, the regularization parameter is chosen such that it is half of the mesh spacing: 
.
See this note about improving the visual quality of the animation.
The simulation begins with an undisturbed domain where 
. With a layer heat source 
placed in the domain, thermal energy is generated and spreads out in all directions. The speed of the heat transfer depends on the heat conductivity 
and the heat capacity 
of the material.
Anisotropic and Orthotropic Heat Transfer
In the previous sections, it was assumed that heated medium is isotropic; that is, the rate of heat transfer is independent of its direction given the same the temperature gradient 
. In reality, however, a medium may be anisotropic. This means that heat diffuses in different directions at a different rate. The diffusion term (19) is then rewritten as:
where 
is the thermal conductivity tensor. 
and 
are called the principal conductivity coefficients and off-diagonal conductivity coefficients, respectively.
Based on Onsager's [20] principle of the thermodynamics of irreversible processes, the off-diagonal conductivity must obey the reciprocity relation [21]:
Orthotropic heat transfer is a special case of anisotropic heat transfer. Here, the thermal conductivity of a material is symmetric along the principal directions 
, 
and 
. This means values along the principal direction are nonzero but unequal to each other. The off-diagonal conductivity coefficients are zero. This behavior can be seen, for example, in fiber composite materials. Then the thermal conductivity tensor becomes:
As an example, consider a 2D composite material with layered-like structure:
This case is an orthotropic case where the heat transfer is more efficient horizontally. The thermal conductivity tensor can then be described by:
See this note about improving the visual quality of the animation.
Unlike the example shown in the previous section, in this case the heat transfer is faster in the horizontal direction, resulting in a higher temperature zone within 
.
Nonlinear Heat Transfer
In the previous sections, it was assumed that the PDE coefficients, namely the density 
, the heat capacity 
and the thermal conductivity 
, are independent of the temperature field 
. In reality, however, these parameters might change significantly with temperature, especially for pure metals.
As an example, consider a 1D heat transfer model with an initial temperature field at 
and a temperature-dependent thermal conductivity 
:
Equation (22) is a nonlinear heat transfer model, since the conductivity coefficient 
in the PDE model now depends on the temperature 
itself.
To heat up the domain, a constant heat flux 
is applied on the left-hand boundary.
and an initial temperature field 
:To understand the effects of the nonlinearity compare to a linear heat transfer PDE.
The simulation begins with an undisturbed domain where 
. With a constant heat flux 
applied on the left side, thermal energy is then transferred across the boundary and heats up the domain.
For the nonlinear model, as the temperature 
increases, the thermal conductivity 
will increase correspondingly, which further speeds up the heat transfer and results in a flatter temperature field.
Temperature-Dependent Heat Capacity
What follows is a second nonlinear example.
In general, the properties of a heat-conducting material, such as the density 
, the heat capacity 
and the thermal conductivity 
, can depend on the temperature 
. Since temperature is also a dependent variable, this additional temperature dependency makes the PDE model nonlinear. In the following section, a model is presented where the heat capacity 
is temperature dependent such that 
, and thus is a nonlinear model.
To show the effect of the nonlinear specific heat capacity, the result of the nonlinear specific heat capacity model will be compared with that of a linear specific heat capacity model.
Consider a 2D heat transfer model with an initial temperature field at 20[°C] and a temperature-dependent heat capacity 
:
The specific heat capacity 
is a nonlinear function due to its dependency on the variable 
. Specific heat capacity has units of [
], while temperature has units of [
]. This means that there is an implicit factor of 1 [
], which multiplies the temperature function to match the units of the specific heat capacity.
Next, the simulation needs to run sufficiently long such that the system comes to thermal equilibrium with the ambient. To detect the time when the system is in thermal equilibrium, a WhenEvent is used. The use of a WhenEvent has the distinct advantage that NDSolve has special mechanisms built in to detect a specified event during the time integration. More information on events can be found in the section Heat Transfer with Events.
For this specific simulation, the event for stopping the time integration is specified by choosing a tolerance value for the difference in the current temperatures at a specific point and the ambient temperature:
:NDSolve still needs an end time for the simulation. An arbitrary value that is larger than the time needed to reach equilibrium is chosen.
Next, to better understand the effects of the nonlinearity in the specific heat, a comparison is made to a linear heat transfer PDE where the specific heat is set to the constant ambient temperature used as the initial condition.
The linear model reached the equilibrium point earlier than the nonlinear model.
See this note about improving the visual quality of the animation.
The simulation begins with an undisturbed domain where 
=20[°C]. With convective heat transfer taking place from the ambient to all sides except to the bottom boundary, which is insulated, thermal energy is transferred across the boundaries and heats up the domain. As time progresses, the temperature of the domain equilibrates with the ambient.
For the nonlinear model, as the temperature 
increases, the heat capacity 
will increase correspondingly, which slows down the heat transfer, resulting in a relatively lower temperature profile compared to that of a linear model. This leads to a slightly faster equilibration of the system with the ambient for the linear model case. However, the gradients in the temperature field inside the domain ultimately subside as time progresses, and even the nonlinear model system equilibrates with the ambient. At steady state, no differences in the linear and nonlinear model system temperature profiles are seen.
Heat Transfer with Events
A common topic in heat transfer modeling is to simulate dynamic or pulsed thermal loads, that is, a heat flux that turns on and off under different conditions. In such cases, a WhenEvent can be used to construct and efficiently solve the heat transfer model.
Consider a 1D room model where heat continuously flows out of the room at both sides. A heater is placed in the domain to warm up the room, but will only turn on when the temperature at the center 
drops below the threshold temperature 
. To prevent overheating, the heater will turn off when 
is above 
.
To model the conditional heating of the heater, a volumetric heat source 
is applied using WhenEvent. The heat source 
is switched on and off when the central temperature 
reaches 
or 
.
To model the cooling of the room, a constant cooling flux 
is applied on both ends of the domain.
Next, the heat transfer PDE is defined with the conditional heat source 
. Since the heater is off at the beginning, the initial value of the heat source is set at 
.
To solve this heat transfer PDE, 
is specified as a discrete variable, which means it only changes at discrete times during the temporal integration by NDSolve.
The simulation begins with a uniform temperature at 
. With a constant cooling flux 
applied on both sides, heat continuously flows out of the domain and brings down the room temperature. The heat source 
is then turned on and off when the central temperature 
reaches 
or 
, resulting in an oscillating temperature field.
While it would also be possible to make use of an If statement to model the dynamic or pulsed heat source, the use of WhenEvent has the distinct advantage that NDSolve has special mechanisms built in to detect the events during the time integration. This mechanism may not be available when modeling the pulsed heat source with an If or similar statement.
Details about modeling heat pulses are presented in the appendix section Possible Issues and Workarounds for Modeling Heat Pulses.
Heat Transfer in Porous Media
Porous media are multiphase objects with a solid skeleton portion and a porous region that is filled with a fluid. Due to their special thermal and mechanical properties, porous materials have been widely used in many industrial applications, such as vibration suppression, heat insulation and sound absorption.
To model heat transfer within a porous medium, one approach, called a direct approach, is to build a coupled PDE with two heat equations with material coefficients suitable for each phase. One equation describes the temperature field 
in the solid region, and the other equation models the temperature 
in the fluid region:
Here, the subscript 
denotes the parameters of the solid phase and 
parameters of the fluid phase. The two equations are coupled by the volume fraction 
of each phase, and the heat exchange between two phases is accounted for explicitly by an additional heat source/sink term 
on the right-hand side.
However, to make use of this approach, it is required to reproduce the entire pore structure of the domain. Due to the geometric complexity of the microscopic porous structure, a fine, finite element mesh may be required to resolve the geometry accurately and thus a significant computational cost to solve the model can be expected.
An alternate approach is to model pores on a macroscopic scale. In this case, the heat transfer model uses an average temperature field to describe the entire porous structure of both phases. To do so, volume-averaged effective thermal properties are used, and the model can be expressed by a single heat equation:
Here 
is the effective volumetric heat capacity and 
is the effective heat conductivity, which are computed based on the volume fraction 
and the properties of each phase:
Assuming the porous medium to be fully saturated, the volume fraction of both phases can be related by 
. Then the equation (23) becomes:
Note that the fluid's volume fraction 
is also known as the porosity of the medium. For the sake of simplicity, its subscript 
will be omitted in the following section.
As an example, consider a 2D heat transfer model of a porous medium. A fluid flow is passing through the domain 
with a width of 
and a height of 
, and a constant heat flux is applied on the left surface to heat up the domain. To understand the effects of the porosity on the heat transfer, the model will be solved with three different porosity values of 
, 
and 
(purely fluidic medium).
Equation (24) has an effective volumetric heat capacity in front of the time derivative and at the same time uses material parameters for the fluidic phase. To accommodate for that, the equation is generated in part with HeatTransferPDEComponent and an additional time derivative term. For that, the model variables include time as a dependent variable but not the time variable itself. The time derivative with its effective volumetric heat capacity will be added manually.
To heat up the domain, a constant heat flux 
into the domain is applied on the left surface at 
.
and a right-going flow 
:
Note that the temperature field is symmetric in the 
direction.
To study the effects of the porosity on the heat transfer, the temperature evolution along the 
axis can be compared at three different porosity values of 
, 
and 
(purely fluidic medium).
Since the heat capacity value in the solid phase is smaller than in the fluid phase: 
, the medium with a larger amount of solid (i.e. lower porosity 
) will be more susceptible to the temperature change. This can be verified by computing the effective volumetric heat capacity 
at 
, 
and 
.
Note that the effective heat capacity 
is the smallest for 
:
Heat Transfer Model with Mixed Dimensions
The following section demonstrates how to model heat transfer phenomena defined in different spatial dimensions. As an example, consider a 2D ceramic strip with a uniform initial temperature of 
. At the left surface, the ceramic is cooled down by attaching it to a thin pad with a cooling flux through it. The cooling flux 
is proportional to the temperature difference between the pad and the ceramic, with a heat transfer coefficient 
.
A constant heat flux 
is applied on the top of the ceramic. The right and bottom surfaces of the ceramic are assumed to be thermally insulated.
To solve for the temperature field 
of the ceramic, a first idea might be to build a single 2D system and model the cooling pad as a heat flux boundary condition. In this example, however, the cooling flux 
depends not only on 
but also on the pad temperature 
. In order to determine the value of the pad temperature 
, it is necessary to also model the heat transfer in the thin cooling pad.
Assuming the cooling pad to be much thinner than the width of the ceramic, the temperature variation of the pad in the 
direction can be neglected. The cooling pad can thus be modeled as a 1D region, while the ceramic strip is modeled as a 2D region. In other words, a mixed-dimensional model can be used.
For simplicity, the heat transfer coefficient 
and the thermal conductivity 
, density 
and heat capacity 
for both the ceramic strip and the pad are set to one.
The temperature field 
within the ceramic strip is described by a 2D heat equation:
The pad is described by the 1D domain along the 
axis 
, which coincides with the left surface of the ceramic strip. A 1D heat equation can be used to model the pad temperature 
along the 1D pad region as:
However, to solve a coupled PDE system, NDSolve requires all the dependent variables to have the same spatial dimensionality. For this reason, a "fictitious dimension" has to be introduced for the pad temperature 
in the 
direction, yielding:
That is, the pad temperature 
will be solved not only in the 1D pad region, but also the entire 2D domain of the ceramic strip. However, since 
is just a fictitious dimension for the pad, there is no physical meaning to a 
value within the ceramic region 
, which means the resulting pad temperature field 
will only be valid along the 1D pad region on the 
axis 
.
Next, the heat exchange between the cooling pad and the ceramic strip needs to be considered. From the perspective of the ceramic, heat is lost through the left ceramic boundary to the pad and can be modeled with a heat flux boundary condition.
From the perspective of the cooling pad, heat is gained from the ceramic strip over the entire pad domain. This can be modeled by a heat source term 
in the equation (25).
Due to the law of energy balance, the heat source 
of the pad should have the same magnitude but an opposite sign from the cooling flux 
of the ceramic.
A constant heating flux 
is applied on the upper surface of the ceramic.
See this note about improving the visual quality of the animation.
Next, the pad temperature 
within the 1D pad region 
is inspected and compared to the ceramic temperature 
on the left ceramic surface. A custom function, TwoAxisPlot, is defined and applied to rescale the temperature range in the plot for better visualization.
Since the heat gain of the pad depends on the ceramic temperature 
, the pad temperature 
follows a similar pattern with 
along the left ceramic surface 
.
Introducing a fictitious dimension in the coupled PDE system enables solving a mixed-dimensional model involving a 1D and a 2D heat equation. This technique can also be applied in dimensions 1, 2 or 3D in a similar manner.
Heat Transfer in Multi-material Media
Heat transfer in multi-material domains is a fundamental problem of interest to many industrial applications. The varying material properties of the composite along the spatial coordinates of the geometry under consideration leads to a more complicated heat transfer mechanism than seen in an isotropic uni-material medium. The following example illustrates the setup and solution of heat transfer through a composite material with suitable boundary conditions.
Consider a 1D transient heat transfer multi-material model with an initial temperature field at 0 [
]:
The 1D composite is made up of three different materials of thicknesses 0.25, 0.114 and 0.04 [
] with varying properties, such as thermal conductivity of 8, 1.8 and 44 [
], density of 3100, 2100 and 7800 [
] and heat capacity of 1050, 1100 and 540 [
]. The left end of this multi-material medium is subjected to a high temperature of 1700 [
], and at the right end, heat is getting dissipated by the combined effect of convection and radiation. The initial temperature of the multi-material domain is maintained at 297 [
]. This problem is solved by specifying the region-dependent property values in the governing heat transfer model, along with the specified boundary and initial conditions.
To specify the properties of multiple materials in a given region, region markers are used, and they are specified as a part of the ToElementMesh. More information about specifying region markers in the mesh can be found in the section on markers in the "Element Mesh Generation" monograph.
As a next step, Dirichlet boundary conditions are set up at the left boundary and convective and radiative boundary conditions on the right boundary of the multi-material domain. The left boundary is maintained at a high temperature of 
, while heat gets dissipated at the right boundary by the combined effect of convection and radiation to the ambient.
To analyze the heat conduction between the three layers, the model is solved from 
to 
.
The temperature of the system is initially maintained at 
throughout. When the temperature of the left end is switched to 
suddenly and heat is dissipated at the other end by convective and radiative flux that depends on the time-dependent temperature at that boundary, the temperature in each division of the multi-material domain decreases with a discontinuous gradient at the junctions that depends on the material properties of the adjacent domains. As time progresses, the temperature profile in each division attains steady state with a constant gradient.
The Contactless Anemometer or the Heat Conduction in a Multilayer Sphere application models show further cases of multi-material heat transfer.
Heat Transfer with Phase Change
A phase change in thermodynamics denotes the phenomenon where materials transit from one state (solid, liquid, gas, plasma) to another, which only occurs at certain temperature and pressure and when sufficient energy is added or removed from the system. The energy associated with the phase change, which is known as the latent heat 
, is used to alter the molecular structure instead of creating a temperature change of the material. A phase change is sometimes also called a phase transition. The mathematical term Stefan problem is also common.
Ice-to-Water Solidification
As an example, consider the following 1D model that describes the ice-to-water phase change along a 
bar of ice. The rod has an initial temperature of 
, and a constant heat flux 
is applied at the left end to melt the rod. At the right end, the rod is assumed to be thermally insulated.
In the heat transfer model, instead of simulating the phase transition exactly at the phase change temperature 
, assume that the transition occurs in a temperature interval: 
to 
. The material phase during the transition is then described by a smoothed step function 
, which denotes the ratio of the original phase to the new phase within the material.
Within the temperature interval: 
, the equivalent density 
and conductivity 
are given by:
To make things more general, the smoothed step function is also an argument to the equivalent density 
and conductivity 
.
The equivalent specific heat capacity 
, however, should include an extra term 
to account for the latent heat 
required for the phase transition. Here 
denotes the distribution of the latent heat 
during the phase change, and is approximated by a regularized delta function around the phase change temperature 
:
Note that the integral of 
equals the latent heat 
required for the phase change:
The equivalent heat capacity 
is then given by:
Note that 
depends on three parameters: 
, 
and 
, the distribution of the latent heat. 
is a function that was specified above. In a later experiment, 
will be replaced with a different function.
A default thermally insulated boundary condition is implicitly applied at the right end of the rod.
To better understand the effects of the latent heat 
on the phase change, the above result will be compared with a solution that neglects the latent heat. For this case, 
is set to a function that returns 0 for all input.
As mentioned above, a function for 
is now being set up that is always 0, in order to neglect the effect of latent heat. To do so, a function is created that returns 0 for any input.
To answer the question of how the choice of the temperature interval 
influences the simulation, another experiment is conducted. Here a temperature interval of 
is used.
The absolute temperature difference in the solutions between using 
and 
is about 
. Compared to the overall temperature of about 
, this is small.
Freezing of Liquid in a Pipe
This second example will model the freezing of liquid water in a pipe. The water in the pipe is initially at a temperature of 
. The right end of the pipe is maintained at the same temperature throughout. The left end of the pipe is cooled to 
for a certain amount of time. Then, the temperature goes up to 
again.
The basic concept, functions and material parameters are the same as above.
As a next step, Dirichlet boundary conditions are set up at the left and right boundaries of the domain. As indicated earlier, the right boundary is always maintained at a temperature of 5[°C], while the left boundary is subjected to a change in temperature.
In the time span from 
, a smoothed step function is used to change from 5[°C] to -5[°C]. Then, in the time span from 
, the temperature remains at 
. In the time span from 
, the left-hand side is heated to 
, at which temperature it remains for the remainder of the simulation.
The usage of a smoothed function for the transition from one temperature to another is to have a physical transition between them. A jump from one temperature to another without intermediate values is an unphysical behavior and should be avoided.
The temperature of the system is initially maintained at a temperature of 
throughout. When the temperature of the left end is switched to 
very quickly, the formation of ice at that end can be seen. This ice layer penetrates more into the system toward the right end until the time when the temperature of the left end is switched back to 
. At this point, the melting of ice from both the ends is seen until all the ice phase in the system disappears and gets converted to the water phase. Once the entire system is occupied by the water phase, the system gets heated up more rapidly as the temperature rise only needs the supply of sensible heat. Ultimately, the system reaches the steady-state temperature of 
.
Heat Transfer with Model Order Reduction
Sometimes one wants to rerun the same heat transfer simulation with different initial data. In this case, model order reduction can be of help. The idea behind model order reduction is to make use of the discretization of the PDE and find a similar discretization that can be time integrated much more efficiently. In order to perform a model order reduction, access to the discretization NDSolve makes is necessary. This is currently not possible on the NDSolve level but requires a bit of programming and is explained in the section "Model Order Reduction of Transient PDEs with Stationary Coefficients and Stationary Boundary Conditions" in the "Finite Element Programming" tutorial.
Multiphysics Heat Transfer
Heat transfer is often combined with other fields of physics. What follows is a list of multiphysics application examples that make use of heat transfer:
Boundary Conditions in Heat Transfer
The most common boundary conditions in heat transfer modeling can be modeled with DirichletCondition, NeumannValue and PeriodicBoundaryCondition and can be categorized in the following four types:
- Dirichlet-type boundary conditions. This type of boundary condition specifies the temperature
at the boundary and can be modeled with DirichletCondition.
- Neumann-type boundary conditions. This type of boundary condition specifies the heat flux
at the boundary and can be modeled with NeumannValue.
- Robin-type boundary conditions. This type of boundary condition specifies the relation between the temperature
and its normal derivatives at the boundary and can be modeled with a NeumannValue, since Robin-type boundary conditions are technically generalized Neumann boundary conditions.
- Periodic boundary conditions. This type of boundary condition specifies the temperature
at one part of the boundary to be the same at another part and can be modeled with PeriodicBoundaryCondition.
Under these four types, the following boundary conditions are introduced:
The following section describes several physical boundaries commonly encountered in heat transfer and how they can be modeled with the use of DirichletCondition, NeumannValue and PeriodicBoundaryCondition. For this purpose, the boundary condition currently being discussed is always on the left-hand side of the simulation domain. In some examples, additional boundary conditions are specified on the right-hand side to better demonstrate the behavior of the boundary condition on the left-hand side.
Surface Temperature Boundary Condition
Purpose
The purpose of a surface temperature boundary condition is to set a specific temperature on some part of the boundary.
Formulation
With a specified temperature 
on the boundary 
, the temperature surface condition is given by:
Derivation
A temperature boundary condition exists when the surface temperature 
is prescribed on a boundary. The surface temperature 
can be either a constant or time-dependent value and is set with a DirichletCondition in the heat transfer PDE model.
To model, for example, a heating wall that sends thermal energy into the domain, a transient surface temperature 
can be set up at the left end. Note that a Neumann zero condition is implicitly applied at the right end as a thermal insulated boundary.
Here a smoothed step function is used to describe the profile of the surface temperature 
from 
to 
. The parameters 
and 
are arbitrarily chosen to simulate the heating process.
The simulation begins with an undisturbed domain where 
. As the surface temperature 
increases at the left boundary, the excess thermal energy is then transferred to the right and brings up the temperature throughout the domain. The speed of the heat transfer depends on the heat conductivity 
and the heat capacity 
of the material.
Heat Flux Boundary Condition
Purpose
The purpose of a heat flux boundary condition is to model the rate of thermal energy flowing into or out of some part of the boundary per unit area.
Formulation
With a prescribed heat flux 
on the boundary 
, the heat flux boundary condition is given by:
Derivation
A boundary where the heat flux 
normal to the boundary is specified and not equal to zero is called a heat flux boundary:
By convention, a negative sign is added in front of 
to indicate that the heat flux is specified opposite to the outward normal 
. Therefore, a positive value of 
denotes the inward heat flux where the thermal energy enters the domain, and a negative 
denotes an outward flux.
Fourier's law of thermal conduction (26) relates the heat flux 
with the temperature gradient 
:
Inserting (27) into (28), the heat flux boundary condition can be written as:
Note that the unit of heat flux 
depends on the dimension of the boundary. In 1D (
), 2D (
) and 3D domain (
), 
has a unit of 
, 
and 
, respectively.
In the following example, a transient heat flux 
is applied on the left boundary when 
to heat up the domain.
The profile of the heat flux is defined as:
With the heat flux 
applied on the left boundary, thermal energy flows across the boundary and gradually heats up the domain. The heat flux is turned off at time 
. At this point, the character of the boundary condition changes from a heat flux to an insulation. The uneven temperature field is then smoothed out over time by the internal heat transfer.
When looking at the animation, it appears that there is a jump in the solution at 
. This is an artifact of the number of frames computed.
Note that the value of the heat flux is related to the temperature gradient 
by Fourier's law: 
. That is, the heat flux 
directly controls the temperature gradient normal to the boundary.
Thermally Insulated Boundary Condition
Purpose
The purpose of a thermally insulated boundary condition is to model a boundary where there is no heat flux across it.
Formulation
A thermally insulated boundary condition is given by:
If on some part of the boundary no boundary condition is set, an implicit Neumann zero boundary condition is used.
Derivation
A thermally insulated condition denotes a boundary where there is no heat flux across it:
Inserting (29) into the heat flux boundary condition (30), then the thermally insulated boundary condition can be written as:
In the following example, an insulated boundary is placed on the left-hand boundary, and a constant heat flux 
is added on the right end to serve as a heat source.
If no boundary condition is specified on any part of the boundary, then by default a Neumann zero boundary condition is implicitly used. This implies that the thermally insulated boundary is the default boundary condition used if no boundary condition is specified at a given boundary.
With a constant heat flux enforced on the right boundary, the temperature gradually increases within the domain. On the insulated boundary at the left end, however, the temperature gradient 
remains at zero at all times.
Since the temperature gradient 
is related to the heat flux 
by Fourier's law: 
, a zero temperature gradient implies a zero heat flux on the boundary (i.e. thermal insulated boundary).
Symmetry Boundary Condition
Purpose
A symmetry boundary condition is used when the computational domain and the expected temperature field have mirror symmetry along an axis of the simulation domain.
Formulation
The symmetry boundary condition is given by:
If on some part of the boundary no boundary condition is set, an implicit Neumann zero boundary condition is used.
Derivation
A symmetry boundary condition is used to reduce the extent of the computational domain to a symmetric subdomain of the full physical model geometry. This allows for a faster solution process with a lower memory requirement.
Consider the case of solving the temperature field of a 1D system from 
to 
. If the temperature pattern is expected to have a mirror symmetry along 
, the simulation domain 
can be efficiently constructed with only the left half of the system. Then a symmetry boundary condition should be applied at 
.
Due to the symmetry, the temperature gradient at the symmetry boundary will remain at zero at all times, which implies a zero heat flux across the boundary. Therefore, a symmetry boundary condition is equivalent to a thermally insulated boundary condition.
Outflow Boundary Condition
Purpose
If the heat transfer occurs in a fluid medium where the flow velocity 
, then an outflow boundary condition is used to model an outlet where heat is transferred out of the domain by fluid flow.
Formulation
When modeling heat transfer in a fluid medium, the outflow boundary condition at the outlet 
is given by:
If on some part of the boundary no boundary condition is set, an implicit Neumann zero boundary condition is used.
Derivation
When modeling heat transfer with a fluid flow, the diffusion heat flux 
is set to zero at the flow outlet boundary. This condition means that the temperature field of the flow outside the domain is assumed to have no impact on the flow inside the modeling domain 
.
The outflow boundary condition can only be applied on fully developed flows. That is, at the flow outlet, the velocity profile 
is unchanging in the flow direction.
In a case where there is recirculation through the outlet boundary, which often happens for turbulent flow, the reentering flow will affect the temperature field of the flow inside the domain and break the zero diffusion flux assumption. In this situation, the outflow boundary condition is no longer applicable.
Since the outflow boundary condition is essentially a Neumann zero condition, it will be implicitly applied if no boundary condition is specified at a given boundary.
Convective Boundary Condition
Purpose
The purpose of a convective boundary condition is to model thermal energy transferred across a boundary induced by a flow adjacent to the same part of the boundary.
Formulation
Given the profile of an external temperature 
and a heat transfer coefficient 
on the boundary 
, the convective boundary condition is given by:
Derivation
With the existence of a fluid flowing outside of the domain, adjacent to the boundary surface, part of the thermal energy will be transferred across the boundary through the movement of fluid particles, which is known as convective heating or cooling.
In 1701, Newton found that the rate of convective heat transfer between two media is proportional to their temperature difference. The convective heat flux is therefore defined as:
Here 
denotes the temperature of the external fluid, and 
in units of 
is the convective heat transfer coefficient. The heat transfer coefficient is determined experimentally and depends on material properties like density 
, the thermal diffusivity 
and the flow situation of the external fluid such as viscosity 
and Rayleigh number 
.
The approximate range of the convection heat transfer coefficients 
is presented in the following table:
Several empirically based formulas are also built [31, 32] to estimate the heat transfer coefficient 
in different situations.
Inserting (33) into the heat flux boundary condition (34), the convective boundary condition can be written as:
In the following example, a constant external flow at 
is applied on the left boundary to heat up the domain 
of an initial temperature field at 
. The convection heat transfer through the boundary is modeled by the convective boundary condition with a given heat transfer coefficient 
.
With the external flow flowing outside the domain, left of the left boundary, excess thermal energy is transferred across the boundary into the domain. Since the convective heat flux is proportional to the temperature difference across the boundary: 
, the temperature gradient 
at the left boundary will gradually decrease as the temperature 
at the boundary approaches the external temperature 
.
Thermal Radiation Boundary Condition
Purpose
The purpose of a thermal radiation boundary condition is to model heating or cooling through radiation on some part of the boundary.
Formulation
Given an ambient temperature 
, the surface emissivity 
and the Stefan–Boltzmann constant 
on the boundary 
, the thermal radiation boundary condition can be written as:
Derivation
All bodies with a temperature above the absolute zero (i.e. 
) will constantly emit thermal energy through electromagnetic radiation. The amount of radiation depends on both body temperature and surface condition. For a perfect thermal radiator, a black body, the Stefan–Boltzmann law states that the emitting heat flux is proportional to the fourth power of the body's absolute temperature:
where 
is the Stefan–Boltzmann constant.
In practice, however, the actually emitted heat flux is less than that of the black body radiation by a fraction 
known as the "surface emissivity". The value of the emissivity is 
and depends on factors such as physical properties and surface condition of the radiative body.
The Stefan–Boltzmann law can be rewritten as:
Based on (35), a radiation boundary condition can be formulated by inspecting the net radiative heat flux 
across a boundary:
The emitted outgoing radiative flux 
from a boundary depends on its surface temperature 
and the emissivity 
by:
At the same time, a boundary will absorb the radiation coming from the environment. This absorbed incoming radiative flux 
is given by:
Here 
is the ambient temperature and 
is the emissivity of the environment. Note that an extra term 
, the surface absorptivity factor, is multiplied on the right-hand side to account for the absorptivity of the boundary. This term denotes the ratio of the actual absorbed radiative flux to the total arriving flux.
Therefore, the net radiative heat flux across the boundary is given by:
To satisfy the thermodynamic equilibrium, for an arbitrary body the absorptivity 
should be equal to its emissivity 
. This is known as Kirchhoff's law of thermal radiation [36].
Equation (37) is then simplified to:
Assuming that the ambient surroundings behave as a black body with the emissivity 
, the equation can be further simplified as:
Inserting (38) into the heat flux boundary condition (39), then the radiative boundary condition is given by:
Note that the above derivation is performed based on the absolute temperature. That means the temperature terms in (40) have the unit of Kelvin 
.
To apply a radiation boundary condition in Celsius 
, an unit conversion should be done in (41):
Here 
denotes the temperature of absolute zero.
As an example, consider the ambient temperature of 
and a surface emissivity 
at the left-hand boundary. The net radiative heat flux across the boundary is modeled by the radiation boundary condition.
Due to the lower ambient temperature 
on the left end, the net radiative heat flux flows out of the system from the left-hand boundary, which gradually cools down the domain.
Similar to the convective boundary condition, the temperature gradient 
on the left boundary depends on the temperature difference across the boundary and can be calculated by Fourier's law: 
.
The temperature gradient on the left end is shown to be 
.
Periodic Boundary Condition
Purpose
The purpose of a periodic boundary condition is to map the temperature from one part of a boundary to another in order to model periodicity of the domain.
Formulation
Given a function 
that maps the temperature 
from the periodic boundary 
to the targeted boundary 
, the periodic boundary condition can be written as:
More information on how to specify PeriodicBoundaryCondition can be found on the reference page.
Derivation
A periodic boundary condition is applied to compute heat transfer in spatially periodic domains. Given a targeted boundary 
, the temperature 
on a periodic boundary can be mapped to the temperature 
on the targeted surface by a prescribed function 
. The boundary condition is set by the PeriodicBoundaryCondition in the heat transfer PDE model.
As an example, a ring heater is modeled, where a thin film is inserted to serve as a heat source. It is possible to perform the simulation with a 1D domain by using the periodic boundary condition.
The ring heater is converted into a 1D model with the length 
, which is the perimeter of the ring. To simulate the heating process, the temperature 
of the heating film is specified at the left boundary. At the right-hand boundary 
, a periodic boundary condition is applied to map the boundary temperature 
to the target boundary 
.
Here a smoothed step function is used to prescribe a temperature profile 
on the heating film. The parameters 
and 
are chosen arbitrarily to fit the parameters of the heat transfer model.el.
The simulation begins with an undisturbed domain where 
. During the heating process, the temperature 
on the periodic boundary has been mapped to the left-hand boundary, while the temperature 
of the heating film increases. The thermal energy is then transferred from both ends and gradually heats up the domain.
As an alternative, in this example, the PeriodicBoundaryCondition could be replaced by also using the DirichletCondition from the left boundary at the right boundary.
It is also worth noting that if the left-hand side had been any other boundary condition, the PeriodicBoundaryCondition would have projected that to the right. In other words, the PeriodicBoundaryCondition projects whatever boundary condition it finds where the mapping function 
points to, to the target that is specified as the predicate.
Appendix
Special Cases of the Heat Equation
Stationary Case
If the temperature field is at steady state, the transient term in (42) vanishes and the heat equation simplifies to:
Heat Equation in Cylindrical Coordinates
When modeling a heat transfer problem, sometimes it is not convenient to describe the model in Cartesian coordinates 
. In such cases, the heat equation may also be expressed using a cylindrical or spherical coordinate system.
A graphic showing cylindrical coordinates:
In the cylindrical coordinate system 
, 
and 
denote the radial, azimuthal and vertical directions, respectively. In terms of the Cartesian coordinates 
, the cylindrical coordinates are defined by:
By inserting the coordinate relations (43) into (44), the heat equation can be expressed in cylindrical coordinates as:
If the heat transfer within a model is rotationally symmetric about the 
axis, the resulting temperature field 
will be invariant in the 
direction. The equation (45) then simplifies to:
In that case, a 3D heat transfer problem can be modeled in a 2D domain by making use of this symmetric property. This type of model is known as an axisymmetric model.
HeatTransferPDEComponent can generate the axisymmetric form of the heat transfer equation (46) by specifying the parameter “RegionSymmetry” and setting it to “Axisymmetric”.
Examples that make use of the cylindrical coordinate system with a 2D axisymmetric model can be found in the separate Heat Transfer Verification Tests notebook: one time-independent 2D example and one time-dependent 2D example. Several other examples can be found on the HeatTransferPDEComponent reference page.
The Radial Effects in a Tubular Reactor application model also makes use of an axisymmetric model.
Heat Equation in Spherical Coordinates
A graphic showing spherical coordinates:
In the spherical coordinate system 
, 
and 
denote the radial, azimuthal and polar directions, respectively. In terms of the Cartesian coordinates: 
, the spherical coordinates are defined by:
By inserting the coordinate relations (47) into (48), the heat equation can be expressed in the spherical coordinate system as:
The Smoothing Characteristic of the Diffusion Equation
The basic behavior of a diffusion equation is smoothing. To see the effect, set up initial conditions that are discontinuous at 
.
The predicate True is used to specify the boundary condition on all boundaries. This may include internal boundaries, if present.
In a short time, the sharp tip of the discontinuous initial condition at 
has been smoothed out. This smoothing effect is a major characteristic of the diffusion equation and therefore also of the heat equation.
Modeling Heat Sources Using Element Markers
The most accurate and efficient method to deal with heat sources is by making use of element markers, as then the mesh will have a specific subregion for the heat source, which will result in an accurate solution. This is mostly because the finite element method can do a better job of solving the PDE when individual mesh elements do not cross material boundaries. Element markers and their use in meshes is explained in detail in the section "Element Marker" in the "Element Mesh Generation" tutorial.
In the following 2D example, a rectangular heat source 
is introduced to heat up the domain.
Next, region markers with the "RegionMarker" option for ToElementMesh are specified. To do so, a coordinate within the heat source should be given, as well as an integer marker. Optionally, an additional maximum cell measure can be specified to refine the source region.
An alternative to generate the mesh is to make use of the Boolean region functions and specify that no region holes should be inserted. The same region marker process from above is used.
See this note about improving the visual quality of the animation.
This result agrees with the one that did not make use of element makers, which is shown in the section Volumetric Heat Source. For complicated geometries, the use of region element markers will be easier to set up and compute the solution more efficiently.
Conservation Laws with Discontinuous PDE coefficients
To explain why the heat transfer PDE holds for discontinuous density/velocity fields, start with the mass conservation:
Assume that the density 
and the flow velocity 
are discontinuous at the interface ℐ:
Since there is no mass creation or destruction on the interface ℐ, the mass flux should still be the same on both sides of ℐ:
That is, the mass flux 
should be continuous throughout the domain. Then the divergence theorem [49] can be applied on the equation (50) leading to:
Since the domain 
is completely arbitrary, the integral can be discarded and the mass conservation equation in the differential form can be yielded:
That means the mass conservation equation (51) also holds true for the domain with discontinuous density/velocity fields.
The heat equation is essentially an energy conservation equation and can be derived in a similar way by substituting the density 
in (52) with the internal energy 
. Therefore, the heat transfer model presented in the tutorial can be applied on both continuous and discontinuous density/velocity fields.
Possible Issues and Workarounds for Modeling Heat Pulses
The section Heat Transfer with Events mentions that it is possible to use an If statement to model heat pulses with a time-dependent heat transfer model. However, when pulse durations become very small, the default time-stepping algorithm of NDSolve may fail to detect the prescribed heat pulse.
The following section demonstrates the issue with an example, followed by two workarounds for modeling heat pulses.
Consider a 1D time-dependent heat transfer model with periodic heat pulses applied in the middle part of the domain. During the time span 
, there is a total of five pulses with a duration of 
.
To model the cooling of the domain, the temperature at both ends are fixed at 
.
For demonstration purposes, the thermal conductivity 
, the density 
and the heat capacity 
are set to one.
Note that without using options, NDSolve missed the last two heat pulses at 
and 
.
Next, two possible workarounds to deal with this issue are presented.
Method 1—Reduce MaxStepFraction
The first method is to specify a smaller time step for the time integration process by NDSolve. The option "MaxStepFraction" allows the user to set the maximum fraction of the total time range to cover in a single step.
With a smaller time step, NDSolve successfully captured all five heat pulses between 
.
Method 2—Using WhenEvent
Another, better approach is to use WhenEvent to specify the heat pulses. The use of WhenEvent has the distinct advantage that NDSolve has special mechanisms built in to detect the pulses during the time integration. This mechanism may not be available when modeling heat pulses with an If or similar statement.
This case uses an event that is triggered when either of the expressions is zero. Each trigger will "toggle" the discrete variable 
.
More details about the usage of WhenEvent and its event detection method can be found here.
In this approach, the heat pulses are treated as a discrete variable 
by NDSolve, that is, it only changes at discrete times during the time integration. Note also that time integration step-size reduction is not necessary to obtain a solution.
With the usage of WhenEvent, NDSolve successfully captured all five heat pulses between 
.
Nomenclature
References
1. Bilbao, S. and Hamilton, B. Directional Source Modeling In Wave-Based Room Acoustics Simulation. IEEE, 2017.
2. Peskin, C. The Immersed Boundary Method. Cambridge University, 2002.
3. Churchill, S. W. and Chu, H. H.S. "Correlating Equations for Laminar and Turbulent Free Convection from a Vertical Plate." International Journal of Heat and Mass Transfer, 18 (11): 1323–1329, 1975.
4. Sukhatme, S.P. A Textbook on Heat Transfer (4th ed.). Universities Press. pp. 257–258, 2005.
5. Riedl, M. Optical Design Fundamentals for Infrared Systems (2nd ed.). SPIE Press, Bellingham, WA, 2001.
6. Holman, J. P. Heat Transfer Tenth Edition, McGraw-Hill. p. 111, Example 3–10 (2008).
7. Weisstein, E. W. Divergence Theorem, MathWorld-A Wolfram Web Resource. http://mathworld.wolfram.com/DivergenceTheorem.html.
8. Onsager, L. Physical Review. 37, 405–426, (1931); 38, 2265–2279, (1931).
9. Casimir, H. B. G. Reviews of Modern Physics. 17, 343–350, (1945).
10. Hahn, D. W. and Özişik, M. N. Heat Conduction, John Wiley & Sons, Inc., ch. 15 (2012).