Published: 05/30/2012, Last Updated: 05/30/2012

**By Dr. Michael J. Gourlay**

Download Fluid Simulation for Video Games (part 10) [PDF 1.1MB]

Download Fluid Simulation Source Code (part 10) [787KB]

Heat transfers by contact during collisions at the molecular level, and air expands when it heats. The expanded air occupies more volume but has the same mass; it therefore has lower density. Lower-density air buoys above higher-density air. Thus, hot air rises.

A previous article (Part 8) already explained how a vortex-based fluid simulation can handle variable density in a fluid. So, to get hot air to rise, the simulation must incorporate thermal effects: conduction, diffusion, and expansion.

This article-the tenth in a series-describes how density varies with temperature, how heat transfers throughout a fluid, and how heat transfers between bodies and fluid. Part 1 summarized fluid dynamics; Part 2 surveyed fluid simulation techniques, and Part 3 and Part 4 presented a vortex-particle fluid simulation with two-way fluid-body interactions that runs in real time. Part 5 profiled and optimized that simulation code. Part 6 described a differential method for computing velocity from vorticity, and Part 7 showed how to integrate a fluid simulation into a typical particle system. Part 9 described how to approximate buoyant and gravitational forces on a body immersed in a fluid with varying density.

This article introduces features to simulation code presented in previous articles: Now, bodies immersed within the fluid heat or cool surrounding fluid (and vice versa), fluid can transfer heat within itself, and fluid density depends on its temperature. These new features augment how visual effects have a two-way interaction with physical objects in the simulation and will segue into the next article, on combustion, in which a chemical process will heat fluid.

As mentioned in Part 8, this simulation employs the Boussinesq approximation to neglect compressibility effects-essentially, skipping sound waves, which are computationally expensive to simulate and do not usually affect fluid flow enough to worry about them. Although it assumes that fluid does not compress much because of pressure, the Boussinesq approximation does allow fluid density to change according to its temperature.

**Note:** *As mentioned in Part 8, the Boussinesq approximation has many limitations. To learn about improving that approximation, see "Further Reading" for a link to information about the anelastic and pseudo-incompressible approximations.*

*Thermal Expansion*

Hot air rises, because that air is less dense. The simulation needs to quantify how density relates to temperature.

Fluid density, , depends on temperature, *T*, according to this formula:

or

where is a reference density when the fluid is at some reference temperature . Solving this equation for leads to . Here, is the coefficient of thermal expansion, defined as

*Thermal expansion* describes how a material's volume changes with temperature, as depicted in Figure 2. Note that this formula is a definition (it's not derived from anything) and that is in general a function of temperature (not a constant).

To derive a formula for thermal expansion coefficient for an ideal gas, apply that formula to the ideal gas law,

(where is the amount of substance and *R* is the universal ideal gas constant) to get . Combined, this equation leads to the approximation . For the simulation code, use ambient fluid density for . (See part 8 for details about this so-called *ambient density*.)

This formula lets the simulation reuse density instead of having to keep track of fluid temperature as a separate property; temperature is a simple function of density and vice versa. Later in the article, you'll see how to implement this formula in the simulation code.

*Thermal Diffusion*

The fluid simulation presented so far has concerned itself with mechanical energy. In contrast, temperature is tied to the "internal energy," such as vibration and rotation of molecules. So now, the simulation needs an equation that expresses what happens to internal energy when temperature varies in a fluid. In particular, the simulation needs to know how heat moves around, even when the fluid does not (Figure 3).

The internal energy equation (a.k.a. *thermal energy equation*, a.k.a. *heat equation*) for a fluid describes how heat varies with time and space:

Here, *e* is the internal energy and is the heat flux-that is, the rate at which heat enters and leaves a tiny volume. And is viscous heating: When fluid molecules collide, they lose some of their kinetic energy to internal energy in the form of molecular vibration, rotation, and translation, also known as *heat*.

