Convective Heat Transfer in Uniformly Accelerated and Decelerated Turbulent Pipe Flows

Abstract

This study presents a detailed investigation of the temporal evolution of the Nusselt number ($Nu$) in uniformly accelerated and decelerated turbulent pipe flows under a constant heat flux using direct numerical simulations. The influence of different acceleration and deceleration rates on heat transfer is systematically studied, addressing a gap in the previous research. The simulations confirm several key experimental findings, including the presence of three distinct phases in the Nusselt number temporal response—delay, recovery, and quasi-steady phases—as well as the characteristics of thermal structures in unsteady pipe flow. In accelerated flows, the delay in the turbulence response to changes in velocity results in reduced heat transfer, with average $Nu$ values up to $48\%$ lower than those for steady-flow conditions at the same mean Reynolds number. Conversely, decelerated flows exhibit enhanced heat transfer, with average $Nu$ exceeding steady values by up to $42\%$ due to the onset of secondary instabilities that amplify turbulence. To characterize the $Nu$ response across the full range of acceleration and deceleration rates, a new model based on a hyperbolic tangent function is proposed, which provides a more accurate description of the heat transfer response than previous models. The results suggest the potential to design unsteady periodic cycles, combining slow acceleration and rapid deceleration, to enhance heat transfer compared to steady flows.

1. Introduction

Understanding the heat transfer characteristics of unsteady turbulent pipe flows is crucial for optimizing a wide range of industrial and technological processes. These include exhaust systems in internal combustion engines; heating, ventilation, and air conditioning (HVAC) systems; aerospace propulsion; chemical reactors; and biomedical technology applications. In recent decades, there has been growing interest in the potential of unsteady periodic pipe flows to enhance convective heat transfer compared to steady flows at the same Reynolds numbers. In particular, pipe flow driven by pulsations has attracted significant attention. Numerous experimental studies (Dec et al. [1], Habib et al. [2], Barker and Ffowcs Williams [3], Elshafei et al. [4], Patel and Attal [5], Simonetti et al. [6], Brahma and Singh [7]) and computational fluid dynamics (CFD) analyses (Wang and Zhang [8], Elshafei et al. [9], Nishandar et al. [10]) have investigated how pulsation parameters, such as amplitude, frequency, and mean Reynolds number, affect heat transfer. Despite extensive research, the effects of these parameters remain unclear, and contradictory findings continue to challenge the development of reliable predictive models.

Investigations have shown that pulsation can increase, decrease, or have no effect on heat transfer relative to steady conditions, depending on the pulsation control parameters. It is generally accepted that enhancing heat transfer requires the pulsation amplitude to exceed the mean velocity, leading to flow reversal (Dec et al. [1], Patel and Attal [5], Simonetti et al. [6], Brahma and Singh [7], Wang and Zhang [8]). However, it is still uncertain whether this condition alone is sufficient or if additional factors play a role in determining the heat transfer enhancement. When the pulsation amplitude is below the mean velocity, researchers have observed slight reductions (Habib et al. [2], Elshafei et al. [4,9]) in heat transfer, as well as cases where heat transfer remains similar to steady-flow conditions (Barker and Ffowcs Williams [3], Elshafei et al. [4,9], Nishandar et al. [10]).

Unlike amplitude, there is no consensus on how pulsation frequency affects heat transfer. Some studies indicate that heat transfer increases with higher pulsation frequencies (Dec et al. [1], Patel and Attal [5]), while others report that significant changes only occur at optimal or resonant frequencies (Habib et al. [2], Simonetti et al. [6], Wang and Zhang [8]). In contrast, numerous studies have found that varying the pulsation frequency has little to no significant effect (Barker and Ffowcs Williams [3], Elshafei et al. [4], Brahma and Singh [7], Elshafei et al. [9], Nishandar et al. [10]). There is currently no theory that explains the differing results reported in these experiments and simulations.

Similarly, the influence of the mean Reynolds number on heat transfer remains contentious. While most studies have not detected significant effects from variations in the Reynolds number (Barker and Ffowcs Williams [3], Elshafei et al. [4], Brahma and Singh [7], Elshafei et al. [9], Nishandar et al. [10]), some have observed that heat transfer enhancement due to pulsation increases with the mean Reynolds number up to a threshold, beyond which further increases yield no significant changes (Patel and Attal [5], Wang and Zhang [8]).

While the effects of pulsatile flows on heat transfer have been extensively studied, much less attention has been given to non-periodic unsteady flows, such as uniformly accelerated and decelerated flows. Understanding the heat transfer characteristics of these flows is important not only due to their broad industrial and civil engineering applications but also because unsteady periodic flows are composed of alternating acceleration and deceleration phases. Thus, investigating the heat transfer behavior during these individual phases could provide valuable insights into unresolved questions regarding heat transfer in pulsatile flows.

To date, only the laboratory experiments conducted by Shiibara et al. [11] and Nakamura et al. [12] have explored the effects of acceleration or deceleration on heat transfer in turbulent pipe flow. Using high-speed infrared thermography, these studies reported on the temporal response of the Nusselt number and thermal structures when flow speed changed abruptly (Shiibara et al. [11]). Based on these observations, a model was proposed to characterize this temporal response (Nakamura et al. [12]). Their results indicated a time delay in the heat transfer response to changes in velocity, resulting in lower heat transfer during acceleration and higher heat transfer during deceleration compared to steady conditions. However, the experimental setup did not allow control over the acceleration and deceleration rates (denoted by the parameter $\alpha$), so they were unable to explore the effects of varying $\alpha$.

To address this gap, the present work uses direct numerical simulations (DNS) to systematically investigate the impact of $\alpha$ on the temporal evolution of the Nusselt number in turbulent pipe flow under constant heat flux. DNS captures all relevant scales of motion and avoids the approximations of turbulence models, allowing for a detailed examination of flow structures and heat transfer mechanisms. However, due to its high computational cost, DNS has not been previously used to study heat transfer in unsteady flows.

This study contributes to the field in two significant ways. First, the simulations offer the first numerical verification of the experimental findings, confirming the existence of three phases in the Nusselt number evolution—delay, recovery, and quasi-steady phases—as well as the characteristics of the thermal structures during these phases. The results further show that, for the parameters investigated, accelerated flows may reduce heat transfer by up to $48\%$ compared to the steady conditions, while decelerated flows may increase it by up to $42\%$. These findings suggest that unsteady periodic flow cycles, composed of alternating slow acceleration and rapid deceleration phases, could be designed to achieve greater heat transfer efficiency compared to steady flows.

Second, it is found that the experimental model for the temporal response of the Nusselt number is accurate only for high values of $\left|\alpha \right|$. For moderate or small $\left|\alpha \right|$, the model fails to capture the response accurately. To address this limitation, a new model based on a hyperbolic tangent function is proposed, which accurately describes the temporal response of the Nusselt number across a broad range of $\alpha$ values.

2. Materials and Methods

2.1. Governing Equations and Dimensionless Parameters

