%20-%20https://sketch.com%20--%3e%3ctitle%3e3ECBBF50-FD08-412A-943D-2CFC5E025A16%3c/title%3e%3cdesc%3eCreated%20with%20sketchtool.%3c/desc%3e%3cg%20id='Stylesheet'%20stroke='none'%20stroke-width='1'%20fill='none'%20fill-rule='evenodd'%3e%3cg%20transform='translate(-926.000000,%20-668.000000)'%20id='BUTTONS'%3e%3cg%20transform='translate(100.000000,%20597.000000)'%3e%3cg%20id='BTN-/-RIGHTS'%20transform='translate(816.000000,%2065.000000)'%3e%3cg%20id='ICON-/-Rights'%20transform='translate(10.000000,%206.000000)'%3e%3cpath%20d='M13.5752841,0%20L1.34943182,0%20C0.59375,0%200,0.586666667%200,1.33333333%20L0,22.6666667%20C0,23.4133333%200.620738636,24%201.34943182,24%20L17.6505682,24%20C18.40625,24%2019,23.3866667%2019,22.6666667%20L19,5.36%20L13.5752841,0%20Z%20M13.5752841,1.89333333%20L17.0838068,5.36%20L13.5752841,5.36%20L13.5752841,1.89333333%20Z%20M1.34943182,22.6666667%20L1.34943182,1.33333333%20L12.1988636,1.33333333%20L12.1988636,5.36%20C12.1988636,6.10666667%2012.7926136,6.69333333%2013.5482955,6.69333333%20L17.6235795,6.69333333%20L17.6505682,22.6666667%20L1.34943182,22.6666667%20Z'%20id='Shape'%20fill='%23595959'%20fill-rule='nonzero'%3e%3c/path%3e%3ccircle%20id='Oval'%20stroke='%23595959'%20stroke-width='1.5'%20fill='%23595959'%20cx='9.5'%20cy='13.5355339'%20r='5.5'%3e%3c/circle%3e%3cpolyline%20id='Path-2'%20stroke='%23FFFFFF'%20stroke-width='1.5'%20transform='translate(9.500000,%2013.035534)%20rotate(45.000000)%20translate(-9.500000,%20-13.035534)%20'%20points='8%2016.0355339%2011%2016.0355339%2011%2010.0355339'%3e%3c/polyline%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/g%3e%3c/svg%3e)
1. Introduction
Vertical heated pipe flow is widely used in engineering applications, e.g. in geothermal energy capture, nuclear reactor cooling systems and fossil-fuel power plants, to transfer heat from one device to another. The important difference from iso-thermal pipe flow is that buoyancy, caused by the expansion of the fluid near the heated wall, can partially or even fully drive the flow, referred to as mixed or natural convection. Mixed convection has been widely researched, due to the interesting and significant effects of buoyancy on the dynamics of the flow and on heat-transfer performance. Buoyancy plays a different role in downward versus upward flow. In the former, buoyancy acts in the opposite direction to the flow and always enhances heat transfer. In an upward flow, buoyancy first deteriorates the heat transfer (Ackerman Reference Ackerman1970), then heat transfer recovers only when buoyancy is strong enough. Three typical regimes for the heat-transfer characteristics are identified in upward heated flow, namely shear turbulence, laminarised flow and convective turbulence (Parlatan, N.E. & M.J. Reference Parlatan, Todreas and Driscoll1996; Yoo Reference Yoo2013; Zhang et al. Reference Zhang, Xu, Liu and Dang2020). (See also figure 1 a.)
Extensive research has been conducted to understand the mechanism of heat-transfer deterioration in upward heated flow (Hall & Jackson Reference Hall and Jackson1969; Steiner Reference Steiner1971; Carr, Connor & Buhr Reference Carr, Connor and Buhr1973; Polyakov & Shindin Reference Polyakov and Shindin1988; Satake et al. Reference Satake, Kunugi, Mohsen Shehata and McEligot2000; You, Yoo & Choi Reference You, Yoo and Choi2003; Bae et al. Reference Bae, Yoo, Choi and McEligot2006; Jackson Reference Jackson2013; He, He & Seddighi Reference He, He and Seddighi2016). Hall & Jackson (Reference Hall and Jackson1969) and Jackson & Hall (Reference Jackson and Hall1979) proposed that the reduced shear stress in the buffer layer caused by buoyancy leads to a reduction of turbulence production, suppressing turbulence and even laminarising the flow, consequently deteriorating the heat transfer. More recently, He et al. (Reference He, He and Seddighi2016) successfully reproduced the laminarisation phenomenon by modelling the buoyancy with a radially dependent body force added to the isothermal flow. They noticed that the body force causes little difference to the key characteristics of turbulence (in particular, the turbulent viscosity), and proposed that laminarisation is caused by the reduction of the ‘apparent Reynolds number’, which is calculated based only on the pressure force of the flow (i.e. excluding the contribution from the body force). A similar laminarisation phenomenon is also found in isothermal pipe flow, where it has attracted much attention due to its implications for drag reduction (Hof et al. Reference Hof, de Lozar, Avila, Tu and Schneider2010; He et al. Reference He, He and Seddighi2016; Kühnen et al. Reference Kühnen, Song, Scarselli, Budanur, Riedl, Willis, Avila and Hof2018; Marensi, Willis & Kerswell Reference Marensi, Willis and Kerswell2019). Kühnen et al. (Reference Kühnen, Song, Scarselli, Budanur, Riedl, Willis, Avila and Hof2018) examined the phenomenon of laminarisation from the perspective of the self-sustaining process of shear turbulence (Hamilton, Kim & Waleffe Reference Hamilton, Kim and Waleffe1995) and suggested that the decay of turbulence can be triggered by modifying the flow such as to reduce transient growth. Marensi et al. (Reference Marensi, Willis and Kerswell2019) investigated this phenomenon using nonlinear stability analysis (Pringle & Kerswell Reference Pringle and Kerswell2010), and found that nonlinear stability is enhanced in the presence of a body force that flattens the velocity profile. Recently, for the vertical heated pipe, Marensi, He & Willis (Reference Marensi, He and Willis2021) systematically studied the flow regimes and found evidence that heat-transfer deterioration and laminarisation are caused by weakened streamwise vortices.
Enhanced heat transfer means more effective and efficient energy conversion or cooling, and thus there have been many interesting investigations aimed at improving the heat transfer in fluid systems. Strategies can generally be classified as active, passive and compound remedies (Webb & Bergies Reference Webb and Bergies1983; Liu & Sakr Reference Liu and Sakr2013; Kumar & Kim Reference Kumar and Kim2015; Suri, Kumar & Maithani Reference Suri, Kumar and Maithani2018). Active methods (Ohadi, Nelson & Zia Reference Ohadi, Nelson and Zia1991; Wang, Zhou & Sun Reference Wang, Zhou and Sun2020; Yuan et al. Reference Yuan, Zou, Yang and Chen2023) require an external power input to improve heat transfer. For example, Ohadi et al. (Reference Ohadi, Nelson and Zia1991) studied the effect of corona discharge on forced-convection heat transfer in a tube. Wang et al. (Reference Wang, Zhou and Sun2020) and Yuan et al. (Reference Yuan, Zou, Yang and Chen2023) proposed a method of vibrating the boundary layer to enhance the heat transfer. Passive methods include a curving or twisting flow geometry, adding extended surfaces and so on. Compound methods (Gau & Lee Reference Gau and Lee1992; Naphon et al. Reference Naphon, Wiriyasart, Arisariyawong and Nualboonrueng2017; Kareem & Gao Reference Kareem and Gao2018) adopt both active and passive techniques. These strategies have significantly improved heat transfer in many systems.
Many techniques for heat-transfer enhancement have been developed empirically. Here, we seek a mathematical strategy that is ‘optimal’ with respect to a constrained magnitude of an applied body force. In principle, maximisation of the heat transfer can be solved by variational methods. However, the heat transfer, measured by the Nusselt number
$\textit {Nu}$
, is a local field variable, depending on the gradient of the temperature evaluated at the wall, which leads to awkward delta functions in a variational approach. A quantity is required for variational formulations that is spatially global, but which measures the heat transfer at the wall. A suitable quantity is the scalar dissipation, which has a history stemming back to its use in the study of mixing of passive scalars (Batchelor Reference Batchelor1959). In the thermal community, this has recently become known as the entransy dissipation (Guo, Zhu & Liang Reference Guo, Zhu and Liang2007; Kostic Reference Kostic2017), where it was derived via an analogy between heat conduction and electrical conduction. This quantity has been successfully used to optimise heat transfer in various thermal systems, e.g. in heat exchangers (Guo et al. Reference Guo, Liu, Tao and Shah2010; Guo & Xu Reference Guo and Xu2012) and heat exchanger networks (Chen et al. Reference Chen, Chen, Li and Guo2009).
For pipe flow, Meng, Liang & Li (Reference Meng, Liang and Li2005) have sought a steady velocity field that maximises heat transfer. Although the Navier–Stokes equation was not prescribed as a constraint, it was shown that the velocity field must satisfy a similar equation, subject to a particular force called the synergy force, which produces a velocity field that tends to align with the temperature gradient (Guo Reference Guo2001). Jia et al. (Reference Jia, Liu, Liu and Nakayama2014) did a similar optimisation but set power consumption as a constraint condition. They also found that longitudinal swirl flow with multi-vortex structure can enhance heat transfer greatly, and the number of vortices of the optimal velocity field increases for a larger power consumption. Wang, Liu & Liu (Reference Wang, Liu and Liu2015) proposed a similar criterion for heat-transfer optimisation, exergy destruction minimisation, and a similar optimal velocity field was found. Such heat-transfer optimisations in pipe flow have motivated several heat-transfer enhancement designs, e.g. the alternating elliptical axis tube (Meng et al. Reference Meng, Liang and Li2005), discrete double-inclined ribs tubes (Li, Meng & Guo Reference Li, Meng and Guo2009) and many other interesting modifications (Liu & Sakr Reference Liu and Sakr2013; Sheikholeslami, Gorji-Bandpy & Ganji Reference Sheikholeslami, Gorji-Bandpy and Ganji2015). However, the above calculations have assumed a steady laminar flow, while it is common to find that the flow is turbulent. There have been efforts to construct variational equations based on the Reynolds-averaged Navier–Stokes turbulence description (Chen, Ren & Meng Reference Chen, Ren and Meng2007), but this approach does not capture the detailed dynamical characteristics of the flow under heating conditions, the self-sustaining mechanisms of the flow and transitions between the flow regimes, that we wish to retain and optimise here. Motoki, Kawahara & Shimizu (Reference Motoki, Kawahara and Shimizu2018) also adopted a variational method to find the optimal steady velocity field for plane Couette flow with the largest Nusselt number. They found the optimal flow state is composed of streamwise-independent rolls at Reynolds number
${\textit{Re}} \sim 10^1$
, but there appear smaller-scale hierarchical quasi-streamwise vortex tubes near the walls in addition to the large-scale rolls at
${\textit{Re}} \geqslant 10^2$
. Although Motoki et al. (Reference Motoki, Kawahara and Shimizu2018) performed calculations up to
${\textit{Re}}=10^4$
, their analysis assumes a time-independent velocity field.
Optimisations need to be extended to include both the momentum equation for the velocity field, in which the buoyancy term may affect the flow regime, and time-dependent flows, such as turbulence. This requires a new framework that includes the dynamical effects of the flow on the mean heat transfer. The fully nonlinear variational method has been used in isothermal pipe flow by Pringle & Kerswell (Reference Pringle and Kerswell2010) to find initial flow perturbations that grow maximally. The smallest perturbation which triggers transition is called the ‘minimal seed’. For a review of wider applications, see Kerswell (Reference Kerswell2018). In the context of pipe flow, this framework has been successfully employed to find the minimal seed under various conditions affecting the flow (Pringle & Kerswell Reference Pringle and Kerswell2010; Pringle, Willis & Kerswell Reference Pringle, Willis and Kerswell2012; Marensi et al. Reference Marensi, Willis and Kerswell2019) and extended to induce transition ‘the other way’, i.e. with the aim to construct an optimal ‘baffle’ that destabilises turbulence, thereby causing transition from turbulence to the laminar state (Ding et al. Reference Ding, Marensi, Willis and Kerswell2020; Marensi et al. Reference Marensi, Ding, Willis and Kerswell2020). The minimal seed for transition in vertical heated pipe flow has been calculated using the model of § 2 (Chu, Willis & Marensi Reference Chu, Willis and Marensi2024). Here, we extend this new nonlinear variational framework to maximise the heat transfer. While previous optimisations in this geometry have identified optimal stationary velocity fields that maximise the heat transfer, here, we seek to optimise a time-independent body force that modifies the time-dependent flow, which is subject to the full governing equations. Although hard to accurately reproduce body forces in an engineering application, this is a step towards guiding such an approach. It should also be noted that the present study mainly focuses on improving heat transfer and understanding the physical mechanism – changes in the pumping power are not considered, although they are indirectly limited by the amplitude limit we apply to the force.
The plan of the paper is as follows. In § 2 we present our model of vertical heated pipe flow used for direct numerical simulation (DNS) and the variational equations used in the optimisation. In § 3, we first show behaviour for preliminary optimisations, including the features of optimal force, and the effects of variation in the target time. Optimisations are then performed in the laminarisation regime, shear-turbulence regime and convective-turbulence regime. Finally, the paper concludes with a summary in § 4.
2. Formulation
Consider flow upwards through a vertical pipe that passes through a hot room or chamber. The fluid will reach the ambient temperature exponentially along the pipe, but will typically do so over sufficiently long distances that the mean temperature can be approximated as linear along subsections of the pipe. The local temperature gradient depends on the rate of heat transfer, and hence on whether the flow is laminar or turbulent.
2.1. The heated pipe flow model
We begin with the model of Chu, Marensi & Willis (Reference Chu, Marensi and Willis2025), which assumes a linear mean temperature over the section of pipe considered. Axial periodicity of the velocity and temperature fluctuations is also assumed, so that the model should be applied to a section downstream of a bend or inlet effects. A particular feature of the model is that the mean axial temperature gradient is allowed to vary in time, reflecting the flow-dependent nature of heat transfer, e.g. the flow can transition from the shear-driven turbulent state to the convective state causing a significant drop in the heat flux. Good correspondence with previous models and experimental results under statistically steady conditions are shown in figure 1(b), where
$Bo$
is an empirically determined buoyancy parameter for which collapse over a large range of parameter sets is observed (Jackson Reference Jackson2006).
$(a)$
Regions of laminar flow (L), shear turbulence (S) and convective turbulence (C); SL and SC indicate that the flow may be found in either of the two states and
${\textit{Re}}$
and
$C$
are the Reynolds number and buoyancy parameters.
$(b)$
Heat flux measured by
$\textit {Nu}$
, normalised by the value in the isothermal limit (
$C\to 0$
), as a function of
$Bo = 8\times 10^4\,(Nu\,Gr)/({\textit{Re}}^{3.425}Pr^{0.8})$
where
$\textit{Gr}$
and
$\textit{Pr}$
are the Grashof and Prandtl numbers. Present data from simulations at
${\textit{Re}} = 5300$
,
${\textit{Pr}} = 0.7$
and various
$Gr = 16\,{\textit{Re}}\,C$
. Upper and lower branches correspond to shear turbulence and convective turbulence respectively (reproduced from Chu et al. Reference Chu, Marensi and Willis2025).
Using cylindrical coordinates
$\boldsymbol{x}=(r,\phi ,z)$
for a pipe of radius
$R$
, the total temperature is decomposed as
\begin{equation} T_{\textit{tot}}(\boldsymbol{x},t)=T_w(z,t) + T(\boldsymbol{x},t) - 2 T_b ,\end{equation}
where the wall temperature is given by
$T_w(z,t)=a_{\textit{tot}}(t)z+b$
, with temperature gradient
$a_{tot}(t)$
and constant reference temperature
$b$
, and
$T(\boldsymbol{x},t)$
carries the temperature fluctuations. Heat enters at the hot wall, then exits via a sink term proportional to
$a_{\textit{tot}}(t)$
(see (2.3) below). Let
$\langle {\boldsymbol{\cdot }}\rangle$
denote integration over the volume
$V$
and
$(1/V)\langle {\boldsymbol{\cdot }}\rangle$
be the spatial average. Temperature fluctuations in the model have a fixed positive spatial average
$T_b=(1/V)\langle T\rangle$
, maintained by adjustments in the temperature gradient
$a_{\textit{tot}}(t)$
. As
$T_b$
is fixed, the temperature difference between the wall and bulk is constant in the model. The constant factor
$-2T_b$
has been inserted in (2.1) so that the fluctuations
$T(\boldsymbol{x},t)$
are positive and are largest at the heated wall, where
$T|_{r=R}=2T_b$
. The heat flux through the wall per unit area is
$q_w=\lambda \overline {\left .({\partial T}/{\partial r})\right |_{r=R}}$
, where
$\lambda$
is the thermal conductivity of the fluid and the overline denotes the time average. For laminar flow, radial heat flux is purely conductive, as the flow is purely streamwise. The key quantity that measures the observed heat flux through the wall relative to the value for laminar flow is the Nusselt number
\begin{equation} \textit {Nu} = \frac {2R\,q_w}{\lambda \,(T|_{r=R}-T_b)} \, . \end{equation}
Whilst this definition for
$\textit {Nu}$
is conventional for pipe flows, it is not precisely
$1$
for the laminar case, mostly due to a geometric scale factor. (Also,
$T(r)$
and thereby
$q_w$
for the laminar case are not trivially determined, hence the practical dimensional-based approximation.) However,
$\textit {Nu}$
will always be presented in the ratio with the value observed in the absence of the body force, denoted
${\textit {Nu}}_{F=0}$
, so that such factors drop out. We will sometimes use the instantaneous
$\textit {Nu}(t)$
, where the time average is dropped.
Non-dimensionalisation using the temperature scale
$2T_b$
, the length scale
$R$
and velocity scale
$2U_b$
, where
$U_b$
is the mean flow speed, and assuming the Boussinesq approximation, leads to the dimensionless governing equations
\begin{align} \frac {\partial T }{\partial t}+({\boldsymbol{u}_{\textit{tot}}}\boldsymbol{\cdot } \boldsymbol{\nabla }) T &=\frac {1}{{\textit{Re}}\,Pr}{\nabla}^{2}T -{\boldsymbol{u}_{\textit{tot}}}\boldsymbol{\cdot }{\hat {\boldsymbol{z}}}\,a_{\textit{tot}}(t) , \end{align}
\begin{align} \frac {\partial \boldsymbol{u}_{\textit{tot}}}{\partial t}+(\boldsymbol{u}_{\textit{tot}}\boldsymbol{\cdot }\boldsymbol{\nabla }) {\boldsymbol{u}_{\textit{tot}}}&=-\boldsymbol{\boldsymbol{\nabla }} p_{\textit{tot}} + \frac {1}{{\textit{Re}}}{\nabla}^{2} {\boldsymbol{u}_{\textit{tot}}}+\frac {4}{{\textit{Re}}}(1+\beta '(t)+C T )\hat {\boldsymbol{z}} , \end{align}
\begin{align} \boldsymbol{\nabla } \boldsymbol{\cdot } \boldsymbol{u}_{\textit{tot}}&=0 , \end{align}
where
$\boldsymbol{u}_{\textit{tot}}(\boldsymbol{x},t)$
is the velocity field. The dimensionless boundary condition for the temperature is then
$T=1$
and the no-slip condition is applied to the velocity,
$\boldsymbol{u}_{\textit{tot}}=\boldsymbol{0}$
, at
$r=1$
. Axial periodicity over a distance
$L$
is assumed for the temperature fluctuations and velocity. The dimensionless fixed bulk temperature and flow rate are, respectively,
$(1/V)\langle T\rangle =\textstyle {1/2}$
and
$(1/V)\langle \boldsymbol{u}_{\textit{tot}}\boldsymbol{\cdot }\hat {\boldsymbol{z}}\rangle =\textstyle {1/2}$
, and these two conditions determine values for
$a_{\textit{tot}}(t)$
and the excess pressure fraction
$\beta '(t)$
via the spatial averages of the respective governing equation. The Reynolds and Prandtl numbers are
${\textit{Re}}=2U_bR/\nu$
and
${\textit{Pr}}=\nu /\kappa$
, where
$\nu$
and
$\kappa$
are the kinematic viscosity and thermal diffusivity respectively. The third dimensionless parameter
\begin{equation} C=\frac {Gr}{16\,{\textit{Re}}} \, , \end{equation}
measures the buoyancy force relative to the pressure gradient force that drives laminar isothermal flow,
$\textrm{d}p_{\textit{tot}}/\textrm{d}z=4/{\textit{Re}}$
, where
$Gr={\gamma g(T|_{r=R}-T_b)(2R)^{3}}/{\nu ^{2}}$
is the Grashof number,
$\gamma$
is the coefficient of volume expansion and
$g$
is gravitational acceleration. The parameters
$Gr$
and the Richardson number
$Ri=Gr/{\textit{Re}}^2$
are often used in systems involving buoyancy, but in this this application the pressure gradient force scales inversely with
${\textit{Re}}$
, and the parameter
$C$
is more readily When
$C=1$
, buoyancy forces are expected to be of similar magnitude to the pressure gradient force, and fundamental changes in flow are observed for
$C=O(1-10)$
. Figure 1(a) shows the regimes of flow for the model. The switch at
$C\approx 4$
is due to a convection-driven linear instability to helical modes.
For this configuration, the laminar solution does not have a simple analytic form, and must be computed numerically. We consider perturbations from the laminar state, subscripted with
$0$
, and decompose variables as
${\boldsymbol{u}_{\textit{tot}}}=\boldsymbol{u}_0+\boldsymbol{u}$
,
$p_{\textit{tot}}=p_0+p$
,
$1+\beta '(t)=1+\beta _{0}+\beta (t)$
,
$T=\mathit{\varTheta }_0+\mathit{\varTheta }$
,
$a_{\textit{tot}}(t)=a_{0}+a(t)$
. The perturbations then satisfy
\begin{align} \frac {\partial \mathit{\varTheta } }{\partial t}+u_0\frac {\partial \mathit{\varTheta } }{\partial z}+u_r\frac {d\mathit{\varTheta }_{0}} {\textit{dr}} + (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }) \mathit{\varTheta } &=\frac {1}{{\textit{Re}}{\textit{Pr}}}{\boldsymbol{\nabla }} ^{2}\mathit{\varTheta } -u_z a_{0}-(u_0+u_z)a(t) , \end{align}
\begin{align} \frac {\partial \boldsymbol{u}}{\partial t}+u_0\frac {\partial \boldsymbol{u}}{\partial z}+u_{r}\,\frac {du_0}{\textit{dr}}\, \hat {\boldsymbol{z}}+(\boldsymbol{u}\boldsymbol{\cdot } \boldsymbol{\nabla }) \boldsymbol{u} &= -\boldsymbol{\nabla } p + \frac {1}{{\textit{Re}}}{\boldsymbol{\nabla }} ^{2}\boldsymbol{u}+\frac {4}{{\textit{Re}}}(C \mathit{\varTheta }+\beta (t))\hat {\boldsymbol{z}}, \end{align}
\begin{align} \boldsymbol{\nabla } \boldsymbol{\cdot } \boldsymbol{u}&=0 , \end{align}
where
$\boldsymbol{u}=(u_r,u_{\phi },u_z)$
. Further details on the numerical model can be found in Chu et al. (Reference Chu, Willis and Marensi2024) and Chu et al. (Reference Chu, Marensi and Willis2025).
2.2. Variational optimisation of a body force for heat transfer
In the following, we suppose that a body force
$\boldsymbol{F}(\boldsymbol{x})$
is appended to the right-hand sides of the Navier–Stokes equations (2.4) and (2.8), then seek to optimise the form of
$\boldsymbol{F}(\boldsymbol{x})$
such that it maximises
$\textit {Nu}$
, subject to a constraint on the magnitude of
$\boldsymbol{F}$
. The scalar dissipation
${({1}/{2})}\lambda (\boldsymbol{\nabla } T)^2$
will be used as a proxy for the heat transfer (Batchelor Reference Batchelor1959). (The relationship between
$\textit {Nu}$
and the scalar dissipation is shown in Appendix A.) The term ‘entransy dissipation’ is also being used in thermal science (Guo et al. Reference Guo, Zhu and Liang2007; Kostic Reference Kostic2017). Here, the time-averaged quantity
\begin{equation} J=\frac {1}{\mathcal T}\int _{0}^{\mathcal T}\frac {1}{2} \left\langle ({\boldsymbol{\nabla }} T)^2\right\rangle \, \textrm{d}t, \end{equation}
is used as our objective function. The way in which the scalar dissipation is used depends upon the thermal boundary condition (Appendix A). As the fixed temperature boundary condition is applied, maximal
$J$
corresponds to maximised ‘dissipation’ of
$T^2$
and larger heat flux. (For the fixed heat flux boundary condition, minimal
$J$
corresponds to minimal thermal resistance within the body, and hence maximum heat flux.)
In this work, we seek for the first time an optimal that is subject to the full governing equations for the velocity field. A Lagrangian is constructed as follows:
\begin{align} {\mathcal L} &=\frac {1}{N}\sum _{i=1}^{N}J_i -\lambda _0\left (\left \langle \frac {1}{2}(\boldsymbol{F})^{2}\right \rangle -A_{0}\right )-\sum _{i=1}^{N}\int _{0}^{\mathcal T}\langle \boldsymbol{v}_i\boldsymbol{\cdot }(\textrm{NS}(\boldsymbol{u}_i))\rangle \textrm{d}t \nonumber\\&\quad - \sum _{i=1}^{N}\int _{0}^{\mathcal T}\langle \varPi _i (\boldsymbol{\nabla }\boldsymbol{\cdot } \boldsymbol{u}_i)\rangle \textrm{d}t - \sum _{i=1}^{N}\int _{0}^{\mathcal T}\langle \pi _i (\textrm{Tem}(\mathit{\varTheta }_i))\rangle \textrm{d}t \nonumber\\&\quad - \sum _{i=1}^{N}\int _{0}^{\mathcal T}\varGamma _i\langle (\boldsymbol{u}_i\boldsymbol{\cdot } \hat {\boldsymbol{z}}) \rangle \textrm{d}t-\sum _{i=1}^{N}\int _{0}^{\mathcal T}Q_i\langle (\mathit{\varTheta }_i) \rangle \textrm{d}t. \end{align}
As the initial condition for the velocity field could be turbulent, to improve the robustness of the results we apply the optimisation to
$N$
initial velocity fields. The variables
$\lambda _0$
,
$\varPi _i$
,
$\pi _i(\boldsymbol x,t)$
,
$\varGamma _i(t)$
,
$Q_i(t)$
and
$\boldsymbol{v}_i(\boldsymbol{x},t)=(v_{r, i},v_{\phi , i},v_{z, i})$
are Lagrange multipliers. The first term, the ensemble average of the time-averaged scalar dissipation, is the objective function to be maximised. The second term fixes the amplitude of the body force; the spatial average of
$|\boldsymbol{F}|$
is given by
$\sqrt {2A_0/V}$
. Next, the velocity perturbation
$\boldsymbol{u}$
is constrained to satisfy the Navier–Stokes equation
$\textrm{NS}(\boldsymbol{u})$
and the continuity equation, and the temperature perturbation satisfies the temperature equation
$\textrm{Tem}(\mathit{\varTheta })$
, each over the period from
$t=0$
to
$t=\mathcal T$
. The last two terms ensure that the velocity satisfies the fixed mass flux and that the bulk temperature is fixed.
Taking variations of
$\mathcal L$
with respect to each variable and setting them equal to zero, we obtain the following set of Euler–Lagrange equations. The adjoint Navier–Stokes, temperature equation and continuity equations are
\begin{align} &\quad\frac {\partial \mathcal L} {\partial \boldsymbol{u}_i} = \frac {\partial \boldsymbol{v}_i}{\partial t}+u_0\frac {\partial \boldsymbol{v}_i}{\partial z}-v_{z,i}{u_0}'\boldsymbol{\hat {r}}+ \boldsymbol{\nabla }\times (\boldsymbol{v}_i\times \boldsymbol{u}_i)- \boldsymbol{v}_i\times \boldsymbol{\nabla } \times \boldsymbol{u}_i +\boldsymbol{\nabla } \varPi _i\nonumber \\ &\qquad\quad\,\,\, + \frac {1}{{\textit{Re}}}{\nabla}^{2}\boldsymbol{v}_i-\pi _i\mathit{\varTheta }_{0}'\hat {\boldsymbol{r}}-\pi _i\boldsymbol{\nabla }\mathit{\varTheta }_i-\pi _i (a(t)+a_0(t))\hat {\boldsymbol{z}}-\varGamma _i \hat {\boldsymbol{z}}=0, \end{align}
\begin{align} &\frac {\partial \mathcal L}{\partial \mathit{\varTheta }_i}=\frac {\partial \pi _i}{\partial t}+u_0\frac {\partial \pi _i}{\partial z}+\frac {4}{{\textit{Re}}} v_{z,i} C+\boldsymbol{u}_i\boldsymbol{\cdot } \boldsymbol{\nabla }\pi _i + \frac {1}{{\textit{Re}} {\textit{Pr}}}{\nabla}^{2}\pi _i-Q_i-\frac {1}{\mathcal T }{\nabla}^{2} T_i=0 , \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad\quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}_i=0 . \end{align}
The compatibility conditions (terminal conditions for backward integration of (2.12) and (2.13)) are given by
\begin{align} \frac {\delta \mathcal L}{\delta \boldsymbol{u}_i(\boldsymbol x,\mathcal T)}&=-\boldsymbol{v}_i(\boldsymbol x,\mathcal T)=0, \end{align}
\begin{align} \frac {\delta \mathcal L}{\delta \mathit{\varTheta }_i(\boldsymbol x,\mathcal T)}&=-\pi _i(\boldsymbol x,\mathcal T)=0 ,\end{align}
and the optimality condition is
\begin{equation} \frac {\delta \mathcal L}{\delta \boldsymbol{F}}=-\lambda _0 \boldsymbol{F} +\frac {1}{N} \sum _{i=1}^{N}\int _{0}^{\mathcal T}\boldsymbol{v}_i\,\textrm{d}t =0. \end{equation}
For an arbitrary initial
$\boldsymbol{F}$
and set of initial conditions
$\boldsymbol{u}_i$
, the force
$\boldsymbol{F}$
is incrementally updated to produce a maximum in
$\mathcal{L}$
where
$ \delta \mathcal L/\delta \boldsymbol{F}$
should vanish. An iterative algorithm similar to that in Pringle et al. (Reference Pringle, Willis and Kerswell2012) is applied. The update for
$\boldsymbol{F}$
at
$(j+1)\textrm{th}$
iteration is
\begin{equation} \boldsymbol{F}^{(j+1)}=\boldsymbol{F}^{(j)}{+}\epsilon _0\frac {\delta \mathcal L}{\delta \boldsymbol{F}^{(j)}}, \end{equation}
where
$\epsilon _0$
is a small value, controlled using a procedure described in Pringle et al. (Reference Pringle, Willis and Kerswell2012) and
$\lambda _0$
is determined by the condition that
$\langle [\boldsymbol{F}(\boldsymbol{x})^{(j+1)}]^2\rangle =2\,A_0$
.
2.3. Numerical methods
Calculations are carried out using the open-source code Openpipeflow (Willis Reference Willis2017). Variables are discretised in the domain
$\left \{ r,\phi ,z \right \}=[0,1]\times [0,2\pi ]\times [0,2\pi /\alpha ]$
, where
$\alpha =2\pi /L$
, using Fourier decomposition in the azimuthal and streamwise direction and finite difference in the radial direction, e.g.
\begin{equation} \boldsymbol{u}(r_s,\phi ,z)=\sum _{k\lt |K|} \sum _{m\lt |M|} \boldsymbol{u}_{skm}\textrm{e}^{\textrm{i}(\alpha kz+m\phi )} \, , \qquad s = 1,\ldots , S ,\end{equation}
where the radial points
$r_s$
are clustered towards the wall. Temporal discretisation is via a second-order predictor–corrector scheme, with Euler predictor for the nonlinear terms and Crank–Nicolson corrector. To keep the nonlinear optimisations manageable, a Reynolds number
${\textit{Re}}=3000$
and
${\textit{Pr}}=0.7$
are adopted with a domain of length
$L=10$
radii. We use mesh resolution of
$S=64, M=48, K=42$
and the size of the time step is
$\Delta t=0.01$
. This resolution is sufficient to maintain a drop-off in the amplitude of the coefficients by three to four orders magnitude, which experience has shown to be sufficient for accurate simulation of shear-driven turbulence. For the values of
$C$
considered here, the convective state is less computationally demanding to simulate.
3. Results
In order to assess the parameters of the optimisation method, we first show preliminary optimisations in § 3.1. Then, we optimise the body force to maximise the heat transfer in the three typical flow regimes of vertical heated pipe flow, i.e. the laminarisation regime (§ 3.2), shear-turbulence regime (§ 3.3) and convective-turbulence regime (§ 3.4). (Parameter regimes for this model are shown in figure 1).
3.1. Preliminary optimisation
For the laminar case, the state is unique and we require only one initial velocity field,
$N=1$
. We start with the unforced laminarised flow at
${\textit{Re}}=3000$
,
$C=3$
and take random fields for the initial force (such as a turbulent velocity field). Results from a preliminary optimisation with
$A_0=5 \times 10^{-7}, \mathcal T=400$
are shown in figure 2. The instantaneous Nusselt number,
$\textit {Nu}(t)$
, normalised by the mean for the unforced flow,
$\textit {Nu}_{F=0}$
, is shown in figure 2
$(a)$
for each iteration. The final value increases by more than 80 % over the unforced laminar case. Figure 2
$(b)$
shows the objective function
$J$
(2.10), the (time-averaged) Nusselt number
$\textit {Nu}$
(2.2) and the final value of the instantaneous Nusselt number
${\textit {Nu}}(\mathcal T)$
, versus iteration, normalised by values for the unforced case. Changes in these quantities show good agreement, indicating that the global (volume integrated) scalar dissipation quantity effectively captures the local (boundary) heat-transfer behaviour measured by the Nusselt number. Figure 2
$(c)$
shows the residual of the calculation, which drops quickly in first
$50$
iterations then more gradually. (Spikes are related to the method that seeks to increase
$\epsilon _0$
as much a possible, which affects the magnitude of the residual via
$\lambda _0$
.) Usually, the optimisation is stopped when the change in the Nusselt number drops below
$10^{-5}$
. Figure 2
$(d)$
tracks the amplitude of the three components of the body force versus iteration. The amplitude of the streamwise component drops significantly, while the amplitude of the cross-stream components increases. This suggests that the cross-stream components of the body force play a dominant role in enhancing heat transfer, whereas the contribution of the streamwise component is nearly negligible.
$(a)$
Instantaneous Nusselt number
$\textit {Nu}(t)$
as the iteration proceeds, where
${\textit {Nu}}_{F=0}$
refers to the value of the unforced case.
$(b)$
Dimensionless time-averaged scalar dissipation, (time-averaged)
$\textit {Nu}$
and
$\textit {Nu}(\mathcal T)$
versus iteration, normalised by their values at the zeroth iteration.
$(c)$
The residual
$\langle ({\delta \mathcal L}/{\delta \boldsymbol F})^2\rangle$
versus iteration.
$(d)$
Magnitude of the components of body force
$({1}/{2})\langle F_r^2\rangle , ({1}/{2})\langle F_\phi ^2\rangle , ({1}/{2}) \langle F_z^2\rangle$
versus iteration. The optimisation is run at
$A_0=5\times 10^{-7}, \mathcal T=400, C=3, {\textit{Re}}=3000$
.
If the time horizon
$\mathcal{T}$
is long enough, then
$\textit {Nu}$
is optimised for the steady response to the force. Figure 2(a) suggests that
$\mathcal{T}$
should be greater than
$200$
. Indeed, the form of the optimal is found to change when increasing
$\mathcal{T}$
from
$50$
to
$100$
, and again to
$200$
, but upon further increase to
$600$
, the optimal is essentially the same up to a rotation, as shown in figure 3. Interestingly, the optimal body force optimised for a short target time has perfect rotational symmetry. That the force induces rolls is consistent with the steady velocity fields computed by Meng et al. (Reference Meng, Liang and Li2005). The time horizon
$\mathcal{T}$
is an extra parameter here, and as it is increased, the azimuthal wavenumber
$m$
decreases. This is consistent with smaller-scale vortices growing more rapidly (Schmid Reference Schmid2007), thereby increasing the heat transfer within a shorter time. With a longer target time, however, the larger-scale mode is more effectively amplified for the given magnitude of force.
It is found that the distribution of the body force is almost uniform in the streamwise direction. Although streamwise dependence of the force was included as a possibility in the preliminary calculations, it was not expected – structures sweep rapidly through the domain, leading to streamwise independence in the time-integrated update to
$\boldsymbol{F}^{(j)}$
in the optimisation. In the isothermal calculations of Marensi et al. (Reference Marensi, Ding, Willis and Kerswell2020), streamwise-dependent forces were not identified even with domains 5 times longer. Therefore, to simplify the form of the body force and to accelerate convergence, we constrain the body force to be streamwise-independent in the optimisations of the following sections by zeroing coefficients of the streamwise-dependent Fourier modes.
Optimised body forces for target time horizons
$(\textit a)$
$\mathcal T=50, (\textit b)$
$\mathcal T=100, (c)$
$\mathcal T=200, (\textit d)$
$\mathcal T=600$
at
$C=3, A_0=5e-7, {\textit{Re}}=3000, N=1$
. Colours indicate the streamwise component, while arrows represent the cross-stream components of the body force. Panels
$(\textit {a})$
and
$(\textit {b})$
share the first colour bar,
$(\textit {c})$
and
$(\textit {d})$
share the second. The largest arrow has magnitude
$2.38 \times 10^{-4}$
in
$(\textit a)$
,
$2.36 \times 10^{-4}$
in
$(\textit b)$
,
$2.93 \times 10^{-4}$
in
$(\textit c)$
and
$3.21 \times 10^{-4}$
in
$(\textit d)$
.
Instantaneous Nusselt number normalised by the value for unforced flow for several force amplitudes, starting from a laminar initial condition at
$C=3, {\textit{Re}}=3000$
. The vertical dashed line indicates the optimisation target time
$\mathcal T=600$
.
3.2. Optimisation in the laminarisation regime
3.2.1. Optimal force for the laminar state
Having examined properties of the parameters necessary for optimisation, in this section we consider optimisation in the laminarised regime at
$C=3,{\textit{Re}}=3000$
in more detail. In particular, we examine the dependence of the rotational symmetry on
$A_0$
and the presence of local optimals (dependent on the initial guess for the force).
For the laminar initial condition, as we have assumed streamwise independence for the force, a small streamwise-dependent perturbation must be added to the initial velocity so that transition to turbulence may be triggered if the resulting two-dimensional flow is unstable. We add a perturbation of magnitude
$E_0=({1}/{2})\langle \boldsymbol{u}^2\rangle \approx 10^{-7}$
and set a longer
$\mathcal{T}=600$
to allow the occurrence of transition. Figure 4 shows the instantaneous Nusselt number,
$\textit {Nu}(t)$
, for several
$A_0$
. Optimisation improves the heat transfer substantially: for
$A_0=10^{-7}$
heat transfer is almost 50 % greater than that of the unforced flow, and for
$A_0=5\times 10^{-7}$
is almost doubled. For slightly higher
$A_0=6\times 10^{-7}$
,
$\textit {Nu}(t)$
experiences a sudden increase near the end of the optimisation target time and fluctuates thereafter, indicating the onset of turbulence. (At
$A_0=5\times 10^{-7}$
, transition is observed very late, at around
$t=1000$
, and interestingly, the transition does not lead to a larger
$\textit {Nu}$
. This phenomenon will be discussed later.)
The typical amplitude-dependent form of the forces obtained from optimisations are given in figure 5(a–d). (As the axial component is at least an order of magnitude smaller, it is not shown.) At small
$A_0$
, the body force has a single pair of rolls. At increased force amplitude, figure 5(b) illustrates how the vortex structure gradually approaches the wall, reducing the spatial scale in both the radial and spanwise directions. This is actually found to be a local optimal for this
$A_0$
, as two pairs of rolls may be squeezed in to increase
$\textit {Nu}$
a little further. At larger
$A_0$
, more rolls are seen in figure 5(c–d). For the largest
$A_0$
, turbulence is triggered within
$\mathcal{T}$
, and the optimisation struggles to converge to a well-organised optimal force. However, a preference for roll structures of larger
$m$
is clear. The form and increase in wavenumber is consistent with the calculations of Meng et al. (Reference Meng, Liang and Li2005), Jia et al. (Reference Jia, Liu, Liu and Nakayama2014) and Wang et al. (Reference Wang, Liu and Liu2015) for steady flow.
Optimised body forces at
$(\textit a)A_0=10^{-8}, (\textit b)A_0=10^{-7}, (\textit c)A_0=5\times 10^{-7}, (\textit c)A_0=6\times 10^{-7}$
, starting from a laminar initial condition at
$C=3, \mathcal T=600, {\textit{Re}}=3000$
. The largest arrow has magnitude
$4.85 \times 10^{-5}$
in
$(\textit a)$
,
$1.87 \times 10^{-4}$
in
$(\textit b)$
,
$4.05 \times 10^{-4}$
in
$(\textit c)$
and
$4.26 \times 10^{-4}$
in
$(\textit d)$
. As the axial component is at least an order of magnitude smaller, it is not shown.
We have observed that optimal body forces display rotational symmetry of different azimuthal wavenumbers
$m$
in figures 3 and 5. We let
$O_2$
denote an optimal with
$2$
-fold rotational symmetry. This optimal may have non-zero Fourier coefficients only for
$m=0,2,4,6,8,\ldots$
, but note that this does not exclude a force with only non-zero modes
$m=0,4,8,\ldots$
, which corresponds to an optimal
$O_4$
with
$4$
-fold rotational symmetry. To examine the influence of the azimuthal periodicity, and to simplify the optimisation further, we consider optimal forces restricted to azimuthal Fourier modes of wavenumbers
$0$
and
$m$
only, and denote them
$O_{Fm}$
. Note that rotational symmetry is imposed only on the force, and not on the velocity field.
Optimisations have been computed for
$O_{Fm}$
, then used as starting forces for optimisations in the full space
$O_1$
. In this way, we examine the dependence on
$m$
to determine which rotational symmetry is the global optimal. The Nusselt numbers of the final states,
$\textit {Nu}(\mathcal{T})$
, as a function of iteration are shown in figure 6(
$\textit {a}$
) starting from the
$O_{Fm}$
forces. Three optimisations starting from random initial forces are also shown. Several observations can be made. Firstly, there are multiple local optimals
$O_m$
, of which
$O_2$
(figure 5
c) is the global optimal for this
$A_0$
. Secondly, the optimal of type
$O_m$
does not produce much greater
$\textit {Nu}$
than the optimal
$O_{Fm}$
, of the reduced Fourier space, used as the starting force. Thirdly, if the starting force is quite perturbed, such as for the random initial forces, it is most likely to end up at the global optimal. Similarly, the optimisation starting from
$O_{F4}$
(light green line) appears to pass close to
$O_1$
, but ends at the global optimal
$O_2$
.
(a) Nusselt number of final state
$\textit {Nu}(\mathcal{T})$
versus iteration for different initial forces, with
$A_0=10^{-7}, C=3, {\textit{Re}}=3000$
. The legend indicates the initial force, and the resulting local optimal forces are labelled on the curves. (b) Normalised (time-averaged) Nusselt numbers at
$C=3, {\textit{Re}}=3000$
versus
$O_{Fm}$
for force amplitudes
$A_0$
the legend. The global optimal is highlighted with a dashed circle. For dashed lines, the forced flow state remains laminar. For solid lines, the forced state is turbulent and the values are time averages.
As we have observed that
$O_m$
does not produce much greater
$\textit {Nu}$
than
$O_{Fm}$
, for small
$A_0$
at least, we directly compare
$\textit {Nu}$
of the flow forced by
$O_{Fm}$
for several
$m$
to determine the global optimal. Figure 6(b) shows the Nusselt number, calculated using averages over
$5000$
time units for each simulation. For the dashed lines, the final state is still laminar and good convergence is easily achieved. For the solid lines at larger
$A_0$
, the final state is turbulent, which renders convergence difficult. For the latter case, and for convenience in this section, the force has been calculated using an artificially stabilised two-dimensional flow, by putting
$K=1$
in (2.19) (equivalent to adding no three-dimensional perturbation,
$E_0=0$
). The optimal force leads to an artificially stabilised optimal velocity field, as for the optimal steady velocity fields reported by Meng et al. (Reference Meng, Liang and Li2005), Jia et al. (Reference Jia, Liu, Liu and Nakayama2014) and Wang et al. (Reference Wang, Liu and Liu2015). The force is then applied, here resulting in a fully three-dimensional time-dependent turbulent simulation, from which
$\textit {Nu}$
is calculated. (In principle, the stabilisation of the flow during optimisation may render the force no longer optimal. This is examined further for each initial flow regime, and is found to be a good approach for intermediate
$C$
. We will show this for the convective case in § 3.4). As
$A_0$
increases, the rotational symmetry
$m$
of the (global) optimal increases, from
$m=1$
at
$A_0=10^{-8}$
, to
$m=2$
at
$A_0=10^{-7}$
and is
$m\approx 5$
at substantially larger
$A_0=10^{-5}$
. However, it should be noted that
$\textit {Nu}$
is not strongly dependent on
$m$
.
3.2.2. The path to transition and effect on heat transfer
Figure 7(a) shows the instantaneous Nusselt number
$\textit {Nu}(t)$
for flows forced by
$O_{F1}$
at several force amplitudes. At small
$A_0$
, the flow is reshaped into a two-dimensional forced laminar state. With an increase of force amplitude to
$A_0=6\times 10 ^{-7}$
, the flow does not transition to turbulence directly – instead, the two-dimensional state quickly forms, then transitions later to a travelling wave solution. Figure 8 shows the two-dimensional reshaped laminar solution and a travelling wave solution found in the flow forced by
$O_{F1}$
at
$A_0=6\times 10^{-7}$
(
$t\approx 500$
and
$t\approx 1500$
of figure 7(a) respectively). In isothermal flow, forces have been used to find travelling wave solutions, via homotopy (Faisst & Eckhardt Reference Faisst and Eckhardt2003; Wedin & Kerswell Reference Wedin and Kerswell2004). The travelling wave solution has larger
$\textit {Nu}$
compared with the two-dimensional reshaped laminar solution. With further increase of force amplitude, the flow transitions from the two-dimensional state to a (mildly) chaotic three-dimensional state, and finally also converges to a travelling wave state at a later time (not shown here). At
$A_0=10^{-5}$
, the flow directly transitions to a strong chaotic three-dimensional state, along with a greatly increased
$\textit {Nu}$
. However, this is not always the case, and in fact figure 7(b) shows that at larger
$A_0$
the transition from two-dimensional flow (stabilised by setting
$K=1$
) to the chaotic three-dimensional state leads to a decrease in
$\textit {Nu}$
, except for the case
$m=1$
.
Enhanced
$\textit {Nu}$
of forced flows at
$C=3, {\textit{Re}}=3000$
.
$(\textit a)$
Instantaneous
$\textit {Nu}(t)$
for flows forced by
$O_{F1}$
.
$(\textit b)$
Values of
$\textit {Nu}$
(time-averaged) for flows forced by
$O_{Fm}$
in the turbulent state (solid lines) and
$\textit {Nu}(\mathcal{T})$
for the stabilised streamwise-independent state (dashed lines).
Isosurfaces of streamwise velocity for
$(\textit a)$
the streamwise-independent state at
$t=500$
and
$(\textit b)$
the travelling wave solution at
$t=1500$
, for flow forced by
$O_{F1}$
with
$A_0=6\times 10^{-7}$
. Red/yellow are at 20 % of the min/max streamwise velocity.
Time series of (a) instantaneous
$\textit {Nu}(t)$
and (b) instantaneous energy of rolls
$E_{roll}=E(u_r)+E(u_\phi )$
,
$E(u_r)=({1}/{2})\langle u_r^2 \rangle$
,
$E(u_\phi )=({1}/{2})\langle u_\phi ^2 \rangle$
, for flow forced by
$O_{F4}$
at
$C=3, {\textit{Re}}=3000$
at two different
$A_0$
. Solid lines are for the forced turbulent state, and dashed lines are for the streamwise-independent (artificially stabilised) state.
Contours of the radial temperature gradient
${\partial T}/{\partial r}$
at the boundary at
$C=3, {\textit{Re}}=3000$
, normalised by the unforced laminar radial temperature gradient. The flow is forced by
$O_{F4}$
with
$(\textit {a,b})$
$A_0=10 ^{-6}$
and
$(\textit {c,d})$
$A_0=10 ^{-5}$
, the flow states are taken in
$(\textit {a,c})$
from the forced laminar state (
$t=1000$
) and
$(\textit {b,d})$
forced turbulent state (
$t=1000$
).
It is interesting that the more chaotic state does not necessarily lead to an improvement in heat transfer. For the case where the flow is forced by
$O_{F4}$
, the instantaneous Nusselt number
$\textit {Nu}(t)$
and the roll energy
$E_{roll}(t)=E(u_r)+E(u_\phi )$
are shown in figure 9 for two amplitudes
$A_0$
. At
$A_0=10^{-6}$
heat transfer increases after transition, but at
$A_0=10^{-5}$
it is reduced after transition. In both cases, the energy of rolls increases, but in the latter case by not as much. Stronger rolls are typically associated with enhanced heat transfer, but this shows that higher roll energy alone does not necessarily correspond to more efficient heat transfer. Figure 10 shows the contours of the radial temperature gradient
${\partial T}/{\partial r}$
evaluated at the boundary, for the forced laminar and forced turbulent states, normalised by the mean radial temperature gradient of the unforced flow. When the flow is forced by
$O_{F4}$
with
$A_0=10^{-6}$
, the laminar state (figure 10
a) exhibits distinct regions of strong and weak heat transfer. After the transition (figure 10
b), these regions are not fixed, and the regions of higher heat flux widen and intensify a little. At
$A_0=10^{-5}$
(figure 10
c–d), the regions of higher heat flux also widen after the transition, but are less organised and weaker than for the steady flow, despite the slightly increased roll energy. This suggests that the streamwise vortices in the forced turbulent states are not as efficient as those in the forced laminar states. The rolls are unsteady in forced turbulent states and move further from the wall intermittently due to the waving of low-speed streaks, leading to a non-uniform heat-transfer distribution. Heat transfer is enhanced on the sides where the rolls are positioned close to the wall, but weakens as the rolls move further away from the wall. Overall, the transition to the forced turbulent state involves two main competing effects on heat transfer: the first is the enhancement by the rolls, which increases heat transfer by facilitating better mixing, the second is unsteadiness of the rolls, which at larger
$A_0$
can reduce heat transfer due to inconsistency in their position near the wall. The case
$O_{F1}$
is an exception in figure 7(b), where the transition to a chaotic state leads to an increase in
$\textit {Nu}$
– only one pair of rolls is inefficient and the unsteadiness can lead to the creation of additional rolls.
In this section, we have shown that optimisations in the laminarised regime can substantially enhance heat transfer, with the optimal forcing inducing axially aligned rolls of increasing azimuthal symmetry – though with only a weak dependence of
$\textit {Nu}$
on
$m$
. Interestingly, triggering turbulence does not necessarily lead to an increase in
$\textit {Nu}$
, as the resulting unsteady rolls are less efficiently positioned for heat transfer.
3.3. Optimisation in the shear-turbulence regime
Returning to smaller
$C$
, the shear-turbulence case is the most challenging, as the flow state remains highly chaotic. As turbulence is already effective in enhancing heat transfer relative to the laminarised case, it is not obvious that optimisation should be able to improve heat transfer substantially.
We focus on the case
$C=1$
at
${\textit{Re}}=3000$
. For the highly chaotic flow, it is difficult to apply the method with a large target time (Pringle & Kerswell Reference Pringle and Kerswell2010; Pringle et al. Reference Pringle, Willis and Kerswell2012; Marensi et al. Reference Marensi, Willis and Kerswell2019), and convergence was found to fail for
$\mathcal{T}$
even as low as
$100$
. However, reasonably good convergence was found for
$\mathcal{T}=50$
. Although this is not sufficient time to capture the statistics of the end state, it is sufficient time for a response to be observed, so that it is reasonable to examine whether the heat transfer has been pushed in the right direction.
Instantaneous Nusselt number for flows forced by
$O_1$
at
$A_0=10^{-6}$
, where the legend indicates the number of initial velocity fields used in the optimisation. The time horizon
$\mathcal{T} =50$
is marked by the vertical dashed line.
Values of
$\textit {Nu}/\textit {Nu}_{F=0}$
for flows forced by optimal forces indicated. Here,
$N$
is the number of initial velocity conditions used in the optimisation.
Figure 11 shows the instantaneous Nusselt number
$\textit {Nu}(t)$
for flows forced by the full space optimal
$O_1$
with
$A_0=10^{-6}$
. The optimisation, although with a short target time, still increases the Nusselt number significantly. Surprisingly, table 1 shows that the (subsequent time-averaged)
$\textit {Nu}/\textit {Nu}_{F=0}$
value does not change significantly with more initial conditions
$N$
, despite that the short
$\mathcal{T}$
might suggest greater dependence on the initial condition, nor does
$\textit {Nu}$
vary significantly with rotational symmetry. Structures of the optimised forces are shown in figure 12 for
$N=1$
and
$3$
(the cases for
$N=9$
are almost identical to
$N=3$
). Although
$\textit {Nu}$
varies very little with
$N$
, for
$O_1$
the structure of the optimal force does change – going from
$N=1$
to
$3$
initial conditions, the force develops towards a structure more like the
$O_3$
optimal, but larger
$N$
slightly improves convergence. Optimising for
$O_3$
itself, the structure clearly becomes more regular as
$N$
is increased. For
$O_{F3}$
(third column),
$O_{5}$
and
$O_{F5}$
(not shown), the structure of the force does not change for larger
$N$
, so that only one initial velocity field is sufficient. This is reasonable for larger
$m$
, as the angular section
$[0,2\pi /m]$
of the force is determined by
$m$
angular subsections of the rotationally unconstrained velocity field.
Forces
$O_1$
,
$O_3$
and
$O_{F3}$
optimised using
$N$
initial velocity fields at
$A_0=10^{-6}$
,
$C=1$
. The largest arrow has magnitude
$6.95\times 10^{-4}$
in
$(\textit {a})$
,
$7.52\times 10^{-4}$
in
$(\textit {b})$
,
$5.76\times 10^{-4}$
in
$(\textit {c})$
,
$8.14 \times 10^{-4}$
in
$(\textit {d})$
,
$6.60\times 10^{-4}$
in
$(\textit {e})$
,
$5.76\times 10^{-4}$
in
$(\textit {f})$
.
Although we have computed optimals only for a short
$\mathcal{T}$
, our results strongly suggest that inducing rolls remains optimal even for flow already turbulent, in the shear-turbulence regime. As before, it is interesting to examine whether or not there is a strong dependence on the rotational symmetry
$m$
. Figure 13(
$\textit {a}$
) shows the Nusselt number for flows forced by
$O_{Fm}$
(solid) and
$O_{m}$
(dashed) for several
$m$
, using
$N=1$
for
$O_{Fm}$
and
$N=3$
for
$O_m$
. For small
$m=1,2$
, the force
$O_m$
appears to produce larger
$\textit {Nu}$
than
$O_{Fm}$
, but this is because the optimisation for
$O_1$
actually found a structure closer to that for
$O_3$
, previously seen in figure 12 (cf. (d) and (e)) and similarly, the
$O_2$
optimal is structurally more like
$O_4$
, seen in figure 14 (cf. (a) and (c)). Constraining the number of rolls strictly, by using the single Fourier mode, the rolls of
$O_{F3}$
are a little stretched (figure 14
e) relative to those of
$O_3$
(figure 14
b), but for
$m\gt 3$
there is essentially no visible difference between
$O_m$
and
$O_{Fm}$
. For
$m\geqslant 3$
,
$\textit {Nu}$
is relatively insensitive to the wavenumber (figure 13
$\textit {a}$
).
An interesting observation is the occurrence of laminarisation of the shear turbulence when forced by
$O_{F1}$
at low amplitudes
$A_0=10^{-7}$
. A very similar observation was reported by Willis, Hwang & Cossu (Reference Willis, Hwang and Cossu2010), where there the roll force was beneficial for drag reduction in isothermal flow. Here, as the turbulent flow enhances heat transfer, the laminarisation can lead to a reduction in the heat transfer, i.e.
$\textit {Nu}/\textit {Nu}_{F=0}\lt 1$
, see figure 13(
$\textit {b}$
).
Optimisation for a steady laminar flow can be imposed by setting
$K=1$
in (2.19), as applied in § 3.2. Figure 13(
$\textit {b}$
) compares optimisation with the short
$\mathcal{T}$
(solid lines) with steady laminar optimisation (dashed lines) with
$\mathcal{T}=600$
. Particularly for larger
$A_0$
, it should be noted that including the time dependence of the flow in the optimisation does improve the resulting
$\textit {Nu}$
over the steady assumption, despite the short
$\mathcal{T}$
. There is an exception for the
$m=1,2$
case at the largest
$A_0$
, however, where the short-time optimal force has an unusual structure, see, e.g. figure 14(d). For the structure of the optimal forces, when optimised for the steady two-dimensionalised state (figure 14
g–i) these forces may induce flows that are similar to those of previous calculations (Meng et al. Reference Meng, Liang and Li2005; Jia et al. Reference Jia, Liu, Liu and Nakayama2014; Wang et al. Reference Wang, Liu and Liu2015), although there the optimisations were independent of the flow regime.
Overall, it is found that for the shear-turbulent regime, rolls remain optimal, but the best choice of azimuthal wavenumber may differ from that in the laminarised flow regime, i.e. there is dependence on
$C$
through the governing equations. Including time dependence in the optimisation, the rolls are notably closer to the wall (figure 14
e,f). It is possible that this difference is linked to the flattened turbulent mean velocity profile in the shear-turbulence state, which leads to a lift-up process more localised towards the near-wall region.
Normalised Nusselt numbers for flows at
$C=1, {\textit{Re}}=3000$
subject to optimal forces in the full and reduced rotational symmetries.
$(\textit a)$
Comparison between
$O_{Fm}$
(solid lines) and
$O_m$
(dashed lines).
$(\textit b)$
Comparison between
$O_{Fm}$
optimised for shear turbulence (solid line) and a steady two-dimensional laminar state (dashed line).
Optimal forces at
$A_0=10^{-6}, C=1$
. The largest arrow has magnitude
$6.01 \times 10^{-4}$
in
$(\textit {a})$
,
$6.60\times 10^{-4}$
in
$(\textit {b})$
,
$5.50\times 10^{-4}$
in
$(\textit {c})$
,
$5.15\times 10^{-4}$
in
$(\textit {d})$
,
$5.76\times 10^{-4}$
in
$(\textit {e})$
,
$5.40\times 10^{-4}$
in
$(\textit{f})$
,
$6.34 \times 10^{-4}$
in
$(\textit {g})$
,
$5.55\times 10^{-4}$
in
$(\textit {h})$
,
$5.99\times 10^{-4}$
in
$(\textit {i})$
.
$N=1$
initial velocity condition is used for
$O_{Fm}$
and
$N=3$
initial velocity conditions are used for
$O_{m}$
. In the ‘steady’ case, axial independence is imposed on the velocity which results in a steady response velocity to the force.
3.4. Optimisation in the convective-turbulence regime
Here, optimisation was first considered for a weakly convective turbulent state at
$C=4$
. As the velocity fields in the convective state are time-dependent, optimisations were initially performed using several initial velocity fields (
$N\gt 1$
) at this
$C$
and random initial forces. However, even for small
$A_0=10^{-7}$
it was found that the flow is rapidly laminarised by the force, it approaches a two-dimensional streamwise-independent state with enhanced
$\textit {Nu}$
. Therefore in this case, like at
$C=3$
for the laminarisation regime,
$N=1$
is sufficient. Also like the laminarised case, the optimal Nusselt number shows substantial improvement compared with the unforced case, 50 % at
$A_0=10^{-7}$
, and the structure of the optimal, although calculated for
$O_1$
, looks very similar to that of figure 5(c), close to
$O_{2}$
symmetry. Due to the laminarisation by the new force at
$C=4$
, optimisations for
$O_{Fm}$
exhibit similar behaviour to those for the laminarised case at
$C=3$
.
Flows forced by
$O_1$
stabilised to larger
$C$
. Here
$A_0=10^{-7}$
. The vertical dashed line indicates the optimisation target time
$\mathcal T=600$
.
We therefore optimise heat transfer at larger
$C$
to examine how far this laminarisation phenomenon occurs. Time series for
$\textit {Nu}(t)$
optimised at
$C=4-8, \mathcal T=600, A_0=10^{-7},$
for
$O_1$
are shown in figure 15. With increased buoyancy force, laminarisation by the force still occurs at
$C=5,6$
, but disappears at
$C=7, 8$
. Although the turbulence does not decay, the Nusselt number at
$C=7$
only fluctuates slightly. At
$C=8$
, the amplitude of fluctuations in
$\textit {Nu}(t)$
are similar to those of the unforced flow, but with a higher mean value.
Due to the laminarisation, the optimisations at
$C=5-7$
were found to be well converged, but as turbulence remained stronger for
$C=8$
, convergence was poor. Following the success of the short-time optimisation in the shear-turbulence regime, we performed short-
$\mathcal{T}$
optimisations at
$C=8$
. In this case, an intermediate value of
$\mathcal{T}=200$
was sufficient to achieve good convergence. Based on the observations for optimisations in the shear-turbulence regime, we consider optimisation in the reduced space
$O_{Fm}$
. Figure 16(a) shows
$\textit {Nu}(t)$
for
$O_{F4}$
as an example. The target time
$\mathcal{T}$
is not sufficient to capture the statistics of the end state, but is sufficient to capture the initial response to the force. With an increase in force amplitude, the Nusselt number gradually increases. Figure 16(b) compares
$O_{Fm}$
optimised for the unsteady convective turbulence versus optimisation for the artificially stabilised steady two-dimensional laminar state (setting
$K=1$
as before). Unlike for shear turbulence, this time they show little difference, and the optimal forces are close in structure, similar to figure 14(g–i) – when the flow is not too chaotic, including time dependence in the optimisation does not improve the Nusselt number further.
(a) Instantaneous Nusselt number for flow forced by
$O_{F4}$
at several
$A_0$
at
${\textit{Re}}=3000, C=8$
. The vertical dashed line indicates the optimisation target time
$\mathcal T=200$
. (b) Comparison between
$O_{Fm}$
optimised with unsteady convective turbulence (solid line) and with a two-dimensional steady velocity (dashed line).
Towards more chaotic convective turbulence, a further optimisation was carried out at
$C=16$
. Due to the stronger chaos, the target time was again reduced to
$\mathcal{T}=100$
. Comparison between
$O_{Fm}$
optimised with the short
$\mathcal{T}$
(solid line) and a steady two-dimensional laminar state (dashed line) is shown in figure 17(a). Improvement for the short-
$\mathcal{T}$
optimisation starts to be seen at larger
$A_0=10^{-5}$
. The roll structures of the force optimised with the short
$\mathcal{T}$
at all amplitudes are found to be located closer to the wall than for the corresponding steady calculation.
$O_{F5}$
at
$A_0=10^{-5}$
for the two optimisations are shown in figure 17(b,c). Such roll structures improve
$\textit {Nu}$
at a larger force amplitude, similar to optimisation in the shear-turbulence regime. Similar conclusions were drawn at
$C=32$
.
Further optimisations were also performed at
${\textit{Re}}=5000$
,
$C=1$
and
$40$
, corresponding to the strong shear-turbulence and convective-turbulence regimes (see figure 1). A short target time
$\mathcal{T}=40$
was necessary for optimisations with unsteady flows in both cases, due to very strong chaos. Comparison between
$\textit {Nu}$
for
$O_{Fm}$
at
$A_0=10^{-5}$
optimised with the short
$\mathcal{T}$
(solid line) and when imposing a steady two-dimensional laminar state during the optimisation (dashed line) are shown in figure 17(d). The inclusion of time dependence leads to an improvement in
$\textit {Nu}$
in both cases at
$C=1$
and
$C=40$
. The roll structures of the force optimised using unsteady flow are found to be closer to the wall than when steady flow is imposed via suppression of streamwise dependence (figure 17
e, f).
In this section we have found that the added force can laminarise convective turbulence, where it returns to a streamwise independent state. Pushing
${\textit{Re}}$
and
$C$
to larger values, such that the flow is chaotic, further enhancement in
$\textit {Nu}$
is found when time dependence is included in the optimisation, and rolls are again found to be located closer to the wall.
Optimisation at larger parameters. Top row:
$C=16$
,
${\textit{Re}}=3000$
,
$A_0$
in key to (a),
$O_{F5}$
at
$A_0=10^{-5}$
for steady and unsteady flow in (b) and (c). Bottom row:
${\textit{Re}}=5000$
,
$A_0=10^{-5}$
,
$C$
in key to (d),
$O_{F6}$
at
$C=1$
for steady and unsteady flow in (e) and (f). In (a,d), solid line is for the turbulent state, and dashed line indicates a force optimised with a steady state. The largest arrow has magnitude
$1.4 \times 10^{-3}$
in
$(\textit {b})$
,
$1.9 \times 10^{-3}$
in
$(\textit {c})$
,
$1.7 \times 10^{-3}$
in
$(\textit {e})$
and
$2.4 \times 10^{-3}$
in
$(\textit {f})$
.
4. Conclusions
In this work we have developed a heat-transfer optimisation method, based on a variational technique (Pringle & Kerswell Reference Pringle and Kerswell2010; Marensi et al. Reference Marensi, Ding, Willis and Kerswell2020), designed to identify the optimal body force that maximises heat transfer, in particular, subject to the time-dependent governing equations for the flow and subject to limited amplitude of the force. Focussing primarily on the feasibility and practicality of the method, optimisations have been conducted mostly at
${\textit{Re}}=3000$
with the constant temperature boundary condition, but the method has been applied to flows initially in each of the typical states of vertical heated pipe flow, i.e. shear-turbulence (
$C=1$
) laminarised flow (
$C=3$
), and convective turbulence (
$C\geqslant 4)$
.
Our optimisations confirm that the optimal body force is predominantly characterised by near-wall vortex structures that are uniform in the streamwise direction. This is consistent with optimised steady velocity fields that have been computed, e.g. Meng et al. (Reference Meng, Liang and Li2005). Direct comparison with previous results is not straightforward, however, since the optimisations were for the velocity field, rather than for a force, and the velocity field was independent of a buoyancy term (i.e. independent of
$C$
). Although such optimisations were not subject to a momentum equation for the velocity, the velocity was found to satisfy an equation similar to the Navier–Stokes equation with the addition of the ‘synergy force’, of the form
$\boldsymbol{F}_s=c\,A\,\boldsymbol{\nabla } T+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}$
, in which
$\boldsymbol{u}$
is time-independent. The scalar field
$A(\boldsymbol{x})$
is a Lagrange multiplier for the heat equation, and the constant
$c$
is a free parameter that is inversely proportional to a Lagrange multiplier of the viscous dissipation of the velocity field. The forces added to the Navier–Stokes in the present case are
$(4/{\textit{Re}})C\,T\hat {\boldsymbol{z}}+\boldsymbol{F}$
. Whilst we have made comparisons with optimisations we have obtained for
$\boldsymbol{F}$
with a steady response velocity, it should be cautioned that this is not directly equivalent to comparing with
$\boldsymbol{F}_s$
, the latter of which may not exist if subject to limited magnitude, as only
$c$
can be varied. Viscous dissipation is not considered in our calculations, while unrestricted dissipation corresponds to
$c\to \infty$
in the synergy approach. The novelty and advantage of our approach is that
$\boldsymbol{u}$
is subject to the (steady/unsteady) governing equations, where
$\boldsymbol{u}$
may be significantly affected by the flow regime determined by
$C$
, an optimal
$\boldsymbol{F}$
exists for arbitrarily small amplitude
$A_0$
, which changes with variation in
$A_0$
, and the target time
$\mathcal{T}$
is also an extra parameter in our calculations.
Forces optimised with different target times were found to exhibit different rotational symmetry. Specifically, the short-time optimal forces corresponded to larger azimuthal wavenumbers, while for
$\mathcal{T}\gtrsim 400$
the optimal force is unique with smaller azimuthal wavenumber. This flow pattern aligns with linear optimal perturbations that aim to maximise flow perturbation growth in isothermal flow (Schmid Reference Schmid2007). Forcing such modes efficiently modifies the flow. An implication of this result is that if there is limited time to impel a change in the velocity field, or if we link reduced response time with passage through a short section of pipe, then larger
$m$
might be favoured.
In the laminarisation regime (
$C=3$
), heat transfer increases with force amplitude
$A_0$
, as expected, then increases significantly at the point at which time-dependent flow is triggered by the force. To reduce computational costs considerably, we opted to compare the Nusselt number of flows forced by optimals constrained in the Fourier space,
$O_m$
and
$O_{Fm}$
(the former keeping modes
$0,m,2m,\ldots$
and the latter only keeping modes
$0$
and
$m$
). This method is found to be efficient and permits examination upon dependence on
$m$
. (Optimisation in the full space
$O_1$
leads to an
$O_m$
with
$m$
unspecified.) For very small
$A_0$
, rotational symmetry
$m=1$
is favoured, then
$m$
increases with
$A_0$
. However, rather than the enhancement in
$\textit {Nu}$
spiking for a particular
$m$
, the enhancement is similar over a broad range in
$m$
. Also for the laminarisation regime, there is a jump up in
$\textit {Nu}$
when turbulence is first triggered, but for larger
$A_0$
turbulence does not necessarily lead to an increase in
$\textit {Nu}$
. Visualisations of heat flux reveal that disorganisation caused by unsteadiness of the rolls inhibits heat transfer. An implication if this result is that an approach that seeks to simply trigger turbulence may not be efficient in enhancing
$\textit {Nu}$
. Organised rolls may be better than turbulent rolls, depending on the amplitude of roll that can be achieved.
Optimisation in the shear-turbulence regime (
$C=1$
) is challenging, as the flow is highly time-dependent and chaotic, preventing the long target times
$\mathcal{T}\gtrsim 400$
used in the laminar case. However, the method is found to still be effective for much shorter times, with
$\mathcal T = 50$
. Optimisations for unsteady and steady flow have been compared. (Short
$\mathcal{T}$
for the former, for the latter, imposing axial independence results in a flow that is steady in response to the force that is captured with a long
$\mathcal{T}$
.) It is found that including time dependence results in a force with rolls located closer to the wall and leads to flows with greater
$\textit {Nu}$
. Similarly, in the more chaotic convective state at
$C=16\,,32; {\textit{Re}}=3000$
, and the flow at
$C=1 \,, 40 \,; {\textit{Re}}=5000$
, optimisations with short
$\mathcal{T}$
show roll structures closer to the wall. Therefore, the location of the vortex is another consideration for potential application.
Energy consumption of flow forced by
$O_{F5}$
at
${\textit{Re}}=3000$
and
$O_{F6}$
at
${\textit{Re}}=5000$
, with
$A_0=10^{-5}$
and several
$C$
. Here,
$\dot {W}_P=\langle ({4}/{{\textit{Re}}})(1+{\beta ^{\prime }})u_{tot,z} \rangle$
is the pumping power,
$\dot {W}_F=\langle |\boldsymbol{F}\boldsymbol{\cdot }\boldsymbol{u}_{\textit{tot}}| \rangle$
is the rate of work done by the optimal force. For the ‘Steady’ case, the flow is restricted to be axially independent during the optimisation, which leads to steady flow and a ‘predicted’
$\textit {Nu}$
. The force is then applied in a DNS to calculate the ‘observed’
$\textit {Nu}$
. For the ‘Unsteady’ case, the flow is three-dimensional and time-dependent in the optimisation. The PEC is defined as
$\textrm{PEC}=({Nu/Nu_{F=0}})/{((\dot {W}_P+\dot {W}_F)/\dot {W}_{P,F=0})^{(1/3)}}$
. Flow is laminarised (steady) for the case marked ‘n/a’.
While it is acknowledged that accurately inducing the desired flow in practice is challenging, optimisations under the laminar steady-state assumption, such as those by Meng et al. (Reference Meng, Liang and Li2005), have inspired designs like the alternating elliptical axis tube, discrete double-inclined ribs tubes (Li et al. Reference Li, Meng and Guo2009) and many other applications (Liu & Sakr Reference Liu and Sakr2013; Sheikholeslami et al. Reference Sheikholeslami, Gorji-Bandpy and Ganji2015). Our results suggest that additional factors are worthy of consideration – the distance of the rolls from the wall affects the heat transfer, turbulence can both enhance or disrupt heat transfer. These will also be affected by the flow regime, here influenced by the alignment of buoyancy with the vertical pipe axis. Response time is another parameter, which may affect the choice of
$m$
in this case.
An important consideration that we have not directly included in our optimisation is the associated pumping power required for the flow, nor the power expended by the force. For our calculations, the proportional increase in wall friction and
$\textit {Nu}$
are similar, but it should be noted that simply increasing the pressure gradient does not lead to a proportional increase in heat transfer. Indeed, for laminar flow,
$\textit {Nu}$
is independent of
${\textit{Re}}$
(Su & Chung Reference Su and Chung2000). Table 2 shows the energy consumption and a performance evaluation criterion (PEC) close to the popular measure developed by Webb & Eckert (Reference Webb and Eckert1972), used in Meng et al. (Reference Meng, Liang and Li2005) and the more recent review of Li, Kadam & Yu (Reference Li, Kadam and Yu2023). We have replaced the friction ratio measured by
$\dot {W}_P/\dot {W}_{P,F=0}$
with
$(\dot {W}_P+\dot {W}_F)/\dot {W}_{P,F=0}$
in order to include the work done by the added force. The calculation using steady flow overestimates the enhancement, by 11 % and 31 % for the relatively quiescent cases at
${\textit{Re}}=3000,\,C=3$
and
$C=8$
, and more substantially when the flow is chaotic, by 55 % and 57 % at
$C=1$
(shear turbulence) and
$C=16$
(strong convection-driven turbulence). It is overestimated by 126 % for the case at
${\textit{Re}}=5000,\,C=1$
.
In summary, including the full governing equations for the velocity provides a greater range of parameters that are found to influence the optimal force configuration. Further including time dependence in the optimisation is challenging, but can lead to further enhancement of the heat transfer and significantly improves reliability in prediction of the overall heat-transfer enhancement.
Acknowledgements
This work used the Cirrus UK National Tier-2 HPC Service at EPCC (http://www.cirrus.ac.uk) funded by the University of Edinburgh and EPSRC (EP/P020267/1).
Funding
S.C. acknowledges the funding from Sheffield–China Scholarships Council PhD Scholarship Programme (CSC no. 202106260029).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Scalar dissipation
We start from the dimensionless governing equation for temperature
\begin{eqnarray} \frac {\partial T }{\partial t}+({\boldsymbol{u}_{\textit{tot}}}\boldsymbol{\cdot } \boldsymbol{\nabla }) T &=&\frac {1}{{\textit{Re}}\,{\textit{Pr}}}{\nabla}^{2}T -{\boldsymbol{u}_{\textit{tot}}}\boldsymbol{\cdot }{\hat {\boldsymbol{z}}}\,a_{\textit{tot}}(t) \, . \end{eqnarray}
Multiplying (A1) by the temperature
$T$
and integrating in time from
$0$
to
$\mathcal{T}$
, we obtain a new balance equation for
$T^2$
\begin{eqnarray} \int _{0}^{\mathcal{T}} \frac {1}{2} \left ( \frac {\partial }{\partial t} T^2 + ({\boldsymbol{u}_{\textit{tot}}}\boldsymbol{\cdot } \boldsymbol{\nabla }) T^2 \right )\textrm{d}t & = & \int _{0}^{\mathcal{T}} \left ( \frac {T}{{\textit{Re}}\,{\textit{Pr}}}{\nabla}^{2}T - T\,{\boldsymbol{u}_{\textit{tot}}}\boldsymbol{\cdot }{\hat {\boldsymbol{z}}}\,a_{\textit{tot}}(t)\right )\textrm{d}t \, . \end{eqnarray}
We aim for a statistically steady state, and therefore drop the first term on the left hand side. Integrating (A2) over the entire domain, the second term also vanishes, as it transports
$T^2$
but conserves the total. Applying integration by parts on the first term on the right-hand side leads to
\begin{eqnarray} 0 &\approx & \int _{0}^{\mathcal{T}} \left ( \frac {1}{{\textit{Re}}\, {\textit{Pr}}}\langle \boldsymbol{\nabla } \boldsymbol{\cdot }(T\boldsymbol{\nabla } T) - (\boldsymbol{\nabla } T)^2 \rangle - \langle T\,{\boldsymbol{u}_{\textit{tot}}}\boldsymbol{\cdot }{\hat {\boldsymbol{z}}}\,a_{\textit{tot}}(t) \rangle \right )\textrm{d}t \, . \end{eqnarray}
Note the ‘dissipation’ term involving
$\langle (\boldsymbol{\nabla } T)^2\rangle$
can only lead to a reduction in
$T^2$
. Pulling this term across and applying the divergence theorem gives
\begin{eqnarray} \int _{0}^{\mathcal{T}} \frac {1}{{\textit{Re}}\, {\textit{Pr}}}\langle (\boldsymbol{\nabla } T)^2\rangle \, \textrm{d}t \approx \int _{0}^{\mathcal{T}} \left ( \oint \frac {1}{{\textit{Re}}\, {\textit{Pr}}}T\,\boldsymbol{\nabla }T\boldsymbol{\cdot }\boldsymbol{\textrm{d}S} - \langle T\,{\boldsymbol{u}_{\textit{tot}}}\boldsymbol{\cdot }{\hat {\boldsymbol{z}}}\,a_{\textit{tot}}(t) \rangle \right ) \, \textrm{d}t \,. \end{eqnarray}
The surface integral term can be rewritten, given that
$T$
is constant at the wall (and applying periodicity at the ends)
\begin{eqnarray} \oint \frac {1}{{\textit{Re}}\, {\textit{Pr}}}T\,\boldsymbol{\nabla }T\boldsymbol{\cdot }\boldsymbol{\textrm{d}S} \,=\, T|_{r=1}\,\int \frac {1}{{\textit{Re}}\, {\textit{Pr}}}\boldsymbol{\nabla }T\boldsymbol{\cdot }\boldsymbol{\textrm{d}S}&=&T|_{r=1}\,Q_{\textit{in}} , \end{eqnarray}
where
$Q_{\textit{in}}$
is the total heat flux entering at the wall. For the heat sink term in (A4) we make the further approximation that
$\langle T\,{\boldsymbol{u}_{\textit{tot}}}\boldsymbol{\cdot }{\hat {\boldsymbol{z}}}\,a_{\textit{tot}}(t)\rangle \approx T_b \, \langle {\boldsymbol{u}_{\textit{tot}}}\boldsymbol{\cdot }{\hat {\boldsymbol{z}}}\,a_{\textit{tot}}(t)\rangle = T_b \,Q_{\textit{out}}$
, where
$Q_{\textit{out}}$
represents the total net rate at which heat is carried out of the domain. (Observe that, for a spatially uniform heat sink
$\epsilon$
, as in Marensi et al. (Reference Marensi, He and Willis2021), no approximation is necessary:
$\langle T\,\epsilon (t)\rangle = T_b \, \langle \epsilon (t) \rangle = T_b \,Q_{\textit{out}}$
. While that models a temporal rather than axial change in temperature, very little difference is observed, see figure 1(b). Note that the velocity profile tends to be ‘flatter’, i.e. providing a more uniform heat sink, when the flow is turbulent.) Assuming that on average
$Q_{\textit{in}}=Q_{\textit{out}}$
, (A4) becomes
\begin{eqnarray} \int _{0}^{\mathcal{T}} \frac {1}{{\textit{Re}}\, {\textit{Pr}}}\langle (\boldsymbol{\nabla } T)^2\rangle \, \textrm{d}t\approx \int _{0}^{\mathcal{T}} (T|_{r=1}-T_b) \, Q_{\textit{in}} \, \textrm{d}t\,. \end{eqnarray}
This equation implies the maximisation (minimisation) principle for
$\langle (\boldsymbol{\nabla } T)^2\rangle$
– when the temperature difference is fixed, maximising the left-hand side implies the largest
$Q_{\textit{in}}$
, and hence the largest
$\textit {Nu}$
. (When
$Q_{\textit{in}}$
is fixed, minimising the left-hand side implies the smallest temperature difference, and hence the largest
$\textit {Nu}$
.)