**Note:** *The internal energy equation comes from the first law of thermodynamics, which is an expression of conservation of energy. See "Further Reading" for links to details, including derivations and explanations of this formula.*

Given the thermal expansion formula from above in the form you can modify the volume compression term, as

Assume that the fluid is an ideal gas, with equation of state (You can use different equations of state for liquids, but for now assume an ideal gas.) With judicious application of Maxwell relations, you can prove that this assumption also implies that , where is specific heat capacity at constant pressure and is specific heat capacity at constant volume. (See, for example, Mandl p. 125). As mentioned earlier, also, for an ideal gas, . The internal energy equation for an ideal gas becomes:

Plug this equation back into internal energy equation to get

Viscous dissipation is important for the mechanical energy equation, and having in both places (with opposite signs) shows that kinetic energy lost equals thermal energy gained. (Because is always positive, heat energy flows only one way-from kinetic to thermal.) But for fluid velocities well below the speed of sound, you can neglect its contribution to heat. So drop . This leaves the terms .

**Note:** *Viscous dissipation is proportional to the square of velocity gradients, which can also be expressed as the mean of vorticity dotted with itself, which is called enstrophy. See, it all comes back to vorticity!*

**Fourier's Law of Heat Conduction**

To get any use out of the internal energy equation, you need an expression for heat flux, .

Within the material depicted in Figure 3, the molecules move around, jostling other molecules they contact. Each collision transfers internal energy, and this is called *conduction*.

Fourier's law of heat conduction states that the heat transfers in the opposite direction of the temperature gradient (that is, from hot to cold):

Here, *k* is conductivity (where ) and is the temperature gradient. Plug that into the internal energy equation above to get a diffusion equation:

You've seen this equation before-back in Part 3-except for a vector (vorticity) instead of a scalar (temperature). Naively computing the right-hand-side of that equation using the Laplacian () would require computing second-order spatial derivatives and doing a lot of bookkeeping that would complicate the simulation code. But the simulation will employ a trick-particle strength exchange (PSE)-to simplify the computation. (See the section "Computing Thermal Diffusion Within the Fluid," later in this article.)

**Note:** *Because the Boussinesq approximation relates temperature to density, heat diffusion implies mass diffusion, which in turn implies that mass diffuses the way salinity (or some other solute) would. Read up on Fick's laws of diffusion for more information about how diffusion relates to flux.*

To add the effects of thermal expansion and diffusion, the simulation must compute them. The simulation code requires these modifications:

- Treat vortex particles (vortons) as carriers of heat for the fluid.
- Compute thermal diffusion within the fluid using particle strength exchange (PSE).
- Compute thermal conduction between rigid bodies and the fluid.

**Computing Thermal Expansion**

From the modifications made in Part 8, vortons already carry density for the fluid. The thermal expansion formula derived above lets the simulation reuse density instead of having to keep track of fluid temperature as a separate property; temperature is a simple function of density and vice versa. Add the following methods to the `Particle`

class:

```
float GetTemperature( float ambientDensity ) const
{
return /* ambientTemperature* */ambientDensity / ( ambientDensity + mDensity ) ;
}
void SetTemperature( float ambientDensity , float temperature )
{
float density = ambientDensity /* *ambientTemperature */ / temperature ;
mDensity = density - ambientDensity ;
}
```

Simplify this code a step further by assuming that `ambientTemperature`

is 1 in whatever units the simulation uses, and omit that factor from the formulas in the code above. Technically, not even this change is necessary, because you can compute thermal diffusion and conduction directly using density, but it's illustrative to have these formulae in code.

Remember from Part 8 that the baroclinic generation of vorticity causes fluid with density lower than its surroundings to rise. So by simply lowering density, the fluid will rise "automatically."

**Computing Thermal Diffusion Within the Fluid**

As described earlier, this partial differential equation describes how heat diffuses through a fluid:

As described in Part 3, the simulation can compute diffusion by exchanging heat between nearest neighboring particles-so-called *particle strength exchange*. Doing so efficiently requires a fast spatial partition to find nearby vortons. The simulation already creates that partition for vortex diffusion, so you can reuse it for thermal energy diffusion.

To reuse the vorton partition, first extract the spatial partitioning into its own routine. This isn't new code; it's just in a new place:

```
void VortonSim::PartitionVortons( const float & timeStep , const unsigned & uFrame
, UniformGrid< Vector< unsigned > > & ugVortonIndices )
{
const size_t numVortons = mVortons->Size() ;
for( unsigned offset = 0 /* Start at 0th vorton */ ; offset < numVortons ; ++ offset )
{ // For each vorton...
Vorton & rVorton = (*mVortons)[ offset ] ;
// Insert the vorton's offset into the spatial partition.
ugVortonIndices[ rVorton.mPosition ].PushBack( offset ) ;
}
}
```

The simulation (`VortonSim::Update`

, to be precise) calls `PartitionVortons`

before calling `DiffuseVorticityPSE`

.

Next, the simulation needs to exchange heat, which it does in a routine called `DiffuseHeatPSE`

. The recipe follows exactly that used for `DiffuseVorticityPSE`

. The code is so similar, in fact, that this article only shows the "innards" of the routine, which is called inside an inner loop otherwise identical to that used to diffuse vorticity. (Of course, the accompanying code contains everything, so see it for details.)

```
void VortonSim::ExchangeHeat( const unsigned & rVortIdxHere , Vorton & rVortonHere
, float & rDensityHere , const unsigned & ivThere , const Vector< unsigned > & cell
, const float & timeStep )
{
const unsigned & rVortIdxThere = cell[ ivThere ] ;
Vorton & rVortonThere = (*mVortons)[ rVortIdxThere ] ;
float & rDensityThere = rVortonThere.mDensity ;
const float densityDiff = rDensityHere - rDensityThere ;
const float exchange = 2.0f * mThermalDiffusivity * timeStep * densityDiff ;
rDensityHere -= exchange ; // Make "here" temperature a little closer to "there".
rDensityThere += exchange ; // Make "there" temperature a little closer to "here".
}
```

**Note:** T*his routine directly exchanges density (not temperature). For an ideal gas, whose thermal expansion coefficient varies with temperature, this is technically improper; the quantity being exchanged is inversely proportional to temperature. But the result qualitatively the same: Particles exchange heat until they reach equilibrium. Furthermore, this formulation is correct for materials that have constant thermal expansion coefficients. Either way, for visual effects for video games, the distinction is negligible.*

**Computing Thermal Conduction Between Bodies and Fluid**

To transfer heat between rigid bodies and the fluid, add some code to the simulation that handles collisions between vortons and bodies, and exchange heat analogous to how the simulation transfers angular momentum between vortons and bodies. This code goes inside `FluidBodySim::SolveBoundaryConditions:`

```
// Conduct heat from body to fluid.
const float vortonTemperatureOld = rVorton.GetTemperature( ambientFluidDensity ) ;
const float temperatureDifference = rSphere.mTemperature - vortonTemperatureOld ;
const float heatConduction = temperatureDifference * rSphere.mThermalConductivity;
const float heatExchange = heatConduction /* * timeStep */ ;
const float vortonHeatCapacity = 1.0f ;
const float vortonTemperatureNew = vortonTemperatureOld + heatExchange * vortonHeatCapacity ;
rVorton.SetTemperature( ambientFluidDensity , vortonTemperatureNew ) ;
// Conduct heat from fluid to body.
rSphere.mTemperature += heatExchange * rSphere.mOneOverHeatCapacity ;
```

**Note:** *This routine takes a few shortcuts: The time step is ignored, and the heat capacity for a vorton is independent of its size. The rationale is that, for these simulations, the time step is fixed, and the vorton size is uniform and constant, so their effects are incorporated into mThermalConductivity. The code shows, in comments, where those values would be incorporated if for some reason you wanted to have variable time steps or vorton size, but be aware that such a change would likely require additional changes elsewhere.*

