Heat Transfer

Introduction

This tutorial gives an introduction to modeling heat transfer. 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 behind a heat transfer are 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 is 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 as well as 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 is 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 the this notebook requires. The downside is that the visual quality of the animations will not be as crisp as without it. To obtain high quality graphics remove or comment out the call to Rasterize.

To get high fidelity visualizations comment out the rasterization process.

The symbols and corresponding units used throughout this tutorial are summarized in the Nomenclature section.

Load the finite element package.

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.

Beside 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 may very well depend on the temperature itself. This will then result in a nonlinear heat equation.

A 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 from 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 density and is the internal energy per unit mass. The total energy within the control volume is then equal to the product . The red circle in the middle represents a heat source , which denotes thermal energy that is generated inside the domain . The heat flux represents the net energy that exits through the boundaries.

The unit of the heat flux depends on the dimension of the boundary. In a 3D domain the boundary is represented by 2D surfaces, and the heat flux has a unit of . In a 2D domain, however, the boundary has 1D lines, which makes the unit of become . In 1D the boundary is represented by a dimensionless point and the heat flux has a unit of .

The energy balance within the domain can be described by the following equation:

That is, the change in the total energy is equal to the energy generated inside the domain minus the net 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 proportional 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 Flourier's law of heat conduction:

Here the energy diffusivity is represented by the thermal conductivity , 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 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:

or

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 we consider a 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 part is understood as the accumulation (or loss) of mass in the domain over time, while the divergence term part denotes the difference of the mass inflow versus outflow.

When there is no internal mass generation/destruction, the rate of the mass inflow is equal to the rate of the mass outflow plus the accumulation of mass within the domain . That is, the term sums up to zero and can be removed from the equation (6) leading to:

This leads to a restriction on the usage of the heat transfer model: Equation (7) and the resulting heat transfer PDEs (8) and (9) can not 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:

PDE for heat transfer in a medium with flow.

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.

In a solid medium, however, the internal velocity field is set to zero and the governing PDE simplifies to a pure conductive heat equation:

PDE for heat transfer in a solid medium.

Note that the preceding heat equation (13) is written in a non-conservative form. That is, the term and with it the density and the velocity are outside of the divergence operator. One could imagine that 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.

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 detail 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 an 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
Define a stationary heat transfer model function.
Set up a 1D static heat transfer model.
Define a time dependent heat transfer model function.
Set up a 1D time dependent heat transfer model.

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.

Set up model parameters for the domain and the simulation end time .
To make use of specific material parameters in the equation we extract relevant data from the ThermodynamicData.

In some examples we will be using a smoothed step function 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 .

Define and visualize the smoothed step function .

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.

Our goal is to find the steady-state temperature distribution of the ceramic strip.

Set up a rectangular domain with a width of and a height of .

The thermal conductivity , heat transfer coefficient , density , heat capacity and emissivity of the ceramic trip are given by:

Define the model parameters.
Set up temperature surface boundary conditions at the left and right boundaries.
Set up the convective boundary condition on the top surface.
Set up the thermal radiation boundary condition on the top surface.

A default thermally insulated boundary condition is implicitly applied on the remaining bottom boundary.

Solve the heat transfer PDE model.
Visualize the steady-state temperature distribution.

In the steady-state, the minimum temperature is found on the top surface that is cooled by both convection and radiation. 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.

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 , either by explicitly generating the mesh for them or by making use of the MeshRefinementFunction.

Volumetric Heat Source

A volumetric heat source can be used to model an arbitrarily shaped heat source () or heat sink () within the domain. The corresponding source strength denotes the rate of an internal heating or cooling per unit volume.

The term volumetric heat source may be a bit misleading. This is best seen by looking at the unit of a volumetric heat source that depends on the dimension of the system. In a 1D domain (), 2D domain () and in a 3D domain () has a unit of , and , respectively. So name comes from the 3D incarnation of the heat source but is used in other dimensions as well.

In the following 2D example a rectangular heat source is introduced to heat up the domain. The source strength is fixed at .

Define a 2D domain and a rectangular source region.
Set the source strength .

For more involved shapes of a heat source, the function RegionMember can be used to specify the source region .

Set up the heat source with a RegionMember function.

In this case, however, because the shape of the heat source is simple we can simply specify the formula of the source region by making use of an If statement. This will lead to a more efficient time integration. The reason the If statement is more efficient is that it can be automatically compiled while the RegionMemberFunction can not be compiled.

