Heat transfer in fluid flows (i.e., convective heat transfer) is relevant for a wide range of practical applications including the conventional processing industry, pharmaceutical devices, thermal-storage units for sustainable energy systems, and high-precision equipment. Convective heat transfer is traditionally visualized as scalar maps of temperature fields and investigated with heat transfer coefficients characterizing the heat transfer rate at fluid-solid interfaces [1], [2]. However, this view offers only limited insight into the actual physical transport mechanisms: temperature fields reveal only the thermal distribution at a particular instant, not where heat is moving to or where it originated, and heat-transfer coefficients provide no insight into thermal transport pathways. Meaningful visualization of convective heat transfer would instead reveal the thermal paths themselves, the coherency of transport along them, and the structures that emerge from this transport.

Speetjens proposed an alternative Lagrangian approach to convective heat transfer [3] that relies on the notion that convective heat transfer fundamentally is the transport of thermal energy by a net heat flux encompassing the combined effect of fluid motion and molecular diffusion. Convective heat transfer thus becomes the transport of thermal energy along certain paths (“thermal paths”), in a similar way as fluid motion is the transport of fluid parcels by the flow along fluid paths. This mass-transfer analogy admits visualization of convective heat transfer by the topology of such thermal paths using well-established Lagrangian methods and concepts from flow visualization and analysis of mixing by the topology of fluid paths [4]. However, two essential differences between fluid and thermal transport render the utilization of existing Lagrangian methods for thermal transport non-trivial. Namely, heat transport is (i) non-conservative (i.e., net heat flux is divergent), and (ii) aperiodic, whereas fluid transport is conservative and generally periodic in time. Poincaré maps require periodicity that thermal transport does not have. Lagrangian coherent structure techniques [5] can capture transport structures in aperiodic fluid transport but do not reveal the path along which the transport occurs. Their physical interpretation in the context of thermal transport remains an open question. Existing flow visualization techniques such as Line Integral Convolution and the density-based method of Park et al. [6] can reveal flow behavior along with attracting and repelling structures, but only in steady flows, resulting in a view of the instantaneous structures in a time-frozen unsteady flow rather than the time-evolving transport structures.

In this work, we take a first step toward capturing coherent structures in convective heat transport by presenting a visualization technique for aperiodic divergent transport that simultaneously captures finite-time density structures and coherency in transport paths. The method works by advecting massless particles within a time segment of the spacetime formulation of thermal transport following from the Lagrangian transport formalism. Contributions of particles are then accumulated via a Gaussian kernel along path segments on an accumulation grid, producing both still images and animations revealing coherency of thermal transport routes and density evolution.

1 Background: Lagrangian thermal formalism↩︎

Our work uses the Lagrangian thermal formalism (LTF) proposed by Speetjens [3] to extract the Lagrangian view of the convective contribution of heat transport from standard computational fluid dynamics (CFD) derived temperature evolutions. In this section, we briefly introduce the formalism. The inputs are the flow field \(\boldsymbol{u}(\boldsymbol{x},t)\) and temperature field \(T(\boldsymbol{x},t)\), where \(\boldsymbol{x} \in \mathcal{D}\) denotes position in spatial domain \(\mathcal{D}\) and \(t\) denotes physical time. All subsequent fields are functions of \((\boldsymbol{x},t)\); arguments are dropped for brevity. The first step of the LTF involves the decomposition of the temperature field as \[\begin{align} T(\boldsymbol{x},t) = \widetilde{T}(\boldsymbol{x},t) + T'(\boldsymbol{x},t), \label{Temperature1} \end{align}\tag{1}\] with \(\widetilde{T}\) the temperature evolution in the case of a stagnant fluid, i.e., only due to conduction, \(\boldsymbol{u}=\boldsymbol{0}\). The stagnant temperature evolution is readily available via CFD simulations. Contribution \(T'\), by definition, incorporates the impact of the flow \(\boldsymbol{u}\), i.e., \(T'\neq0\) only if \(\boldsymbol{u}\neq\boldsymbol{0}\). This contribution is governed by the transport equation \[\begin{align} \frac{\partial T'}{\partial t} + \boldsymbol{\nabla}\cdot\boldsymbol{Q}' = F, \label{Convective1} \end{align}\tag{2}\] revealing that \(T'\) emanates from the interplay of two fundamental transport mechanisms: net heat flux \(\boldsymbol{Q}' = \boldsymbol{u}T' - \alpha\boldsymbol{\nabla}T'\) and source \(F = -\boldsymbol{u} \cdot \boldsymbol{\nabla}\widetilde{T}\). The net heat flux \(\boldsymbol{Q}'\), in turn, represents the combined heat transport due to the interplay of convection (\(\boldsymbol{u}T'\)) and conduction (\(- \alpha\boldsymbol{\nabla}T'\)) relative to the reference temperature \(\widetilde{T}\).