The hydrodynamically and thermally developed flow of an incompressible fluid through a straight, horizontal pipe of a constant cross-sectional area is analyzed (Figure 1). The fluid properties, including density $\rho$, kinematic viscosity $\nu$, thermal conductivity $\kappa$, specific heat at constant pressure $C_p$, and thermal diffusivity ${\lambda }_t$, are assumed to be constant. The fluid motion is described by the continuity and momentum equations in cylindrical coordinates $(r,\theta ,z)$. The mean bulk velocity $v_{b,m}$ of the acceleration or deceleration period (definition is given below), the pipe diameter D, the dynamic pressure $\rho v_{b,m}^2$, and the advective time scale $D/v_{b,m}$ are used as characteristic scales for velocity, length, pressure, and time, respectively. These considerations lead to the following non-dimensional forms of the momentum and continuity equations:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial \mathbf{v}}{\partial t}}+(\mathbf{v}·\nabla )\mathbf{v}& =-\nabla p+{\displaystyle \frac{1}{R{e}_m}}{\nabla }^2\mathbf{v},\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill \nabla ·\mathbf{v}& =0,\hfill \end{array}$$

(2)

where $\mathbf{v}=(v_r,v_{\theta },v_z)$ represents the velocity field in cylindrical coordinates, t denotes time, p is the non-dimensional pressure, and

$$R{e}_m={\displaystyle \frac{v_{b,m}D}{\nu }},$$

(3)

is the mean Reynolds number.

Periodic boundary conditions are applied in the axial and azimuthal directions, whereas no-slip conditions at the pipe wall and regularity conditions at the axis are imposed in the radial direction.

This study examines flows that are uniformly accelerated or decelerated. To model this behavior, the non-dimensional bulk velocity $v_b\left(t\right)$ is updated at each time step according to the equation

$$v_{\mathrm{b}}(t+\delta t)=v_{\mathrm{b}}\left(t\right)+\alpha \delta t,$$

(4)

where $\delta t$ is the time step size, and $\alpha$ represents the non-dimensional rate of acceleration or deceleration, defined as

$$\alpha ={\displaystyle \frac{d{v}_b}{dt}},$$

(5)

The updated bulk velocity is enforced using the equation

$$v_{\mathrm{b}}(t+{\delta }_t)={\displaystyle \frac{1}{A}}\int v_{z}dA,$$

(6)

where A denotes the cross-sectional area of the pipe. The variation in the bulk velocity implies that the Reynolds number continuously changes during flow acceleration or deceleration. The characteristic Reynolds number in these cases is calculated using the mean bulk velocity, defined as $v_{b,m}={\displaystyle \frac{v_{b,init}+v_{b,end}}{2}}$, where $v_{b,init}$ and $v_{b,end}$ indicate the bulk velocity at the beginning and end of the acceleration or deceleration phase.

The external surface of the pipe is subjected to a constant and uniform heat flux per unit area, denoted by $q_w$. In this analysis, the effects of gravitational buoyancy are neglected, which allows us to treat the temperature T as a passive scalar. Viscous dissipation and internal heat generation are also neglected. Consequently, the temporal evolution of T is governed by the advection–diffusion equation, given in its dimensional form by

$${\displaystyle \frac{\partial T}{\partial t}}+\mathbf{v}·\nabla T=\kappa {\nabla }^{2}T.$$

(7)

Fluid properties are evaluated at the mean bulk temperature, denoted as $T_b$, which is defined as the volumetric flow rate weighted average temperature across the cross-sectional area of the pipe [13]:

$$T_b={\displaystyle \frac{{\int }_A{\langle v_z\rangle }^t{\langle T\rangle }^{t}dA}{{\int }_A{\langle v_z\rangle }^{t}dA}},$$

(8)

where ${\langle .\rangle }^t$ denotes time averaging. An energy balance conducted over a thin section of the pipe reveals a linear axial variation in $T_b$, characterized by the following rate of change:

$${\displaystyle \frac{d{T}_b}{dz}}={\displaystyle \frac{4q_w}{\rho C_p{v}_{b,m}D}}.$$

(9)

This scenario prevents the use of periodic boundary conditions in the axial direction. In line with previous studies [14,15,16], this limitation is addressed by replacing temperature variable T with the dimensionless temperature difference $\varphi ={\displaystyle \frac{T_b\left(z\right)-T(r,\theta ,z,t)}{T_{ref}}}$, where the reference temperature is defined as $T_{ref}={\displaystyle \frac{q_w}{\rho C_p{v}_{b,m}}}$. This modification eliminates the axial temperature gradient ${\displaystyle \frac{d{\langle \varphi \rangle }^t}{dz}}=0$, thus permitting the use of periodic boundary conditions.

The reformulated equation, expressed in non-dimensional terms using the previously defined scales, is

$${\displaystyle \frac{\partial \varphi }{\partial t}}+\mathbf{v}·\nabla \varphi -4v_z={\displaystyle \frac{1}{R{e}_{m}Pr}}{\nabla }^2\varphi ,$$

(10)

where $Pr$ is the Prandtl number, defined as

$$Pr={\displaystyle \frac{\nu }{{\lambda }_t}},$$

(11)

Periodic boundary conditions are used in the axial and azimuthal directions, with regularity conditions at the pipe axis and a prescribed heat flux, often termed the isoflux condition [14], applied at the pipe wall. The isoflux condition is expressed in non-dimensional form as

$${\displaystyle \frac{\partial \varphi }{\partial r}}{|}_{r=1/2}=-R{e}_{m}Pr.$$

(12)

An additional condition is required to ensure the uniqueness of the solution for Equation (10) [14]. It is imposed that the volumetric flow rate weighted average of $\varphi$ over the entire pipe volume must be zero:

$${\int }_V\varphi v_{z}dV=0.$$

(13)

To quantify convective heat transfer in accelerated or decelerated flows, the instantaneous Nusselt number is employed, which measures the relative importance between convective and conductive heat transport in the fluid:

$$Nu\left(t\right)={\displaystyle \frac{h\left(t\right)D}{\kappa }}={\displaystyle \frac{q_{w}D}{\kappa \langle {(T_w(t,\theta ,z)-T_b\left(z\right)\rangle }^{\theta ,z}}},$$

(14)

where $T_w$ is the wall temperature, and ${\langle .\rangle }^{\theta ,z}$ denotes averaging over the axial and azimuthal directions. With the non-dimensionalization carried out in this study, the instantaneous Nusselt number is expressed as

$$Nu\left(t\right)=-{\displaystyle \frac{R{e}_{m}Pr}{{\langle {\varphi }_w(t,\theta ,z)\rangle }^{\theta ,z}}}.$$

(15)

Table 1. Summary of the dimensionless numbers arising in the formulation of the problem.

Dimensionless Number	Symbol	Equation	Character	Value/s
Mean Reynolds number	$R{e}_m$	(3)	Input parameter	6400
Acceleration/Deceleration rate	$\alpha$	(5)	Input parameter	$0.0015\le \left|\alpha \right|\le 2$
Prandtl number	$Pr$	(11)	Input parameter	$0.7$
Nusselt number	$Nu$	(14)	Output parameter	-

provides a summary of the dimensionless numbers that arise in the mathematical formulation of the problem, distinguishing between the control parameters (input parameters) and those derived from the simulation results (output parameters). For the control parameters, the values used in this study are also presented.

2.2. Methodology and Numerical Code

The results presented in this paper were obtained from direct numerical simulations of the equations and boundary conditions outlined in Section 2.1. These simulations were performed using the open-source code NSPipeFlow [17], which has been recently extended to support non-isothermal flow simulations. The updated version of the code has been validated against previously published data in [14,16]. This enhanced version of the code is made publicly available in conjunction with this paper and is accessible in the repository [18].

The spatial discretization of the flow variables is performed using eighth-order central finite differences on a Gauss–Lobatto–Chebyshev grid in the radial direction, combined with Fourier–Galerkin expansions in the two homogeneous directions, $\theta$ and z. The Fourier–Galerkin expansions for the dimensionless velocity, pressure, and temperature fields are given by

$$\begin{array}{cc}\hfill \mathbf{v}(r,\theta ,z)& =\sum _{l=-L}^L\sum _{n=-N}^N\widehat{\mathbf{v}}(r,n,l)e^{i\left(n\theta +l{k}_{z}z\right)},\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill p(r,\theta ,z)& =\sum _{l=-L}^L\sum _{n=-N}^N\widehat{p}(r,n,l)e^{i\left(n\theta +l{k}_{z}z\right)},\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill \varphi (r,\theta ,z)& =\sum _{l=-L}^L\sum _{n=-N}^N\widehat{\varphi }(r,n,l)e^{i\left(n\theta +l{k}_{z}z\right)}.\hfill \end{array}$$

(18)

Here, n and l represent the Fourier mode numbers in the azimuthal and axial directions, respectively, and $k_z$ is the axial wavenumber, which defines the axial length of the computational domain as $L_z=2\pi /k_z$. The coefficients $\widehat{\mathbf{v}}(r,l,n)$, $\widehat{p}(r,l,n)$, and $\widehat{\varphi }(r,l,n)$ are complex spectral coefficients. The values of N and L specify the spectral resolution in the $\theta$ and z directions, respectively.

Substituting these expansions into Equations (1) and (10) yields a set of $(2N+1)\times (2L+1)$ independent equations, each associated with a specific $(n,l)$ pair. Solving this system enables the determination of the spectral coefficients. The time integration of these equations is carried out using a second-order accurate predictor–corrector scheme based on the Crank–Nicolson method [19]. To facilitate understanding of this algorithm, it is useful to define $N_{\widehat{\mathbf{v}}}=-(\widehat{\mathbf{v}}·\nabla )\widehat{\mathbf{v}}$ y $N_{\widehat{\varphi }}=-(\widehat{\mathbf{v}}·\nabla )\widehat{\varphi }+4\widehat{v_z}$ and to rewrite Equations (1) and (10) as

$$\begin{array}{cc}\hfill ({\displaystyle \frac{\partial }{{\partial }_t}}-{\displaystyle \frac{1}{R{e}_m}}{\nabla }^2)\widehat{\mathbf{v}}& =-\nabla \widehat{p}+N_{\widehat{\mathbf{v}}},\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill ({\displaystyle \frac{\partial }{{\partial }_t}}-{\displaystyle \frac{1}{PrR{e}_m}}{\nabla }^2)\widehat{\varphi }& =N_{\widehat{\varphi }}.\hfill \end{array}$$

(20)

Additionally, a pressure Poisson equation must be included, which is obtained by taking the divergence of Equation (19), along with the incompressibility condition

$${\nabla }^2\widehat{p}=\nabla ·N_{\widehat{\mathbf{v}}}.$$

(21)

In all these equations, the hat symbol indicates that the variables are spectral coefficients. The predictor step provides the initial estimates of the velocity and temperature fields at time step $q+1$ using data from time step q, and is formulated as

$$\begin{array}{cc}\hfill {\nabla }^2{\widehat{p}}_1^{q+1}& =\nabla ·N_{\widehat{\mathbf{v}}}^q,\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill ({\displaystyle \frac{1}{{\delta }_t}}-{\displaystyle \frac{c}{R{e}_m}}{\nabla }^2){\widehat{\mathbf{v}}}_1^{q+1}& =-\nabla {\widehat{p}}_1^{q+1}+N_{\widehat{\mathbf{v}}}^q+({\displaystyle \frac{1}{{\delta }_t}}-{\displaystyle \frac{(1-c)}{R{e}_m}}{\nabla }^2){\widehat{\mathbf{v}}}^q,\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill ({\displaystyle \frac{1}{{\delta }_t}}-{\displaystyle \frac{c}{PrR{e}_m}}{\nabla }^2){\widehat{\varphi }}_1^{q+1}& =N_{\widehat{\varphi }}^q+({\displaystyle \frac{1}{{\delta }_t}}-{\displaystyle \frac{(1-c)}{PrR{e}_m}}{\nabla }^2){\widehat{\varphi }}^q,\hfill \end{array}$$

(24)

where ${\delta }_t$ is the time step size, and c is a constant that sets the implicitness of the scheme (set to $0.5$ in our simulations). These initial estimates, ${\widehat{\mathbf{v}}}_1^{q+1}$ and ${\widehat{\varphi }}_1^{q+1}$, are then iteratively refined through a correction process. During each iteration, the non-linear terms are updated, and the velocity and temperature fields are refined by solving

$$\begin{array}{cc}\hfill {\nabla }^2{\widehat{p}}_{j+1}^{q+1}& =\nabla ·N_{{\widehat{\mathbf{v}}}_j}^{q+1},\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill ({\displaystyle \frac{1}{{\delta }_t}}-{\displaystyle \frac{c}{R{e}_m}}{\nabla }^2){\widehat{\mathbf{v}}}_{j+1}^{q+1}& =-\nabla {\widehat{p}}_{j+1}^{q+1}+c{N}_{{\widehat{\mathbf{v}}}_j}^{q+1}+(1-c)N_{\widehat{\mathbf{v}}}^q+({\displaystyle \frac{1}{{\delta }_t}}-{\displaystyle \frac{(1-c)}{R{e}_m}}{\nabla }^2){\widehat{\mathbf{v}}}^q,\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill ({\displaystyle \frac{1}{{\delta }_t}}-{\displaystyle \frac{c}{PrR{e}_m}}{\nabla }^2){\widehat{\varphi }}_{j+1}^{q+1}& =c{N}_{{\widehat{\varphi }}_j}^{q+1}+(1-c)N_{\widehat{\varphi }}^q+({\displaystyle \frac{1}{{\delta }_t}}-{\displaystyle \frac{(1-c)}{PrR{e}_m}}{\nabla }^2){\widehat{\varphi }}^q,\hfill \end{array}$$

(27)

where $j=1,2,\dots$

The iteration terminates once $\parallel {\widehat{\mathbf{v}}}_{j+1}^{q+1}-{\widehat{\mathbf{v}}}_j^{q+1}\parallel \le {10}^{-6}$ and $\parallel {\widehat{\varphi }}_{j+1}^{q+1}-{\widehat{\varphi }}_j^{q+1}\parallel \le {10}^{-6}$, with convergence typically achieved after one corrector iteration. The additional cost of evaluating advective terms twice per time step is offset by the larger $\delta t$ permitted by this scheme compared to other conventional methods.

To solve the pressure Poisson Equation (21), a homogeneous Neumann boundary condition, ${\displaystyle \frac{\partial \widehat{p}}{{\partial }_r}}=0$, is used at the pipe wall. Consequently, the velocity field obtained from (19) does not satisfy the divergence-free condition. This condition is subsequently enforced through influence matrices, a technique that ensures machine-level accuracy (typically of the order ${10}^{-16}$) and avoids the need for artificial pressure boundary conditions.

The code employs a hybrid MPI-OpenMP approach for parallelization, enabling efficient scaling across thousands of processors. For further details about the parallelization strategy and code functionalities, the reader is referred to [17] and the references therein.

All results presented in this article correspond to uniformly accelerated or decelerated flows, where the initial and final Reynolds numbers are $R{e}_{init}=3200$ and $R{e}_{end}=9600$, respectively, for accelerations, and vice versa for decelerations. The mean Reynolds number in all cases is therefore $R{e}_m=6400$. Within this range of Reynolds numbers, the flow is always turbulent under constant flow rate conditions. Figure 2 displays two examples of the temporal evolution of the instantaneous Reynolds number, $Re\left(t\right)={\displaystyle \frac{v_b\left(t\right)D}{\nu }}$. For accelerated flows (Figure 2a), the instantaneous Reynolds number increases linearly from $Re=3200$ to $Re=9600$. A steady phase is then simulated, where the flow rate remains unchanged to monitor the stabilization of flow characteristics at $Re=9600$. Similarly, for decelerated flows (Figure 2b), the instantaneous Reynolds number decreases linearly from $Re=9600$ to $Re=3200$, followed by a steady phase to examine stabilization at the lower Reynolds number. In all cases, $Pr$ is set to $0.7$. The initial conditions for these simulations are calculated from constant-flow rate simulations at $Re=3200$ and $Re=9600$.

To examine how acceleration or deceleration rates affect convective heat transfer, simulations were conducted with $\left|\alpha \right|$ values ranging from $0.0015$ to 2.

The simulations were carried out in a pipe of an axial length $3\pi D$, using 128 radial nodes, 121 azimuthal Fourier modes, and 400 axial Fourier modes. To confirm the resolution adequacy, additional simulations were performed using twice the number of radial nodes and $1.5$ times the Fourier modes in both the azimuthal and axial directions. These tests revealed no significant quantitative differences, validating the initial grid size. The time step size was set to ${\delta }_t=5\times {10}^{-3}$ for $\left|\alpha \right|<0.01$ and ${\delta }_t={10}^{-3}$ for $\left|\alpha \right|\ge 0.01$. All simulations were conducted on the Picasso supercomputer at the University of Málaga using 128 processors per simulation, resulting in durations ranging from 3 to 6 days depending on the value of $\alpha$.

3. Results

3.1. Temporal Variation in the Nusselt Number in Uniformly Accelerated Flows

This section examines the behavior of convective heat transfer in a uniformly accelerated flow. Figure 3a illustrates the temporal evolution of instantaneous Nusselt number $Nu$ (solid black line) at a moderate acceleration rate of $\alpha =0.02$. It also presents the corresponding quasi-steady Nusselt number values obtained from the empirical Gnielinski correlation (red dashed line). Three distinct phases can be identified in the temporal variation in $Nu$. In the initial phase (Phase 1), which lasts up to $t\approx 27$, $Nu$ remains nearly unchanged. This is followed by a growth phase (Phase 2), which is characterized by three sub-stages with different growth rates: a slow increase from $t\approx 27$ to $t\approx 35$, a more rapid increase from $t\approx 35$ to $t\approx 50$, and finally, a gradual approach to the quasi-steady value, which is reached around $t\approx 61$. In the last phase (Phase 3), $Nu$ fluctuates around the quasi-steady value.

The qualitative behavior of the $Nu$ response shown in Figure 3a is in agreement with experimental observations [11,12], even though those studies were conducted at higher acceleration rates. It is also consistent with the established descriptions of turbulence behavior in unsteady flows [20,21]. These descriptions propose that the response of turbulence to changes in the mean velocity occurs in three stages.

The initial stage, often referred to as the delay phase, corresponds to Phase 1 in Figure 3a. This phase is characterized by the phenomenon of frozen turbulence, where turbulent stresses respond with a delay to changes in mean velocity. Due to this delay, the intensity of the turbulent stresses remains largely unchanged despite the variation in the instantaneous $Re$. Consequently, the convective heat transfer carried by these stresses remains nearly constant, resulting in the almost invariant $Nu$ observed during this initial stage.

Once the turbulent fluctuations respond to the velocity change, the turbulence intensity increases (second stage), leading to the growth in $Nu$ observed in Phase 2. Eventually, the turbulence intensity stabilizes at the quasi-steady condition (third stage), similar to how $Nu$ stabilizes in Phase 3. These latter two stages are typically referred to as the recovery stage and quasi-steady stage, respectively.

A key result of the delayed response of convective heat transport to changes in the mean velocity is that, during acceleration, the instantaneous $Nu$ is always less than or equal to the quasi-steady $Nu$. This delay causes the average $Nu$ for a uniformly accelerated flow to be lower than that of a steady flow at the mean Reynolds number (in this study, $R{e}_m=6400$). For the specific case shown in Figure 3a, the average $Nu$ is $13.7$, while the $Nu$ for a steady flow at $Re=6400$ is $20.32$, indicating a $32.6\%$ reduction from the steady value. A detailed quantitative study on the reduction or increase in $Nu$ relative to the steady case, as a function of the $\alpha$ values, is provided in Section 3.4.

Comparing the temporal response of $Nu$ with that of the friction coefficient $c_f$ for the same case, shown in Figure 3, provides useful insights. The friction coefficient is calculated using the Fanning equation $c_f\left(t\right)={\displaystyle \frac{{\tau }_w\left(t\right)}{2v_b{\left(t\right)}^2}}$, where ${\tau }_w\left(t\right)$ denotes the instantaneous wall shear stress, and $v_b\left(t\right)$ is the instantaneous bulk velocity. The figure also includes quasi-steady values (red dashed line), which, for this range of $Re$ values, are well approximated by the empirical Blasius formula $c_f\left(t\right)=0.079Re{\left(t\right)}^{-0.25}$. The temporal variation in $c_f$ observed in the simulation qualitatively matches the previous findings [22], showing the three stages of the turbulence response in unsteady flows: delay, recovery, and quasi-steady phases. However, two important differences are noted between the temporal responses of $Nu$ and $c_f$.

The first difference occurs during the initial delay phase. Here, $c_f$ initially shows a slight increase above the quasi-steady value caused by the high inertia needed to start accelerating the flow. This is followed by a significant decrease, reaching values approximately $40\%$ lower than the quasi-steady case. This behavior contrasts with $Nu$, which remains constant during this phase. This observation deviates from the widely accepted Colburn analogy [23], which posits a proportional relationship between $c_f$ and $Nu$. While this analogy has been extensively validated in steady flows, the comparison in Figure 3 suggests that the relationship between these two quantities is more complex in unsteady flows.

The second difference is observed during the recovery phase. In this stage, $c_f$ rapidly increases due to the rise in turbulent stresses, surpassing the quasi-steady value and reaching a relative maximum at $t\approx 43$. Subsequently, $c_f$ decreases and stabilizes around the steady value at $t\approx 48$. In contrast, the growth rate of $Nu$ is much more gradual, resulting in a considerably longer recovery phase compared to $c_f$.

An important aspect investigated in the experiments of [11,12] is the spatiotemporal characteristics of heat transfer in unsteady flows. Specifically, they examined the evolution of the spatial distribution of the instantaneous convective heat transfer coefficient near the pipe wall, as the flow was accelerated or decelerated. To compare their observations with the simulation results, Figure 4 illustrates the spatial distribution of the temperature difference relative to the bulk temperature, $\varphi$ (left panels), near the pipe wall at various time instants during acceleration. Notably, this quantity is proportional to the convective heat transfer coefficient, and therefore, its spatiotemporal characteristics are analogous.

In the color scale used to depict $\varphi$ structures, dark blue represents regions of higher temperature, while yellow indicates areas of lower temperature within the section shown. Note that negative values of $\varphi$ indicate a temperature higher than the bulk temperature, which is typically observed near the wall. To aid in the interpretation of the physical processes driving the evolution of these thermal structures, the evolution of the radial velocity, $v_r$ (right panels), is also shown in Figure 4. The color maps for $v_r$ employ a blue-to-red scale, where negative values (blue) correspond to radial flow toward the pipe center, and positive values (red) correspond to radial flow toward the wall.

During the early delay phase, $\varphi$ structures appear as elongated streaks aligned with the flow direction (Figure 4a). These streaks show alternating regions of high and low temperature in the azimuthal direction and have slight modulations that appear to be associated with localized areas of significant radial velocity (Figure 4b). These areas are remnants of the initial steady turbulent state. As the flow accelerates, these regions do not regenerate and gradually dissipate. This characteristic is evident in Figure 4d, where significant radial velocity gradients are only observed in a small area near the outlet section of the pipe.

Without vortices to redistribute momentum and heat, the streaky structures of $\varphi$ lose the weak modulation observed in the early stages of the acceleration, becoming almost parallel and more elongated, often spanning the entire computational domain (Figure 4c). It is important to note that although some structures span the entire pipe length, suggesting that the axial domain used in the simulations may be insufficient to capture the full physics of the problem, additional simulations with extended pipe lengths were conducted and revealed no significant differences from the results obtained with the present configuration.

At the beginning of the recovery phase (point C in Figure 3a), regions of significant radial velocity begin to emerge across a large portion of the section (Figure 4f). The heat transport associated with these fluctuations causes the streaky structures to start oscillating and breaking down into smaller structures (Figure 4e). In areas where the radial velocity remains near zero, the $\varphi$ structures continue to form elongated streaks aligned with the flow direction.

As the recovery phase progresses, and the Nusselt number growth rate increases (point D in Figure 3a), radial velocity structures have spread almost entirely along the pipe’s axial length (Figure 4h). Consequently, the elongated streaks from earlier moments transform into shorter streaks that eventually break down into structures with varied spatial scales (see Figure 4g,i). Toward the end of the recovery phase, the $\varphi$ (Figure 4i) and $v_r$ (Figure 4j) structures become nearly identical to those observed in the steady case (Figure 4k and Figure 4l, respectively).

The magnitude of $\varphi$ substantially decreases during the recovery phase, indicating that as the flow becomes more turbulent, mixing is enhanced near the wall, leading to temperatures closer to the bulk temperature. A smaller value of $\varphi$ corresponds to a higher $Nu$, consistent with the evolution of this parameter shown in Figure 3a. Conversely, the magnitude of $v_r$ increases during the recovery phase, rising from very low values in the delay phase (often referred to as the relaminarization phase due to the low fluctuation intensity) to the typical levels seen in fully developed turbulent flow.

The evolution of the $\varphi$ structures depicted in this figure closely resembles the instantaneous convective heat transfer coefficient patterns reported by Nakamura et al. (see Figure 4 in [12]) for experiments in a similar Reynolds number range, further confirming the high fidelity of the simulations in reproducing the experimental results.

3.2. Temporal Variation in the Nusselt Number in Uniformly Decelerated Flows

This section examines the evolution of the Nusselt number in uniformly decelerated flows. Figure 5a illustrates the temporal response of $Nu$ for a simulation with a deceleration rate of $\alpha =-0.02$, matching the magnitude used for the uniformly accelerated case in Section 3.1. The three phases identified for uniformly accelerated flows are also present in decelerated flows.

Initially, there is a delay phase (Phase 1), during which $Nu$ remains approximately constant despite a decrease in instantaneous $Re$ (see inset in the figure). This phase is significantly shorter than in the uniformly accelerated case. Here, the delay phase extends to $t\approx 6$, accounting for $12\%$ of the deceleration period, whereas it covered nearly $50\%$ of the acceleration period in the accelerated flow.

Following the delay phase, the recovery phase (Phase 2) begins, during which $Nu$ decreases from its initial value to levels consistent with the final condition. Unlike the recovery phase in accelerated flows, where the intensity of turbulent fluctuations increases towards the quasi-steady level, in decelerated flows, the intensity of these fluctuations decreases to match the lower mean flow velocity. However, this decrease occurs more slowly than in the quasi-steady case (dashed red line), leading to greater convective heat transport and consequently higher $Nu$ compared to the quasi-steady case. This implies an increase in the intensity of turbulent fluctuations during the initial part of the recovery phase.

The exact cause of this increase is not entirely clear, but it may be associated with the presence of inflection points in the velocity profile, a characteristic feature of decelerated flows. These inflection points can induce instantaneous linear instabilities [24,25], potentially providing the energy needed for the transient increase in turbulence intensity.

The decrease in $Nu$ continues until $t\approx 110$, extending well beyond the duration of the deceleration period. This is followed by a slight increase, leading to the quasi-steady phase (Phase 3), during which $Nu$ oscillates slightly around a steady value. Notably, there is an offset between this steady value and the quasi-steady value predicted by the Gnielinski correlation. This offset is expected, as the correlation is known to deviate from the experimental values when $Re$ approaches the transitional regime. As in the uniformly accelerated case, the three phases observed in the Nu response to deceleration are consistent with the experimental observations in [11,12].

In contrast to the acceleration case, $Nu$ values during deceleration are always above or equal to those of the quasi-steady case, resulting in a net increase in heat transfer compared to the steady case when the flow is driven at the mean $Re$. For this specific case, the average $Nu$ during deceleration is $24.68$, while, as noted earlier, $Nu$ corresponding to the steady case for Re = 6400 is $20.32$. This leads to a net increase in $Nu$ of $21.45\%$.

A comparison of the temporal evolution of $Nu$ and $c_f$ throughout the deceleration phase (Figure 5b) reveals significant differences in their respective responses. The temporal response of $c_f$ can be divided into four distinct stages.

In the initial phase, $c_f$ undergoes a slight decrease due to the adverse pressure gradient applied to decelerate the flow. This phase is brief, lasting only until $t\approx 1$, which is six times shorter than the delay phase observed in the temporal response of $Nu$. Following this initial decrease, a second phase begins, characterized by a rapid recovery to values exceeding the quasi-steady level. During this phase, $c_f$ follows the same trend as the quasi-steady value but its magnitude remains slightly higher. This behavior supports the hypothesis of a local instability that transiently increases turbulence levels during flow deceleration. The onset of this second phase could therefore be linked to the initiation of this instability.

The third phase starts at the end of the deceleration period ($t=50$) and extends to $t\approx 140$. During this phase, $c_f$ exhibits a pronounced overshoot above the quasi-steady value due to the significant inertia of the fluid. After reaching this maximum, $c_f$ gradually decreases to values below the steady case before stabilizing.

In the fourth and final phase, $c_f$ oscillates around the steady value. Notably, as in the case of uniformly accelerated flow, $c_f$ reaches this final phase before $Nu$ stabilizes.

As discussed in the previous section, the qualitative differences in some phases of the $Nu$ and $c_f$ temporal responses suggest that the physical mechanisms governing the instantaneous values of these parameters in unsteady flows may differ. This observation raises questions about the applicability of the Colburn analogy for unsteady-flow conditions.

The evolution of the $\varphi$ structures near the wall during uniform deceleration (left panels of Figure 6) shows significant differences compared to the case of uniform acceleration depicted in Figure 4. As in Figure 4, the evolution of $v_r$ is shown alongside the $\varphi$ structures, using the same color palettes.

During the delay phase (illustrated in Figure 6a,b for $t\approx 5.5$), both the $\varphi$ and $v_r$ structures remain similar to the initial steady turbulent state. The thermal structures are primarily organized into streaks, alternating regions of high (blue) and low (yellow) temperatures in the azimuthal direction, which coexist with smaller structures.

Comparing the $\varphi$ structures with the $v_r$ distribution reveals that regions with smaller structures align with areas where the radial velocity is more pronounced. These regions are identified in Figure 6b as spatially localized regions with closely spaced high positive (dark red) and negative (dark blue) radial velocities in the azimuthal direction.

A significant change observed during the transition to the recovery phase is the elongation of the $v_r$ structures (Figure 6d). Regarding the $\varphi$ structures, fewer small structures are observed (Figure 6c), and the streaky structures exhibit several clear differences compared to the previous phase: a marked increase in both the axial length and azimuthal width of the structures, and the onset of a certain undulation. This undulation is consistent with the emergence of a secondary instability as previously speculated, which enhances turbulence levels and causes $Nu$ to rise above the quasi-steady value.

As time progresses and the deceleration period nears its end, the width of the $\varphi$ structures continues to grow, and their undulation becomes more pronounced due to fluctuations extracting energy from the secondary instability (Figure 6e). An increase in the magnitude of $\varphi$ is also evident, consistent with the decrease in $Nu$ that results from the diminishing intensity of turbulent fluctuations as the instantaneous $Re$ decreases. This substantial reduction in turbulent fluctuation intensity is clearly visible in Figure 6f. Additionally, this figure shows that the distribution of $v_r$ remains similar to that at the start of the recovery phase, though the azimuthal length of the structures has significantly increased.

The transition between the deceleration phase and the subsequent steady phase is marked by a clear change in the topology of the structures (Figure 6g,h). The elongated streaks observed in earlier stages are replaced by irregular structures with a large azimuthal length, similar to the “mottled structure” observed experimentally in [11,12]. This change likely results from streak collapse due to the local secondary instability during the recovery phase.

Notably, up to this point, the $\varphi$ structures consistently displays negative values near the wall, indicating higher temperatures than the bulk temperature. However, Figure 6g shows that some thermal structures now have positive values, indicating temperatures lower than the average. This change indicates significant heat and momentum transport from the central part of the pipe, where the temperature is lower, to the wall, which is also consistent with the substantial increase in $c_f$ relative to the steady value observed in Figure 5b during the third phase of the temporal response of this parameter. This increase in $c_f$ is also consistent with the substantial rise in the magnitude of $v_r$ observed in Figure 6h.

As $Nu$ approaches the end of the recovery phase, the turbulent fluctuations arising from the secondary instability gradually dissipate, and the $\varphi$ structures revert to streaks aligned with the flow direction but with a significantly larger azimuthal length than during the initial deceleration stages (Figure 6i). This feature is again consistent with the experimental observations in [11,12]. The evolution of the structures during this stage is similar to what occurs during the frozen turbulence phase when the flow accelerates. Initially, the intensity of turbulent fluctuations remains at very low levels for some time (Figure 6j), giving rise to $\varphi$ streaks that extend axially across the entire computational domain, along with smaller streaks that emerge from the breakup of larger streaks in regions where $v_r$ is significant. Eventually, the intensity of turbulent fluctuations increases to adapt to the final steady state (Figure 6l), and the thermal structures take on the characteristic distribution of a steady turbulent flow, exhibiting a pattern of streaks of various sizes alternating high and low temperatures in the azimuthal direction (Figure 6k).

3.3. Characterization of the Temporal Variation in the Nusselt Number as a Function of Acceleration or Deceleration Rate

This section examines the temporal characteristics of the $Nu$ variation as a function of the acceleration or deceleration rate $\alpha$ and introduces a simple model that satisfactorily reproduces the $Nu$ response across a wide range of $\alpha$ values.

Figure 7 presents the temporal evolution of $Nu$ for uniformly accelerated flows, covering $\alpha$ values spanning three orders of magnitude. Two distinct behaviors emerge depending on the magnitude of $\alpha$.

For very small values ($\alpha \le 0.002$, Figure 7a), $Nu$ increases quasi-steadily over time. The change in the mean flow velocity is slow enough for the turbulent fluctuations to adjust almost instantaneously to the evolving flow conditions. As a result, the $Nu$ value corresponding to each instantaneous $Re$ value closely matches that of a steady flow at the same $Re$. In these cases, the Gnielinski correlation (indicated by the red dashed lines in the figure) provides a good estimate of the $Nu$ evolution. Initially, the correlation slightly underestimates $Nu$ because the instantaneous $Re$ is close to transitional values, where this semi-empirical correlation is known to be less accurate. However, as time progresses and the instantaneous $Re$ moves further from the transitional regime, $Nu$ converges with high precision to the value predicted by the Gnielinski correlation.

For $\alpha >0.002$ (Figure 7b), the evolution of $Nu$ follows the three phases described in Section 3.1. The delay phase (demarcated by the brown dashed line) lasts until $t\approx 27$ in all cases, demonstrating that its duration is independent of $\alpha$. However, the range of $Re$ values encompassed during this phase expands as $\alpha$ increases due to the more rapid change in mean velocity, which results in a higher $Re$ by the end of the delay phase. When $\alpha$ exceeds $0.02$, the transition time between the initial and final $Re$ values becomes shorter than the delay phase, causing this phase to extend beyond the acceleration period. Despite this, the qualitative behavior of the $Nu$ evolution is consistent with that observed at lower $\alpha$ values, with the notable exception that the increase in $Nu$ and its approach to quasi-steady values occur while the flow is already being driven at a constant $Re$.

The most significant effect of increasing $\alpha$ is the faster growth in $Nu$ during the recovery phase. As $\alpha$ increases, this increase becomes increasingly sharper until $\alpha \approx 0.1$. Beyond this point, further increases in $\alpha$ have minimal impact on the recovery phase, and the temporal evolution of $Nu$ is practically identical in all cases as observed for $\alpha =0.1$ and $\alpha =0.2$.

The dependence of the $Nu$ temporal response with $\alpha$ in uniformly decelerated flows, illustrated in Figure 8, reveals two significant differences compared to the uniformly accelerated case. The first is that, for low deceleration rates, no quasi-steady variation in $Nu$ is observed. Even at the lowest deceleration rate considered ($\alpha =-0.0015$), shown in Figure 8a, the three phases described in Section 3.2 are still present. After the delay phase, which extends until $t\approx 10$ (see the inset in the figure), $Nu$ decreases to values close to those predicted by the Gnielinski correlation but progressively deviates as time advances. This deviation can be attributed to two factors. First is the presence of a secondary instability that increases turbulence levels and convective heat transfer beyond what would exist in a quasi-steady state. Second, as the instantaneous $Re$ approaches transitional values, the Gnielinski correlation becomes less accurate. The first factor explains the early-stage deviations, while the second factor accounts for deviations later in the deceleration and during the quasi-steady phase that follows.

The second notable difference is the dependence between the delay phase duration and $\alpha$. As the magnitude of $\alpha$ increases, the delay phase shortens (see inset of Figure 8b). However, this trend does not hold across the entire range of $\alpha$ values studied. For rapid decelerations as shown in Figure 8c, the delay phase stabilizes at $t\approx 5$. The variation in the delay phase duration observed for $0.0015\le \left|\alpha \right|\le 0.02$ may be linked to the onset of the secondary instability. As the magnitude of $\alpha$ increases, the instability sets in earlier, causing variations in the turbulent fluctuations level and the associated convective transport. As a result, $Nu$ deviates from its initial value earlier. However, for rapid decelerations, the transition between the initial and final states occurs so quickly that the onset of instability is similar regardless of the value of $\alpha$. This could explain why the duration of the delay phase becomes independent of $\alpha$ at higher deceleration rates.

Similar to accelerated flows, the $Nu$ values during the recovery phase vary more sharply as the deceleration rate increases (Figure 8b), until reaching a limit at $\alpha =-0.04$. Beyond this threshold, further increases in $\alpha$ do not significantly affect the temporal response of $Nu$ (Figure 8c).

The evolution of thermal structures during both acceleration and deceleration is qualitatively similar to the structures described in Section 3.1 and Section 3.2. The exception is the quasi-steady cases for uniformly accelerated flows, where the thermal structures exhibit the characteristic pattern of turbulent flow: streaks of varying sizes coexisting with smaller structures that become finer as the instantaneous $Re$ increases.

Nakamura et al. [12] propose a model to characterize the temporal variation in $Nu$ in unsteady flows upon sudden acceleration (deceleration), assuming an exponential growth (decay) of $Nu$ during the recovery phase. This model is based on two parameters: the delay phase duration $\Delta t_d$ and a parameter $\tau$, which controls the steepness of the exponential growth (decay). According to this model, the instantaneous $Nu$ is given by the following piecewise function:

$$Nu\left(t\right)=\left\{\begin{array}{cc}N{u}_{init}\hfill & \mathrm{if}t\le \Delta t_d\hfill \\ N{u}_{end}-(N{u}_{end}-N{u}_{init})exp\left(-{\displaystyle \frac{t-\Delta t_d}{\tau }}\right)\hfill & \mathrm{if}t>\Delta t_d,\hfill \end{array}\right.$$

(28)

where $N{u}_{init}$ and $N{u}_{end}$ represent the $Nu$ values corresponding to the initial and final $Re$ values in steady-flow conditions. When applied to the simulation data for uniformly accelerated flows, this model accurately reproduces the $Nu$ response for high acceleration rates (Figure 9a). However, for moderate $\alpha$ values (Figure 9b) and quasi-steady cases (Figure 9c), the model fails to adequately predict the $Nu$ growth.

To address this limitation, a new model is proposed, where the $Nu$ growth is modeled using a hyperbolic tangent function. In this model, the temporal variation in $Nu$ is expressed as

$$Nu\left(t\right)=0.5\left(N{u}_{init}+N{u}_{end}\right)+0.5\left(N{u}_{end}-N{u}_{init}\right)tanh\left({\displaystyle \frac{t-t_{inf}}{s}}\right).$$

(29)

This model also introduces two parameters: $t_{inf}$, which marks the inflection point of the $Nu$ growth curve, and s, which controls the steepness of the curve, similar to $\tau$ in the previous model. As shown in the lower panels of Figure 9, the proposed model accurately estimates the $Nu$ response across the entire range of $\alpha$ values. For high $\alpha$ values, the new model matches the accuracy of the exponential model (Figure 9d). However, unlike the exponential model, it also accurately predicts the $Nu$ variation for moderate $\alpha$ values (Figure 9e).

Even for low $\alpha$ values, where $Nu$ varies quasi-steadily, the proposed model provides a reasonable estimate, with only a slight overestimation during the initial phase. In quasi-steady cases, the Gnielinski correlation (shown as a green dashed line in Figure 9e) remains the most accurate predictor of $Nu$ values.

The variation in the model parameters $t_{inf}$ and s with $\alpha$ is presented in Figure 10. The parameter $t_{inf}$, shown in Figure 10a, decreases with increasing $\alpha$ and eventually stabilizes at an approximately constant value for $\alpha >1$. This behavior is well described by an exponential function with three parameters (blue line in the figure):

$$t_{inf}=0.2{\alpha }^{-1.1}+29$$

(30)

The variation in s, displayed in Figure 10b, exhibits two distinct phases. For low to moderate values of $\alpha$ (up to $\alpha \approx 0.02$), s decreases sharply as $\alpha$ increases. However, beyond this threshold ($\alpha >0.02$), the decrease becomes much more gradual. The variation in s across the entire range of $\alpha$ is well approximated by the following function (again represented by the blue line in the figure):

$$s=0.5(0.33{\alpha }^{-1.04}+0.175{\alpha }^{-0.54}-9.76)+0.5(0.175{\alpha }^{-0.54}-0.33{\alpha }^{-1.04}+13.96)tanh\left({\displaystyle \frac{\alpha -t_{inf}}{s}}\right)$$

(31)

Similar conclusions can be drawn when these models are applied to characterize uniformly decelerated flow. The exponential model reasonably estimates the temporal response of $Nu$ for high deceleration rates, but its accuracy decreases as the deceleration rate lowers. This is evident in the upper panels of Figure 11. In Figure 11a, which shows a high deceleration rate simulation, the model satisfactorily reproduces the instantaneous $Nu$ values, except at the end of the recovery phase, where it overestimates the simulation results, and during the initial delay phase, where it predicts an average value and therefore fails to capture the oscillations observed during this phase.

For moderate decelerations (Figure 11b), the model not only overestimates $Nu$ before the quasi-steady phase but also shows slight deviations during the early recovery phase. These deviations increase as the absolute value of $\alpha$ decreases as shown in Figure 11c.

Similar to the uniformly accelerated case, the hyperbolic tangent-based model proposed here satisfactorily estimates the temporal evolution of $Nu$ across the entire range of $\alpha$ values. For high deceleration rates (Figure 11d), it is slightly less accurate than the exponential model, overestimating $Nu$ during the final part of the recovery phase and slightly underestimating the average value during the delay phase. This underestimation persists for all $\alpha$ values. However, as the absolute value of $\alpha$ decreases, the proposed model captures the recovery phase much more accurately than the exponential model, significantly reducing the overestimation of the $Nu$ values during the approach to the quasi-steady phase (Figure 11e,f).

The variation in the parameters $t_{inf}$ and s with $\left|\alpha \right|$ follows a trend similar to that observed in uniformly accelerated flows. For $t_{inf}$ (Figure 12a), a decrease is observed as $\left|\alpha \right|$ increases, which can be accurately fitted with a three-parameter exponential function:

$$t_{inf}=0.51{\left|\alpha \right|}^{-1}+14.91.$$

(32)

For parameter s (Figure 12b), two distinct phases are observed: a sharp decrease up to $\left|\alpha \right|\approx 0.15$, followed by a phase where s remains approximately constant, around $11.16$. This behavior is well captured by the following expression:

$$s={0.5\left(0.41\right|\alpha |}^{-1}{+19.67)+0.5(2.65-0.41\left|\alpha \right|}^{-1})tanh\left({\displaystyle \frac{\left|\alpha \right|-t_{inf}}{s}}\right).$$

(33)

3.4. Heat Transfer Gain or Loss Relative to Steady Flow

This section examines the heat transfer gain or loss in unsteady flows compared to steady flows with the same mean $Re$. As discussed in the introduction, unsteady flows can potentially enhance heat transfer in industrial processes compared to maintaining a constant flow rate. To investigate this possibility, it is essential to quantify how variations in the parameters governing unsteady flows affect $Nu$. This study specifically examines the impact of the acceleration (or deceleration) rate on $Nu$.

To quantify the heat transfer gain or loss compared to the steady case, the following parameter is defined:

$$\Delta Nu\left(\alpha \right)={\displaystyle \frac{N{u}_m\left(\alpha \right)-N{u}_{steady}}{N{u}_{steady}}}\times 100,$$

(34)

where $N{u}_m\left(\alpha \right)$ represents the average Nusselt number during the period of acceleration or deceleration (excluding the steady period that follows these phases in the simulations), $N{u}_{steady}$ is the Nusselt number corresponding to the mean Reynolds number ($R{e}_m=6400$) obtained from the Gnielinski correlation, which provides an accurate estimate of $Nu$ in steady conditions for this $Re$. The values of $\Delta Nu$ as a function of $\alpha$ for uniformly accelerated and decelerated flows are shown in Figure 13a,b, respectively.

As explained in Section 3.1, the significant delay in the turbulence response to changes in the mean velocity results in a reduced heat transfer rate for uniformly accelerated flows compared to steady flows. It is worth noting that the average $Nu$ in cases with quasi-steady behavior is slightly higher (by approximately $3\%$) than the steady $Nu$ value. This deviation is, however, consistent with the expected errors in estimating this parameter using the Gnielinski correlation at $Re=6400$, suggesting that heat transfer remains effectively unchanged from the steady case for these $\alpha$ values. For moderate $\alpha$ values, however, small changes in $\alpha$ can cause significant variations in the average $Nu$. The largest losses, approximately $\Delta Nu=-48\%$ relative to the steady flow, occur for $\alpha$ values, where the delay phase extends beyond the acceleration period.

In contrast, as discussed in Section 3.2, $Nu$ values during deceleration remain consistently above the quasi-steady values. This is partly due to the presence of an instability that increases the intensity of turbulent fluctuations, resulting in a higher average $Nu$ compared to the steady case. Even at the smallest values of $\left|\alpha \right|$, the average $Nu$ shows an approximate 7.5% gain over steady conditions. However, the increase in $Nu$ during deceleration is more gradual than the sharp decrease observed for accelerated flows, reaching maximum values of $\Delta Nu\approx 42\%$ in simulations with large $\left|\alpha \right|$, where the delay phase extends beyond the deceleration period.

These results suggest the potential for designing a periodic unsteady-flow cycle that includes a slow acceleration phase to minimize heat transfer losses, followed by a rapid deceleration phase to significantly enhance heat transfer compared to the steady case. The average Nusselt number for such a cycle would be higher than that of a steady flow, while maintaining the same mean Reynolds number in both cases.

4. Discussion

This paper presents the first study to use direct numerical simulations to systematically analyze the temporal evolution of the Nusselt number ($Nu$) in uniformly accelerated and decelerated turbulent pipe flows under constant heat flux. The findings confirm previous experimental results and provide new insights into how acceleration and deceleration rates affect heat transfer, a parameter not systematically explored in earlier research.

The results confirm the existence of three distinct phases, the delay, recovery, and quasi-steady phases, consistent with the previous experimental observations [11,12]. These phases are present at all acceleration rates ($\alpha$), except when $\alpha$ is very low for accelerated flows, in which case $Nu$ evolves quasi-steadily over time.

In accelerated flows, the delayed turbulence response reduces heat transfer by up to $48\%$ compared with the steady case, while decelerated flows exhibit a maximum 42% increase due to a secondary instability that intensifies turbulence. Evidence of this instability is provided by the observed increase in turbulence intensity and convective transport, as well as changes in the topology of the flow structures. The typical elongated streaks seen in turbulent flows are replaced by more irregular structures with significant azimuthal extent, akin to the ’mottle structure’ reported in [12]. While the exact mechanism driving the instability is unclear, it may be linked to inflection points in the velocity profile as suggested by recent studies [24,25]. Further research is needed to explore this phenomenon in detail.

Another important observation is the distinct temporal behavior of the friction coefficient ($c_f$) compared to $Nu$. It is observed that $c_f$ responds more quickly to velocity changes, challenging the applicability of the Colburn analogy, which is commonly used in steady flows. This suggests that separate models are needed to accurately describe the dynamics of friction and heat transfer in unsteady turbulent flows.

This analysis also introduces a new model based on a hyperbolic tangent function that accurately characterizes the $Nu$ response across various acceleration and deceleration rates. This model is more versatile than previous ones and may serve as a valuable tool in advancing both research and industrial practices related to optimizing heat transfer in unsteady-flow systems. However, it is valid only for the specific conditions investigated, as other parameters such as the initial and final Reynolds numbers, or non-linear acceleration profiles have not been considered. Future work will focus on extending the model to broader conditions.

One potential application of these findings is the design of flow cycles that enhance heat transfer in industrial processes. The results suggest that a periodic unsteady-flow cycle, comprising a slow acceleration phase to minimize heat transfer losses followed by a rapid deceleration phase to boost heat transfer, could outperform steady-state scenarios. Similar cycles have been recently proposed in the literature to reduce friction losses in turbulent pipes while producing net energy savings [26]. However, further research is required to investigate the behavior of $Nu$ during direct transitions between acceleration and deceleration phases, as these may introduce new dynamics not captured in the present study.