Set up a simplified heat source .

The most accurate and efficient method to deal with heat sources is by element makers 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 makers for the heat source set up is presented in the appendix section: Modeling Heat Source by Element Markers.

Define the heat transfer PDE with a volumetric source term and an initial temperature field .
Solve the PDE with NDSolveValue.
Set up a legend bar and ContourPlot options for the temperature field plot.
Make an animation of the solution using Plot and ListAnimate.

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.

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. The corresponding source strength denotes the rate of the internal heating or cooling per unit area.

Here is a Dirac delta function at the source location .

The Dirac delta function, however, poses a problem in numerical simulations as it can not be resolved in the discretized spatial domain. This is because the Dirac delta function is singular at the source location . 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]. In this tutorial we choose:

where is the regularization parameter that controls the support of the regularized delta functions . Typically should have a size comparable to the mesh spacing .

Set up the regularized delta functions centered on the source region.

The units of a layer heat source depend on the dimension of the system, much like the volumetric heat source. In a 3D domain () the layered heat source is presented as a 2D surface, and thus has a unit of . In 2D domain (), however, the geometry of a layer heat source renders into a 1D line, and the unit of becomes . In 1D there is no layer heat source. In that case a point heat source should be used instead.

In the following 2D example a layer heat source is added at to heat up the domain. The source strength is fixed at .

Define a 2D domain and the source region at .

To accommodate the layered heat source in the mesh, a boundary mesh with the line at is generated.

Generate and visualize the mesh that contains the source region.

More information on the mesh generation can be found in the Element Mesh Generation tutorial.

To utilize the regularized delta function , we choose the regularization parameter to be half of the mesh spacing: .

Inspect the mesh spacing and set the regularization parameter .
Set the layer heat source with the strength using the regularized delta function .
Define the heat transfer PDE with a layer heat source and an initial temperature field .
Solve the PDE with NDSolveValue.
Set up a legend bar and ContourPlot options for the temperature field plot.
Make an animation of the solution using Plot and ListAnimate.

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 at , the thermal energy is generated and propagates towards both sides of the domain. The speed of the heat transfer depends on the heat conductivity and the heat capacity of the material.

Note that a layered heat source is not limited to a linear surface () or a straight line (). The following 2D example demonstrates a layer heat source with a curved source region.

Define a 2D domain and a curved source region.
Generate meshes of boundary and the curved source region.
Join the two meshes.
Generate and visualize the mesh that contains the curved source region.

More information on the mesh generation can be found in the Element Mesh Generation tutorial.

Set the layer heat source with a curved source region.
Solve the PDE with NDSolveValue.
Set up a legend bar and ContourPlot options for the temperature field plot.
Make an animation of the solution using Plot and ListAnimate.

See this note about improving the visual quality of the animation.

With a curved layered heat source heading downward, most of the heat is transported toward the lower half of the domain.

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. The corresponding source strength denotes the rate of the internal heating or cooling at the source point in unit power.

Since the point heat source has no spatial extension in all directions, the Dirac delta function should be applied on each dimension (i.e. ) of the modeling domain :

Set up the regularized delta functions centered at the source point .

In the following 2D example a point heat source is added at to heat up the domain. The source strength is fixed at .

Generate the mesh that contains the source point .
Set the point heat source with the strength using the regularized delta function .
Define the heat transfer PDE with a point heat source and an initial temperature field .
Solve the PDE with NDSolveValue.
Set up a legend bar and ContourPlot options for the temperature field plot.
Make an animation of the solution using Plot and ListAnimate.

See this note about improving the visual quality of the animation.

The simulation begins with an undisturbed domain where . With a point 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 we 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 with a different rate. The diffusion term (19) is then rewritten as:

or

where is the thermal conductivity tensor. and are called the principle 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 principle directions , and and the off-diagonal conductivity coefficients are zero. This behaviour 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 the orthotropic case where the heat transfer is more efficient horizontally. Should the thermal conductivity tensor be described by:

Set the thermal conductivity tensor of the orthotropic material.
Solve the heat transfer PDE model.
Set up a legend bar and ContourPlot options for the temperature field plot.
Visualize the heat transfer in an orthotropic material.

See this note about improving the visual quality of the animation.

Unlike the example shown in the previous section, in the case the heat transfer is faster in the horizontal direction, resulting in a higher temperature zone within .

Nonlinear Heat Transfer