Transport equation 2  [1] is the backbone of the visualization and Lagrangian analysis of convective heat transfer enabling its description as the motion of a fluid. Substitution \[\begin{align} T'= \rho,\quad \boldsymbol{Q}'=\boldsymbol{M},\quad F = F_\rho, \label{MassFlux1} \end{align}\tag{3}\] with \(\rho\) the fluid density, \(\boldsymbol{M}\) the mass flux and \(F_\rho\) a volumetric mass source/sink (e.g., due to chemical reactions), translates 2 into the conservation law for mass [1]. The discussion below is in terms of this analogy; “fluid/flow” properties \((\rho,\boldsymbol{M},F_\rho)\) hereby represent their thermal counterparts \((T',\boldsymbol{Q}',F)\) following 3 .

Convective heat transfer in the current LTF corresponds with the motion of fluid particles, density \(\rho\), at velocity \(\boldsymbol{v} = \boldsymbol{M}/\rho\) along trajectories \(\boldsymbol{x}(t)\) described by the kinematic equation and associated formal solution \[\begin{align} \frac{d\boldsymbol{x}}{dt} = \boldsymbol{v}\left(\boldsymbol{x}(t),t\right) \quad\Rightarrow\quad \boldsymbol{x}(t) = \boldsymbol{\Phi}_t(\boldsymbol{x}_0),\nonumber\\ \boldsymbol{\Phi}_t(\boldsymbol{x}_0)\equiv \boldsymbol{x}_0 + \int_0^t\boldsymbol{v}\left(\boldsymbol{x}(\eta),\eta\right)d\eta, \label{KinEq1} \end{align}\tag{4}\] where \(\boldsymbol{x}_0 = \boldsymbol{x}(0)\) denotes the initial position. The thermal equivalent of \(\boldsymbol{v}\) corresponds to \(\boldsymbol{Q}'/T'\) which, due to \(T'\) having both positive and negative regions, is numerically unfavorable to integrate where \(T'\) approaches zero. Kinematic equation 4 has an equivalent formulation