Computing `DiffuseHeatPSE`

costs some CPU time (see below to see exactly how much). To mitigate that cost, use Intel® Threading Building Blocks (Intel® TBB) to parallelize it.

Because the heat diffusion uses basically the same algorithm as for vorticity diffusion, simply reuse the recipe used to parallelize vortex diffusion:

- Write heat diffusion such that it can operate on arbitrary "slices" of data. In this case, slice the data across vortons:

```
void VortonSim::DiffuseHeatPSESlice( const float & timeStep
, const UniformGrid< Vector< unsigned > > & ugVortonIndices
, size_t izStart , size_t izEnd )
{ // Exchange heat with nearest neighbors:
// For each cell in the spatial partition...
// For each vorton in each cell...
ExchangeHeat( ... ) ;
}
```

- Create a functor to compute heat diffusion:

```
class VortonSim_DiffuseHeatPSE_TBB
{
float mTimeStep ; ///< Duration since last time step.
VortonSim * mVortonSim ; ///< Address of VortonSim object
const UniformGrid< Vector< unsigned > > & mUgVortonIndices ;
public:
void operator() ( const tbb::blocked_range<size_t> & r ) const
{ // Compute subset of heat diffusion.
mVortonSim->DiffuseHeatPSESlice(mTimeStep,mUgVortonIndices,r.begin(),r.end());
}
VortonSim_DiffuseHeatPSE_TBB( float timeStep , VortonSim * pVortonSim
, const UniformGrid< Vector< unsigned > > & ugVortonIndices )
: mTimeStep( timeStep )
, mVortonSim( pVortonSim )
, mUgVortonIndices( ugVortonIndices )
{}
} ;
```

- Use Intel TBB to invoke the functor:

```
void VortonSim::DiffuseHeatPSE( const float & timeStep , const unsigned & uFrame
, const UniformGrid< Vector< unsigned > > & ugVortonIndices )
{ // Exchange heat with nearest neighbors
const unsigned & nz = ugVortonIndices.GetNumPoints( 2 ) ;
const unsigned nzm1 = nz - 1 ;
#if USE_TBB
// Estimate grain size based on size of problem and number of processors.
const size_t grainSize = MAX2( 1 , nzm1 / gNumberOfProcessors ) ;
// Compute heat diffusion using threading building blocks
parallel_for( tbb::blocked_range<size_t>( 0 , nzm1 , grainSize )
, VortonSim_DiffuseHeatPSE_TBB( timeStep , this , ugVortonIndices ) ) ;
#else
DiffuseHeatPSESlice( timeStep , ugVortonIndices , 0 , nzm1 ) ;
#endif
#endif
}
```

Figure 4 shows a simulation with a hot ball embedded inside a fluid containing temperature-sensitive dye. Initially, the dye is blue, indicating cold fluid. As the ball gradually heats the fluid, and the heat diffuses through the fluid, and the dye turns red.

Figure 5 shows a continuation of the sequence in Figure 4: The heated fluid expands (so its density decreases); buoyancy therefore causes that fluid to rise.

Figure 6 shows a rendering of a simulation like that shown in Figures 4 and 5 but with only a thin sheet of dye instead of a thick block. This view more clearly reveals the formation of a pair of counter-rotating vortices that cause the fluid to rise. These vortices formed because of baroclinic generation (as described in Part 8).

**Performance**

What does all of this additional computation cost? Very little-less than 1 percent, in fact. Table 1 shows profile measurements.

**Table 1.** *Run durations for processes on a 3.4 GHz Intel® Core™ i7-2600 processor*

Threads | DiffuseHeatPSE (ms) | PartitionVortons (ms) | Total (ms) |

1 | 0.021 | 0.046 | 10.7 |

2 | 0.015 | 0.048 | 7.55 |