In the previous sections we 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.

Set up a heat flux boundary condition at the left end with NeumannValue.
Set the non-linear heat transfer PDE with a temperature-dependent thermal conductivity and an initial temperature field .
Solve the nonlinear PDE using NDSolveValue.

To understand the effects of the nonlinearity compare to a linear heat transfer PDE.

Set up and solve a linear heat transfer PDE with a constant thermal conductivity .
Make an animation of the solutions using Plot and ListAnimate.

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 non-linear 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.

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:

Under these four types, the following boundary conditions are introduced:

The following section describes several physical boundaries common in heat transfer and how they can be modeled with the use of DirichletCondition, NeumannValue and PeriodicBoundaryCondition. For this purpose the boundary condition currently discussed is always on the left hand side of the simulation domain. In some examples additional boundary conditions are on the right hand side to better demonstrate the behaviour of the boundary condition on the left hand side.

Temperature Surface Boundary Condition

Purpose

The purpose of a temperature surface 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:

A temperature surface boundary for the dependent variable modeled with DirichletCondition.

Derivation

We speak of a temperature boundary condition when the surface temperature is prescribed on a boundary. The surface temperature can either be 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 described the profile of the surface temperature from to . The parameters and are arbitrarily chosen to simulate the heating process.

Specify the surface temperature profile with a smoothed step function.
Set up the temperature boundary condition with DirichletCondition.
Set up the PDE with an undisturbed initial temperature field: .
Solve the PDE with NDSolveValue.
Make an animation of the solution using Plot and ListAnimate.

The simulation begins with an undisturbed domain where . As the surface temperature increases at the left boundary, the excess thermal energy is then passed to the right and brings up the 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 amount of thermal energy flowing into or out of some part of the boundary.

Formulation

With a prescribed heat flux on the boundary , the heat flux boundary condition is given by:

A heat flux boundary for the dependent variable modeled with NeumannValue.

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 indicates 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 (23) relates the heat flux with the temperature gradient :

Inserting (24) into (25), the heat flux boundary condition can be written as:

Note that the unit of a 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:

Set up the heat flux boundary at the left end with NeumannValue.
Set up the PDE with an undisturbed initial temperature field: .
Solve the PDE with NDSolveValue.
Make an animation of the solution using Plot and ListAnimate.

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 . The uneven temperature field is then smoothed out with time by the internal heat transfer.

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:

A thermally insulated boundary for the dependent variable modeled with NeumannValue.

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 (26) into the heat flux boundary condition (27), 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.

To model a thermally insulated boundary condition the NeumannValue is set to at the left end.

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.

Set up a constant heat flux on the right hand boundary.
Set up the PDE with an undisturbed initial temperature field: .
Solve the PDE with NDSolveValue.
Make an animation of the solution using Plot and ListAnimate.

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:

A symmetry boundary for the dependent variable modeled with NeumannValue.

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 sub-domain 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 , we can efficiently construct the simulation domain 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 time, 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:

An outflow boundary for the dependent variable modeled with NeumannValue.

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:

A convective boundary for the dependent variable modeled with NeumannValue.

Derivation

With the existence of a fluid flowing 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 are presented in following table:

Several empirical based formulas are also built [28,29] to estimate the heat transfer coefficient in different situations.

Inserting (30) into the heat flux boundary condition (31), 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 .

Define the model parameters.
Set up the convective boundary condition with NeumannValue.
Set up the PDE with an undisturbed initial temperature field: .
Solve the PDE with NDSolveValue.
Make an animation of the solution using Plot and ListAnimate.

With the external flow flowing 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:

A thermal radiation boundary condition for the dependent variable modeled with NeumannValue.

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 (32) we can formulate a radiation boundary condition 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 the Kirchhoff's law of thermal radiation [33].

The equation (34) is then simplified to:

Assuming that the ambient surrounding behaves as a black body with the emissivity , the equation can be further simplified as:

Inserting (35) into the heat flux boundary condition (36), then the radiative boundary condition is given by:

or

Note that the above derivation is performed based on the absolute temperature. That means the temperature terms in (37) have the unit of Kelvin .

To apply a radiation boundary condition in Celsius , an unit conversion should be done in (38) :

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.

Define the model parameters.
Set up the thermal radiation boundary condition in the unit of Celsius .
Set up the PDE with an undisturbed initial temperature field: .
Solve the PDE with NDSolveValue.
Make an animation of the solution using Plot and ListAnimate.

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 as 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: .