\[\begin{align} \frac{d\boldsymbol{x}}{d\xi} = \boldsymbol{M}\left(\boldsymbol{x}(\xi),t(\xi)\right) \quad\Rightarrow\quad \boldsymbol{x}(\xi) = \boldsymbol{\Phi}_\xi(\boldsymbol{x}_0), \label{KinEq2} \end{align}\tag{5}\] with fictitious (\(\xi\)) and real (\(t\)) time relating via \[\begin{align} t = f(\xi) = \int_0^\xi \rho\left(\boldsymbol{x}(\eta),\eta\right)d\eta. \label{TimeAxis1} \end{align}\tag{6}\] This equivalence enables determination of the trajectories described by 4 via 5 , without numerically dividing by \(T'\). Important to note is that 2 and 3 yield \[\begin{align} \frac{d\rho}{dt}\equiv\frac{\partial \rho}{\partial t} + \boldsymbol{v}\cdot\boldsymbol{\nabla} \rho = F_\rho - \rho\boldsymbol{\nabla}\cdot\boldsymbol{v}, \label{NonConservative} \end{align}\tag{7}\] meaning that the generic \(F_\rho\neq 0\) implies non-constant particle density (i.e., \(d\rho/dt\neq 0\)) and, inherently, non-divergence-free flow (i.e., \(\boldsymbol{\nabla}\cdot\boldsymbol{v}\neq 0\)). This renders system 4 non-conservative.

2 Method↩︎

2.1 Overview↩︎

The structures emerging from the LTF remain largely unexplored due to the lack of a suitable visualization mechanism. We propose a method that captures these transport structures and paths in two-dimensional unsteady thermal flows governed by kinematic equation 5 . Our method advects massless particles along the spacetime representation of the thermal flow following from the LTF (illustrated in Figure 1). The inputs are identical to those of the LTF: A time-dependent temperature field with flow \(T\) and without flow \(\widetilde{T}\), from which \(T', \boldsymbol{Q}'\) and \(\boldsymbol{F}\) are derived (see equations 1 and 2 ). Particles are seeded within a temporal window of the domain (see Section 2.3), yielding attracting structures (heat sinks) under forward time integration and repelling structures (heat sources) under backward integration. The particle paths are visualized by accumulating onto a fixed grid via a Gaussian kernel applied to particle path segments (see Section 2.4). This allows finite-time density structures and coherency of transport paths to be visualized either as an image or animation (see Section 2.5).

2.2 Spacetime Integration↩︎

While \(\boldsymbol{Q'}\) is unsteady in the 2D physical domain, lifting \(t\) as a third phase-space coordinate yields an autonomous system \(\boldsymbol{Q}'_*= (Q'_x,\, Q'_y,\, T')\) in the spacetime domain. Massless particles are integrated along this steady thermal flow via fictitious time \(\xi\), i.e., \(\frac{d}{d\xi}(x,y,t) = \boldsymbol{Q}'_*\), enabling visualization of transport along temporal slices of the spacetime domain. The reparametrization produces trajectories equivalent to those of the thermal flow while avoiding integration issues near \(T' = 0\), as shown in Section 1. We integrate numerically using fourth-order Runge-Kutta on a trilinear interpolation of \(\boldsymbol{Q}'_*\). Each particle \(i\) is tracked by its spacetime position \(\mathbf{x_*}^n_i = (\boldsymbol{x}^n_i, t^n_i)\), where \(\boldsymbol{x}^n_i \in \mathcal{D}\) is the spatial position and \(t^n_i\) is the physical time at simulation step \(n\). Forward (\(\Delta\xi > 0\)) and backward (\(\Delta\xi < 0\)) integration reveal attracting and repelling structures respectively. Note that the direction of \(\xi\) does not necessarily coincide with that of physical time \(t\), since \(dt/d\xi = T'\), so regions where \(T' < 0\) cause forward \(\xi\)-integration to decrease \(t\). This signifies a negative heat flux backward in time, which is equivalent to a positive heat flux forward in time. Particles move at varying speeds through \(t\), affecting density evolution. This cannot be resolved by adaptive time steps of \(\Delta\xi\) without reintroducing numerical issues near \(T' = 0\); instead, we address the density bias through our seeding strategy.

2.3 Seeding↩︎

Seeding of particles has three goals: well-distributed particles for path visualization, physically meaningful density evolution, and visual insight into transport structures. These goals conflict. Seeding proportional to \(F\) captures physical density evolution but leaves regions of low \(|F|\) sparse, obscuring transport paths. Uniform seeding distributes particles well initially, but introduces domain-shape bias over time as particles migrate toward transport routes, leaving much of the domain empty and causing artificial sharp transitions to appear [6].

We address these conflicts as follows. Particles are seeded uniformly over the spacetime domain using a Halton sequence for quasi-random seed generation, where each particle \(i\) is assigned weight \(w_i\) initialized to \(F\) at its initial position and time, ensuring well-distributed particles. The weighting provides an initial importance measure used in the subsequent accumulation to improve the visibility of transport structures. To maintain visibility of transport structures over longer integration times, particles are reseeded at a rate proportional to elapsed fictitious time. The constant particle influx and uniform removal stochastically stabilize the particle distribution, preventing particles from collapsing towards stagnation regions. While reseeding sacrifices physical accuracy of the density evolution, it reveals a wider range of transport structures that would otherwise become invisible at long integration times. We specify seeding to occur over a finite time seeding window \([t_0, t_1]\) (see Figure 1) to allow the capture of finite-time structures.

2.4 Path Segment Accumulation↩︎

Particle contributions within time \(\tau\) are accumulated on a two-dimensional regular grid \(A(\boldsymbol{p})\), where \(\boldsymbol{p} \in \mathcal{D}\), using a 2D spatial Gaussian kernel. Similar to Kernel Density Estimation [7] and Smoothed Particle Hydrodynamics [8] techniques, a Gaussian kernel allows for the estimation of a density field from limited samples. The kernel is isotropic with width \(\sigma\), providing a tradeoff between smoothness and detail of the density field.

Accumulating discrete particle positions within \(\tau\) results in gaps when the particle trajectory segments defined by \(\Delta\xi\cdot|\boldsymbol{Q}'|\) exceed the size of a sample grid cell. Other work in steady flow visualization resolves this by normalizing the vector field [6], [9], [10]. However, this is unsuitable here, as the varying divergence of the transport fields is physically meaningful. Instead, we accumulate segments \([\mathbf{x_*}^{n-1}_i, \mathbf{x_*}^{n}_i]\) without modification of the transport field rather than discrete positions. The contribution of a segment to a sample position \(\boldsymbol{p}\) of the regular grid is defined by: \[a(\boldsymbol{p}, \boldsymbol{x}^{n-1}_i, \boldsymbol{x}^n_i) = f(\mu_i) \cdot \left(1 - \frac{|t_i^c - t_\text{slice}|}{\tau}\right)\cdot e^{-\frac{||\boldsymbol{x}^c_i - \boldsymbol{p}||^2}{2\sigma^2}}\] where \(\boldsymbol{x}^c_i\) is the closest spatial point along the segment to \(\boldsymbol{p}\) with the corresponding physical time \(t_i^c\), see Figure 1. The formula consists of three terms: the contribution at \(\boldsymbol{p}\) from the Gaussian kernel centered at the closest segment point \(\boldsymbol{x}^c_i\) ; a time factor that reduces contribution for particle positions proportional to the difference of the physical time \(t_i^c\) from the visualized time slice \(t_\text{slice}\); and a fade-in factor \(f(\mu_i)\), a smoothly increasing function from 0 to 1 over the particle lifetime \(\mu_i\) (e.g., smoothstep), preventing abrupt contributions from newly reseeded particles.

Accumulating segments directly produces overlap artifacts, since the accumulation effectively adds tube segments with varying opacity. The beginning of each segment overlaps with the end of the previous one, resulting in double accumulation at segment boundaries. This can be resolved by subtracting the contribution of the previous segment \([\mathbf{x_*}^{n-2}_i, \mathbf{x_*}^{n-1}_i]\) from that of the current segment. However, full subtraction eliminates contributions from slow-moving particles. While this is beneficial for path visualization, it would reduce visibility of density structures. We therefore subtract proportionally to \(l = 2\sigma^{-1}|\boldsymbol{x}^{n-1}_i - \boldsymbol{x}^n_i|\), reducing the subtraction for slow-moving particles based on segment length relative to the kernel width \(\sigma\). The total accumulation \(A(\boldsymbol{p})\) sums contributions of all particle path segments within \(\tau\), weighted by \(w_i\) as defined in Section 2.3. The accumulation field \(A(\boldsymbol{p})\) is defined as \[A(\boldsymbol{p}) = \sum^{m}_{i=1} w_i \cdot \max \left( a(\boldsymbol{p}, \boldsymbol{x}^{n-1}_i, \boldsymbol{x}^{n}_i) - l \cdot a(\boldsymbol{p}, \boldsymbol{x}^{n-2}_i, \boldsymbol{x}^{n-1}_i), 0\right)\] where \(m\) is the total number of particles. Since the spacetime formulation following from the LTF is autonomous, the accumulation field at a fixed temporal window defined by \(t_\text{slice}\) and \(\tau\) captures the same physical dynamics regardless of fictitious integration time \(\xi\). Blending between successive frames therefore accumulates contributions from particles at different stages of their fictitious integration, all depicting the same physical time slice. This is achieved by applying an exponential decay rate \(d\) at each frame, \[A(\boldsymbol{p}) \gets A(\boldsymbol{p})/(1+d),\] so that older segment contributions fade over time, revealing transport paths. When rendered as an animation, the fading path reveals directionality without altering the underlying physical dynamics.

2.5 Rendering↩︎

The accumulation field \(A(\boldsymbol{p})\) is scaled by a global factor and rendered to a texture by applying a colormap as a transfer function. For performance, segment contributions are computed only over samples within the spatial bounding box of each segment, avoiding iteration over the full accumulation grid. The accumulation field is rendered to a texture at a resolution that is an integer multiple of the accumulation grid, allowing bilinear interpolation to smooth the output without additional computation. The full visualization method is implemented in FlowExplainer flowexplainer?.

3 Results↩︎

3.1 Dataset↩︎

We consider a 2D time-dependent convective heat transfer case introduced in [3]. The domain is \(\mathcal{D} = [0,1] \times [0,\frac{1}{2}]\) with periodic boundary conditions in the \(x\)-axis. Heat transfer is driven by a temperature difference between a “hot” bottom wall and a “cool” top wall, with a time-periodic solenoidal velocity field consisting of two adjacent counter-rotating vortices undergoing horizontal oscillation, with \(\mathrm{Pe}=100\). Consistent with [3], this value balances convection and diffusion in \(\mathbf{Q}'\) without either dominating. Temperature fields are computed using a custom spectral CFD solver. Full details are given in [3], [12].

3.2 Finite-Time Density Structures↩︎

Figure 2 left shows the conventional Eulerian view of temperature evolution at a particular time slice, along with the stagnant case \(\widetilde{T}\) without convection. The right column shows the result of applying our visualization technique with \(t_0=0\), \(t_1=t_\text{slice}\), \(\tau = 0.05\), \(\sigma=0.04\), \(d=0.05\) and \(16{,}000\) particles. The forward integration shows attracting structures where particles concentrate (i.e., in yellow), indicating the regions of convective heating and the transport directions along which heat arrives there. The backward integration reveals the corresponding repelling structures, identifying regions of convective cooling. These structures are not visible in the temperature fields in Figure 2 nor in the fluid velocity field due to the non-trivial relationship between \(\boldsymbol{Q}'\) and \(\boldsymbol{u}\) (see Section 1).

3.3 Coherent Structures in Unsteady Transport↩︎

Beyond visualizing the full heat flux \(\boldsymbol{Q}'\), the method can be applied to individual flux contributions to gain targeted insight. Figure 3 shows the visualization applied to the diffusion contribution \(\alpha\boldsymbol{\nabla}T'\) of two flows with differing degrees of unsteadiness, with \(\tau\) increased to reveal coherency across a wider temporal slice.

The top case is driven by a steady flow, yet the resulting heat transport is unsteady. The largely non-crossing paths indicate coherent transport structure despite this unsteadiness. The bottom case is driven by an unsteady periodic flow, producing more complex aperiodic transport behavior. Despite the heat transport being aperiodic, the visualization reveals what appear to be periodic structures forming. At the blue arrow, transport rotates around an attracting region while moving upward, suggesting a periodic structure, most clearly visible in the provided supplementary video. Regions of consistent directionality are also visible, as seen in the blue rectangle, where uniform transport directions indicate that flow toward attracting structures varies little over the given timeframe.

4 Discussion & Conclusion↩︎

We presented a particle-based visualization technique that reveals finite-time density structures and coherent transport routes in unsteady divergent heat transport, taking a first step toward making coherent structures in unsteady thermal transport visible. Current limitations include the restriction to 2D flows and that large temporal windows \(\tau\) can produce clutter in highly unsteady flows. The method conceptually generalizes beyond the LTF to other divergent transport problems. Future work includes gaining deeper insight into the physical meaning of these structures, extending the method to three-dimensional heat transfer, and applying it to other transport datasets. Overall, the proposed method provides visual insight into convective heat transfer mechanisms that conventional methods cannot reveal.

References↩︎