4 | 0.011 | 0.047 | 6.46 |

Note that `PartitionVortons`

was already running before this addition-it's been there since the first code accompanying any of these articles-but it used to be incorporated into the `DiffuseVortonPSE`

routine. Now, it simply runs in a separate routine. It's a useful benchmark, because it is not parallelized, so the runtimes should be roughly independent of the number of threads (and it is).

In contrast, `DiffuseHeatPSE`

uses Intel TBB and gets significant speed gains as the number of threads increases. But the total cost of `DiffuseHeatPSE`

is less than 0.2 percent of the total, regardless of the number of threads.

What did adding thermal expansion cost the simulation? Trick question: It's "free." Or rather, its price was paid back when fluid buoyancy was added and density incorporated into the simulation. Again, no code was added to make "thermal expansion" happen; you simply had the simulation interpret density as both mass-per-volume and an expression of temperature.

The next part in this series will build on this installment, adding *combustion* (generating heat by chemical processes). That addition will give the simulation the ability to simulate smoldering and burning.

Kundo, P.J. (1990). *Fluid Mechanics*. San Diego: Academic Press. See especially chapter 4, "Conservation Laws," section 17, "Boussinesq Approximation."

Mandl, F. (1988). *Statistics Physics*, 2nd Ed. Chichester: John Wiley and Sons. See especially chapter 5, "Simple Thermodynamic Systems," section 3, "The Difference of Heat Capacities," and section 4, "Some Properties of Perfect Gases."

Fluid Simulation for Video Games (part 1)

Fluid Simulation for Video Games (part 2)

Fluid Simulation for Video Games (part 3)

Fluid Simulation for Video Games (part 4)

Fluid Simulation for Video Games (part 5)

Fluid Simulation for Video Games (part 6)

Fluid Simulation for Video Games (part 7)

Fluid Simulation for Video Games (part 8)

Fluid Simulation for Video Games (part 9)

Fluid Simulation for Video Games (part 10)

Fluid Simulation for Video Games (part 11)

Fluid Simulation for Video Games (part 12)

Fluid Simulation for Video Games (part 13)

Fluid Simulation for Video Games (part 14)

Fluid Simulation for Video Games (part 15)

Fluid Simulation for Video Games (part 16)

Fluid Simulation for Video Games (part 17)

Fluid Simulation for Video Games (part 18)

Fluid Simulation for Video Games (part 19)

Dr. Michael J. Gourlay works as a Senior Software Engineer at Electronic Arts. He currently works as the software architect for the Football Sports Business Unit. He previously worked as a senior lead engineer on Madden NFL, on the procedural animation system used by EA, on Mixed Martial Arts (MMA), and as a lead programmer on NASCAR. He architected the visual effects system used in EA games worldwide and patented algorithms for interactive, high-bandwidth online applications . He also developed curricula for and taught at the University of Central Florida (UCF) Florida Interactive Entertainment Academy (FIEA) , an interdisciplinary graduate program that teaches programmers, producers and artists how to make video games and training simulations. Prior to joining EA, he performed scientific research using computational fluid dynamics (CFD) and the world's largest massively parallel supercomputers. His previous research also includes nonlinear dynamics in quantum mechanical systems, and atomic, molecular and optical physics . Michael received his degrees in physics and philosophy from Georgia Tech and the University of Colorado at Boulder.

Intel's compilers may or may not optimize to the same degree for non-Intel microprocessors for optimizations that are not unique to Intel microprocessors. These optimizations include SSE2, SSE3, and SSSE3 instruction sets and other optimizations. Intel does not guarantee the availability, functionality, or effectiveness of any optimization on microprocessors not manufactured by Intel. Microprocessor-dependent optimizations in this product are intended for use with Intel microprocessors. Certain optimizations not specific to Intel microarchitecture are reserved for Intel microprocessors. Please refer to the applicable product User and Reference Guides for more information regarding the specific instruction sets covered by this notice.

Notice revision #20110804