Inspect the temperature gradient on the left hand boundary at .

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:

A periodic boundary condition for the dependent variable modeled with PeriodicBoundaryCondition.

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 .

Set up the PeriodicBoundaryCondition at the right end with the mapping function .

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 heat transfer model.

Specify the temperature profile of the heating film at the left end.
Set up the PDE with an undisturbed initial temperature field: .
Solve the PDE with NDSolveValue.
Make an animation of the solution using Plot and ListAnimate.

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 replaces by using the DirichletCondition from the left boundary also at the right boundary

Alternative setup to replace the periodic boundary condition.

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 to the target.

Appendix

Special Cases of the Heat Equation

Stationary case

If the temperature field is in a steady state, the transient term in (39) 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 . The heat equation may also be expressed using a cylindrical or spherical coordinate system.

A graphics 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:

or

By inserting the coordinates relation (40) into (41), 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 (42) 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. An example making use of the cylindrical coordinate system can be founded in the separate Heat Transfer Verification Tests notebook.

Heat Equation in Spherical Coordinates

A graphics 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:

or

By inserting the coordinates relation (43) into (44), the heat equation can be expressed in the spherical coordinate system as:

The Smoothing Characteristic of the Diffusion Equation

The basic behaviour of a diffusion equation is smoothing. To see the effect set up initial conditions that are discontinuous at .

The PDE for the temperature field is given by:
Solve the PDE with NDSolveValue.
Visualize the initial conditions and the smoothed out solution.

In a short time the sharp tip of the discontinuous initial condition at x=1/2 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. Element marker 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 we introduce a rectangular heat source to heat up the domain.

Create and visualize a boundary element mesh with an internal boundary separating the heat source from the remaining domain.

Next we specify region markers with the "RegionMarker" option for ToElementMesh. 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.

Create and visualize an element mesh with internal markers.

An alternative to generate the mesh is to make use of the Boolean region functions and specify the no region holes should be inserted. The same region marker process from above is used.

Specify a region.
Set up the mesh.
Set the source strength .
Set up a heat source with the element marker.
Define the heat transfer PDE with a volumetric source term and an initial temperature field .
Solve the PDE with NDSolveValue.
Set up a legend bar and ContourPlot options for the temperature field plot.
Make an animation of the solution using Plot and ListAnimate.

See this note about improving the visual quality of the animation.

This result agrees with the one that didn't make use of element makers, which is shown in the section: Volumetric Heat Source. For large complicates 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 field, we 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. We can then apply the divergence theorem [45] on the equation (46) leading to:

Since the domain is completely arbitrary, we can discard the integral and yield the mass conservation equation in the differential form:

That means the mass conservation equation (47) also holds true for the domain with discontinuous density/velocity field.

The heat equation is essentially an energy conservation equation, and can be derived in a similar way by substituting the density in (48) with the internal energy . Therefore, the heat transfer model presented in the tutorial can be applied on both continuous and discontinuous density/velocity fields.

Nomenclature

References

1.  Bilbao, Stefan and Hamilton, Brian. Directional Source Modeling In Wave-Based Room Acoustics Simulation. IEEE, 2017.

2.  Peskin, Charles. The Immersed Boundary Method. Cambridge University, 2002.

3.  Churchill, Stuart W. and Chu, Humbert H.S. Correlating equations for laminar and turbulent free convection from a vertical plate. International Journal of Heat and Mass Transfer, 18 (11): 13231329, 1975.

4.  Sukhatme, S.P. A Textbook on Heat Transfer (Fourth ed.). Universities Press. pp. 257258, 2005.

5.  Riedl, M. Optical Design Fundamentals for Infrared Systems (Second ed.). SPIE Press, Bellingham, WA, 2001.

6.  Holman, J. P. Heat Transfer Tenth Edition, McGraw-Hill. pp. 111, Example 3-10 (2008).

7.  Weisstein, Eric W. Divergence Theorem, MathWorld-A Wolfram Web Resource. http://mathworld.wolfram.com/DivergenceTheorem.html.

8.  Onsager, L. Physical Review. 37, 405426, (1931); 38, 22652279, (1931).

9.  Casimir, H. B. G. Reviews of Modern Physics. 17, 343350, (1945).

10.  Hahn, D. W. and Özişik, M. N. Heat Conduction, John Wiley & Sons, Inc, ch.15 (2012).