2.1. Waveform parameterisation

The systematic exploration of non-sinusoidal waveforms necessitates a parameterisation methodology capable of representing arbitrary temporal structures whilst maintaining mathematical tractability for gradient-based optimisation. To this end, a novel parameterisation based on Lagrangian-polynomial interpolation is introduced. The spanwise wall velocity signal $w_{\textit{w}}(t)$ is constructed over one complete period $T$ subject to two fundamental constraints: periodicity and anti-symmetry. The anti-symmetry constraint, expressed as $w_{\textit{w}}(t + T/2) = -w_{\textit{w}}(t)$ , guarantees symmetric positive and negative phases. This ensures zero net momentum injection over each complete cycle, a requisite for achieving equilibrium states, which closed-loop strategies often find challenging to attain (Cavallazzi et al. Reference Cavallazzi, Guastoni, Vinuesa and Pinelli2024). Furthermore, the restriction to time-continuous signals facilitates practical implementation compared with discontinuous waveforms.

This mathematical structure leads to the unique definition of the complete waveform via control points specified exclusively over the first half-period, $t \in [0, T/2]$ . The five control points (the three free parameters plus the two boundaries) are positioned at equally spaced temporal intervals $t_i^1 = (i-1)T/8$ for $i \in \{1,\ldots ,5\}$ , as illustrated in figure 1. To ensure continuity at period boundaries, the first and last control points are constrained to zero velocity

(2.1) \begin{align} y_1^1 = 0 \quad (\textrm{at} t=0) \quad \textrm{and} \quad y_5^1 = 0 \quad (\textrm{at } t=T/2). \end{align}

Therefore, the waveform topology is determined entirely by the three interior control point amplitudes ( $y_2^1, y_3^1, y_4^1$ ). The selection of Lagrangian-polynomial interpolation, rather than alternative techniques such as splines, is motivated by its computational efficiency for a small number of evenly spaced control points. The polynomial passing exactly through the design points is constructed directly, thereby circumventing the need to solve complex systems of equations. Furthermore, with only five control points, the primary limitations of high-degree interpolation (such as Runge’s phenomenon) are effectively avoided.

Figure 1.

Illustration of the Lagrangian-polynomial parameterisation methodology. The boundary control points are fixed at zero amplitude to ensure periodicity, leaving three free interior control points ( $y_2^1, y_3^1, y_4^1$ , red crosses) that define the waveform shape over $[0, T/2]$ . The second half-period is obtained through anti-symmetric reflection.

Having established the control point structure and interpolation methodology, attention is now directed towards the parameter mapping strategy itself. To establish a robust parameter mapping whilst ensuring injectivity (whereby each parameter pair $(a_1, a_2)$ corresponds to a unique waveform topology, as detailed in Appendix A), the optimisation process avoids direct manipulation of the interior control point amplitudes. Instead, a reduced parameterisation is adopted, relying on only two parameters, $a_1 \in [0,1]$ and $a_2 \in [0,1]$ . These parameters fully define the normalised waveform shape (prior to amplitude scaling) through determination of the three interior control point amplitudes according to the following transformation:

(2.2) \begin{equation} \begin{cases} y_1^1 = 0, \quad y_5^1 = 0, \\[4pt] y_2^1 = a_1 ,\\[4pt] y_3^1 = a_{2} (1 - a_1) ,\\[4pt] y_4^1 = (1 - a_1) (1 - a_2). \end{cases} \end{equation}

Worth highlighting is that this formulation explicitly enforces the normalisation constraint $\sum _{i=1}^5 y_i^1 = 1$ . As demonstrated rigorously in Appendix A, this constraint plays a fundamental role in guaranteeing the injectivity of the mapping from parameter pairs $(a_1, a_2)$ to waveform topologies. Injectivity, in turn, represents a cornerstone requirement for the effective deployment of gradient-based optimisation algorithms, fulfilling three objectives.

1. (i) Bijective mapping: it ensures that the parameter space defined by $T^+$ , $W^+$ and $(a_1, a_2)$ is bijectively mapped to the space of admissible waveforms, establishing a well-defined correspondence between parameters and flow-control signals.

2. (ii) Elimination of redundancy: it prevents distinct parameter combinations from leading identical performance metrics. Such redundancy would introduce ambiguity in gradient computation and potentially cause the algorithm to oscillate between equivalent parameter sets.

3. (iii) Convergence stability: it guarantees that the PBO algorithm (see § 2.2) can navigate the parameter space without encountering degenerate regions, thereby enabling systematic exploration and stable convergence towards optimal configurations.

The continuous waveform is constructed via fourth-order Lagrangian-polynomial interpolation. Specifically, for any time $t$ within the first half-period (denoted by superscript $k=1$ ), the wall velocity is computed as

(2.3) \begin{equation} P_1(t) = \sum _{i=1}^5 y_i^1 L_i^1(t), \quad \textrm{with} \quad L_i^1(t) = \prod _{\substack {j=1 \\ j \neq i}}^5 \frac {t - t_j^1}{t_i^1 - t_j^1}. \end{equation}

For the second half-period (superscript $k=2$ ), anti-symmetry is enforced by constructing control points with reflected temporal coordinates $t_i^2 = t_i^1 + T/2$ and amplitudes $y_i^2 = -y_{6-i}^1$ , yielding the polynomial $P_2(t)$ . The complete signal over one period is assembled as

(2.4) \begin{equation} w_{\textit{w}}(t) = \unicode {x1D7D9}_{t \leqslant T/2}(t) \, P_1(t) + \unicode {x1D7D9}_{t \gt T/2}(t) \, P_2(t), \end{equation}

where $\unicode {x1D7D9}$ denotes the indicator function. Finally, this pre-scaled waveform is multiplied by a scaling factor $\text{sc} = W^+_{\textit{max}}/\max _t(w_{\textit{w}}(t))$ to ensure the maximum amplitude attains the target value $W^+_{\textit{max}}$ . The complete parameterisation mapping $\mathcal{L}$ thus incorporates $W^+_{\textit{max}}$ as an additional degree of freedom controlling the absolute magnitudes of the interior control points $(y_2^1, y_3^1, y_4^1)$ , whilst $(a_1, a_2)$ govern their relative proportions. This mapping is formally expressed as

(2.5) \begin{equation} \mathcal{L} : \left(a_1, a_2, T^+, W^+_{\textit{max}}\right) \mapsto \text{sc} \times P_1(t). \end{equation}

This methodology enables the representation of diverse waveform topologies, ranging from sinusoidal configurations to plateau-like structures with extended constant-velocity phases, through systematic variation of the parameters $(a_1, a_2)$ and $W^+_{\textit{max}}$ , providing sufficient expressiveness to discover non-trivial temporal structures whilst maintaining a compact parameter space. Before application to waveform optimisation, the Lagrangian parameterisation was validated against high-fidelity DNS benchmarks for sinusoidal oscillations at $T^+=125$ and $W^+=18$ (Quadrio & Ricco Reference Quadrio and Ricco2004), demonstrating that the method reproduces known drag-reduction performance ( $P_{\textit{sav}} = 39.27\,\%$ versus DNS reference of $39.3\,\%$ ) within acceptable tolerances, as detailed in Appendix B.

2.2. Optimisation algorithm: policy-Based Optimisation

The optimisation is performed using the PBO algorithm (Viquerat et al. Reference Viquerat, Rabault, Kuhnle, Ghraieb, Larcher and Hachem2023), a gradient-based reinforcement learning method specifically designed for open-loop control in fluid dynamics applications. PBO is conceived as a hybrid framework that integrates the gradient-based learning efficiency of deep reinforcement learning with the population-based robustness of evolutionary strategies, drawing particular inspiration from the covariance matrix adaptation evolution strategy (CMA-ES) (Hansen Reference Hansen2023). This synthesis permits systematic exploration of high-dimensional parameter spaces whilst maintaining computational tractability through gradient-based updates. Benchmark comparisons have demonstrated that PBO performance is comparable to CMA-ES and substantially exceeds that of simplified single-step methods when minimising complex analytical functions (Viquerat et al. Reference Viquerat, Rabault, Kuhnle, Ghraieb, Larcher and Hachem2023). Furthermore, the algorithm has been successfully applied to sinusoidal SWO parameter optimisation in prior work (Guérin et al. Reference Guérin, Cordier, Flageul, Grieu and Agostini2025), establishing its suitability for turbulent flow-control applications.

In the context of SWO, the control problem is formulated as a single-step decision process. In contrast to closed-loop control strategies, wherein an agent continuously adjusts actions based on time-varying flow states, the present open-loop scenario necessitates the determination of a single optimal parameter vector (comprising $a_1$ , $a_2$ , $T^+$ , $W^+_{\textit{max}}$ ) that is applied throughout the simulation. Consequently, instantaneous state observations are rendered unnecessary, and the optimisation relies exclusively on the terminal reward (e.g. drag reduction or analogy factor) obtained once the flow has attained a statistically steady state.

The PBO agent maintains a stochastic control policy $\pi _{\boldsymbol{\theta }}$ , parameterised as a multivariate Gaussian distribution $\mathcal{N}(\boldsymbol{\mu }, \boldsymbol{\varSigma })$ . This distribution is defined by three distinct neural network components that generate:

1. (i) the mean vector ( $\boldsymbol{\mu }$ ), representing the current estimate of optimal parameters;

2. (ii) the standard deviation vector ( $\boldsymbol{\sigma }$ ), governing the magnitude of exploratory sampling; and

3. (iii) rhe off-diagonal covariance coefficients ( $\boldsymbol{\varSigma }$ ), enabling the algorithm to identify and exploit correlations between control parameters.

This explicit parameterisation of the covariance structure constitutes a principal advantage of PBO relative to simplified reinforcement learning algorithms, such as Proximal Policy Optimization-1 (Ghraieb et al. Reference Ghraieb, Viquerat, Larcher, Meliga and Hachem2021), as it permits the sampling distribution to align with the principal axes of the objective function’s topography. This alignment proves particularly beneficial in landscapes wherein parameters exhibit strong coupling, as is characteristic of turbulent flow-control problems where period and amplitude interactions govern performance.

At the algorithmic level, PBO proceeds as follows. At each generation, a population of $n_{\textit{ind}}$ individuals is sampled from the current policy distribution $\pi _{\boldsymbol{\theta }}$ and evaluated through independent LES. The resulting rewards are utilised to update the policy via gradient ascent on the loss function

(2.6) \begin{equation} L(\theta ) = \mathbb{E}_{a \sim \pi _\theta }[\log \pi _\theta (a)\, \hat {R}(a)], \quad \textrm{with} \quad \hat {R}(a) = \max \!\left (\frac {r(a) - \mu _r}{\sigma _r},\, 0\right ), \end{equation}

where $r(a)$ denotes the reward (e.g. $P_{\textit{sav}}$ or $\overline {A_n}$ ), and $\mu _r$ and $\sigma _r$ represent the reward mean and standard deviation over the current generation. This normalised, clipped formulation reduces variance and discards negative-advantage actions to stabilise the learning process. Three separate neural networks, each consisting of 3 hidden layers with 2 neurons, independently generate $\boldsymbol{\mu }$ , $\boldsymbol{\sigma }$ and the off-diagonal covariance coefficients. The covariance matrix is constructed via the hypersphere decomposition method (Viquerat et al. Reference Viquerat, Rabault, Kuhnle, Ghraieb, Larcher and Hachem2023), which guarantees that valid positive semi-definite matrices are always produced from the network outputs. Each network is trained independently using the ADAM optimiser (Kingma & Ba Reference Kingma and Ba2014), with network-specific learning rates, epoch counts, batch sizes and generation history depths, as summarised in table 1. An exponential decay parameter $\alpha$ is employed to weight recent generations more heavily than earlier ones, thereby balancing off-policy data re-use with responsiveness to the current reward landscape. Complete algorithmic details, including the pseudo-code, sensitivity analyses with respect to population size and network architecture, and benchmark comparisons against CMA-ES, are available in Viquerat et al. (Reference Viquerat, Rabault, Kuhnle, Ghraieb, Larcher and Hachem2023) and Guérin et al. (Reference Guérin, Cordier, Flageul, Grieu and Agostini2025).

Table 1.

The PBO training configuration for waveform optimisation. Here, $n_{\textit{ind}}$ : individuals per generation; $n_{\textit{para}}$ : parallel simulations; $n_{\textit{CPUs}}$ : CPU cores per simulation; $\alpha$ : action sampling parameter; $\lambda _r$ : learning rates for ( $\mu , \sigma , \varSigma$ ); $n_e$ : epochs; $n_b$ : mini-batch sizes; $n_g$ : generation history depth.

The optimisation procedure operates in discrete generations. At each iteration, a population of parameter configurations (termed individuals) is sampled from the current policy $\pi _{\boldsymbol{\theta }}$ and evaluated through independent LES. The collected rewards are subsequently employed to update the policy parameters via gradient ascent, thereby concentrating probability mass within high-performing regions of the parameter space. This population-based evaluation strategy renders PBO particularly well suited to high-performance computing environments, as the computational cost is dominated by parallel fluid evaluations rather than sequential sampling (an inherent advantage relative to methods such as Bayesian optimisation). The complete PBO training configuration is summarised in table 1. Relative to the sinusoidal optimisation study presented in Guérin et al. (Reference Guérin, Cordier, Flageul, Grieu and Agostini2025), wherein the parameter space comprised only two dimensions ( $T^+$ and $W^+$ ) with 8 individuals per generation, the present waveform optimisation involves a four-dimensional parameter space. This increased action-space dimensionality necessitated elevating the population size to $n_{\textit{ind}}=10$ to ensure adequate sampling density, as convergence quality is substantially influenced by the ratio of population size to parameter space dimensionality (Viquerat et al. Reference Viquerat, Rabault, Kuhnle, Ghraieb, Larcher and Hachem2023). The computational configuration was adapted to accommodate the expanded parameter space and available computational resources. The strategy employs serial execution of simulations ( $n_{\textit{para}}=1$ ), evaluating the 10 individuals sequentially rather than simultaneously, with each simulation allocated 256 CPU cores, i.e. representing a fourfold increase relative to the configuration considered in Guérin et al. (Reference Guérin, Cordier, Flageul, Grieu and Agostini2025). These modifications produce contrasting effects on computational efficiency: wall time per generation is reduced by more than 50 % to approximately 2.8 h owing to the elevated per-simulation core allocation, whilst total CPU-time consumption per generation increases to 720 CPU hours owing to serial execution. This configuration reflects a deliberate trade-off, wherein optimisation throughput (reduced wall time) is prioritised over computational resource efficiency (increased CPU time). This strategy proved essential for exploring the higher-dimensional waveform parameter space within practical time scales whilst achieving convergence typically within 40–72 generations depending on problem complexity.

The hyperparameters listed in table 1, with the exception of $n_{\textit{ind}}$ , correspond to the default values provided by Viquerat et al. (Reference Viquerat, Rabault, Kuhnle, Ghraieb, Larcher and Hachem2023), which were demonstrated to lead to robust convergence across benchmark problems of varying dimensionality (from two- to ten-dimensional Rosenbrock functions) without problem-specific tuning. The unique modification relative to the sinusoidal optimisation of Guérin et al. (Reference Guérin, Cordier, Flageul, Grieu and Agostini2025) is the increase in population size from 8 to 10, motivated by the expanded four-dimensional parameter space. This limited tuning requirement is a practical advantage for the broader applicability of the PBO-LES framework to other flow configurations. Regarding initialisation, the PBO algorithm is designed to mitigate sensitivity to the first guess. The initial policy distribution is deliberately centred within the imposed action boundaries with a large variance and a diagonal covariance matrix, thereby ensuring broad exploration during the early generations. As the optimisation progresses, the covariance matrix adapts to the local topology of the reward landscape through gradient-based updates of the correlation network. The population-based nature of PBO, i.e. wherein multiple individuals are evaluated per generation and the entire distribution is updated, provides an intrinsic mechanism for escaping suboptimal initial conditions, analogous to classical evolution strategies.

The choice of PBO over alternative optimisation techniques, such as gradient-based methods relying on finite-difference approximations, is motivated by several considerations. First, the inherent stochasticity of the reward signal (arising from finite-time averaging over limited oscillation cycles) renders finite-difference gradient estimates unreliable, as their accuracy depends on the difference between two noisy evaluations. The PBO, by contrast, extracts gradient information from the collective reward distribution across an entire population, effectively filtering statistical noise. Second, PBO naturally exploits parallel evaluation of independent simulations, avoiding the sequential bottleneck of perturbation-based gradient estimation. Third, the full covariance matrix adaptation provides directional search information analogous to second-order methods, without requiring explicit Hessian computation. Finally, the population-based exploration reduces the risk of convergence to local optima, which is particularly relevant for the DHT optimisation cases wherein the reward landscape is not known a priori.

2.3. Performance metrics

Having established the waveform parameterisation methodology and the optimisation algorithm, it is pertinent to define the performance metrics through which control efficacy is quantified. The assessment of SWO necessitates metrics that capture both the beneficial outcomes and the energetic costs inherent to active flow-control strategies. Throughout this investigation, five distinct metrics are employed, partitioned into two fundamental categories: momentum-transport metrics ( $P_{\textit{sav}}$ , $P_{\textit{req}}$ , $P_{\textit{net}}$ ) characterising drag reduction and energetic efficiency, and thermal-transport metrics ( ${\textit{Nu}}_r$ , $\overline {A_n}$ ) characterising heat-transfer enhancement and thermal-momentum correlation, respectively. Throughout this section, the overbar notation $\overline {(\boldsymbol{\cdot })}$ denotes averaging over time and the homogeneous spatial directions ( $x$ and $z$ ), whilst $\langle \boldsymbol{\cdot }\rangle _{x,z}$ represents spatial averaging over the streamwise and spanwise directions only.

The drag-reduction metric ( $P_{\textit{sav}}$ ) quantifies the relative reduction in skin-friction coefficient achieved through SWO, expressed as a percentage

(2.7) \begin{equation} P_{\textit{sav}} = 1-\frac {\overline {C_{\!f,\textit{ac}}}}{\overline {C_{\!f,\textit{unac}}}}. \end{equation}

The skin-friction coefficient is defined as $C_{\!f} = \tau _w / (({1}/{2}) \rho U_b^2) = 2\nu U_b^{-2} \, \partial \langle u \rangle _{x,z} / \partial y |_{\textit{w}}$ , where $\tau _w$ denotes the wall shear stress. Positive values of $P_{\textit{sav}}$ indicate successful attenuation of wall shear stress, with the magnitude reflecting the extent to which turbulent momentum transport has been suppressed relative to the uncontrolled baseline. This metric serves as the primary objective for validation cases wherein drag reduction is maximised independently of energetic expenditure considerations.

The control-cost metric ( $P_{\textit{req}}$ ) represents the power required to actuate the oscillating walls, normalised by the frictional power dissipation in the unactuated channel (Quadrio & Ricco Reference Quadrio and Ricco2004)

(2.8) \begin{equation} P_{\textit{req}} = \frac {\overline {\tau _{w,z} w_{\textit{w}}}_{\textit{ac}}}{U_b \, \overline {\tau _{w,x}}_{\textit{unac}}}, \end{equation}

where $U_b$ denotes the bulk velocity (maintained constant through the constant flow rate (CFR) strategy) and $\tau _{w,x}=\nu ({\partial \langle u \rangle _{x,z}}/{\partial y})$ and $\tau _{w,z}=\nu ({\partial \langle w \rangle _{x,z}}/{\partial y})$ denote the streamwise and spanwise wall shear stresses, respectively, with summation accounting for contributions from both channel walls. This quantity is invariably positive and increases monotonically with both oscillation amplitude and frequency, thereby constituting the energetic penalty inherent to active flow-control strategies. Whilst drag reduction quantifies the beneficial reduction in pumping power, this metric captures the required investment of energy to maintain the control action itself.

The net energetic gain metric ( $P_{\textit{net}}$ ) represents the comprehensive energetic balance achieved through flow control, reconciling the competing contributions of drag reduction and control cost

(2.9) \begin{equation} P_{\textit{net}} = P_{\textit{sav}} - P_{\textit{req}}. \end{equation}

Positive values indicate that the reduction in pumping power required to maintain CFR exceeds the power expended for wall actuation, thereby rendering the control strategy energetically advantageous for practical applications. This metric proves particularly pertinent for engineering implementations wherein energy efficiency constitutes a primary design constraint, serving as the objective function for validation cases addressing energetically optimal control configurations. The maximisation of $P_{\textit{net}}$ necessitates a delicate balance between achieving substantial drag reduction whilst simultaneously minimising the energetic investment required to produce the oscillatory wall motion.

The heat-transfer enhancement metric ( ${\textit{Nu}}_r$ ) quantifies the relative augmentation of convective heat transfer induced by SWO. The Nusselt number $\textit{Nu}$ is defined as ${\textit{Nu}} = 2H h / \lambda$ , where $h$ denotes the convective heat-transfer coefficient and $\lambda$ the thermal conductivity; in the present dimensionless formulation, this reduces to ${\textit{Nu}} = 2 \, \partial \langle \varTheta \rangle _{x,z} / \partial y |_{\textit{w}}$ , evaluated at the wall. The relative enhancement is thus expressed as

(2.10) \begin{equation} {\textit{Nu}}_r = \frac {\overline {{\textit{Nu}}_{\textit{ac}}}}{\overline {{\textit{Nu}}_{\textit{unac}}}}. \end{equation}

Values exceeding unity indicate enhanced thermal transport relative to the uncontrolled baseline; however, this metric provides no information regarding the associated momentum-transport penalties. It is pertinent to recognise that simultaneous increases in both heat transfer and drag would yield elevated ${\textit{Nu}}_r$ values without constituting desirable performance for practical applications, as thermal enhancement achieved through proportional momentum-transport augmentation offers no engineering advantage. Consequently, the evaluation of DHT performance necessitates a metric that explicitly accounts for the relative magnitudes of thermal- and momentum-transport modifications.

The analogy factor ( $\overline {A_n}$ ) represents the ratio of relative heat-transfer enhancement to relative drag increase, thereby quantifying the extent to which the Reynolds analogy is broken by the control action

(2.11) \begin{equation} \overline {A_n} = \frac {{\textit{Nu}}_r}{C_{\!f,r}} = \frac {\overline {{\textit{Nu}}_{\textit{ac}}}/\overline {{\textit{Nu}}_{\textit{unac}}}}{\overline {C_{\!f,\textit{ac}}}/\overline {C_{\!f,\textit{unac}}}}. \end{equation}

This metric constitutes the primary objective for DHT optimisation, with $\overline {A_n} \gt 1$ indicating preferential enhancement of thermal transport relative to momentum transport. The achievement of $\overline {A_n} \gt 1$ demonstrates that thermal mixing has been augmented to a greater extent than momentum mixing, thereby breaking the classical Reynolds analogy which traditionally couples these transport processes under dynamically similar conditions. Higher values of $\overline {A_n}$ signify more effective decoupling of thermal and momentum transport, representing superior performance for applications wherein heat-transfer enhancement is desired without commensurate drag penalties. The maximisation of this metric thus constitutes a fundamentally different optimisation objective compared with drag reduction or net energetic gain, as it prioritises the selective manipulation of thermal structures whilst minimising momentum-transport modifications.

Mixed boundary conditions (MBC) are employed for the thermal field (Kasagi, Tomita & Kuroda Reference Kasagi, Tomita and Kuroda1992), imposing zero temperature at both walls. The unsteady source term associated with MBC ensures that a constant heat flux enters the domain. This configuration provides a robust baseline Reynolds analogy under uncontrolled conditions, thereby rendering the achievement of DHT ( $\overline {A_n} \gt 1$ ) non-trivial, as both momentum and thermal transport are strongly coupled through the turbulent flow structures. The challenge of breaking this analogy under MBC conditions provides a stringent test for the efficacy of SWO in selectively manipulating thermal-transport mechanisms independently of momentum-transport processes.

The thermal-energy efficiency metric quantifies the thermal enhancement per unit of control energy expended, providing a comprehensive measure of the energetic efficiency with which thermal transport is augmented through active forcing. This metric is expressed as the ratio of relative heat-transfer enhancement to the absolute magnitude of required control power

(2.12) \begin{equation} \textrm{Thermal-energy efficiency} = \frac {\Delta {\textit{Nu}}}{|P_{\textit{req}}|}, \end{equation}

where $\Delta {\textit{Nu}} = {\textit{Nu}}_r - 1$ represents the fractional heat-transfer enhancement relative to the unactuated baseline. This formulation enables direct comparison of thermal performance across configurations operating at different actuation amplitudes and periods, as it explicitly accounts for the energetic investment required to achieve the observed thermal-transport augmentation. Higher values indicate superior thermal benefit per unit of control energy, thereby identifying actuation strategies that maximise heat-transfer enhancement whilst minimising energetic expenditure. This metric proves particularly pertinent for engineering applications wherein both thermal performance and energy efficiency constitute paramount design constraints, such as CSP systems requiring enhanced heat transfer without excessive parasitic power consumption.

2.4. Numerical simulation framework and grid configurations

The optimisation framework described in the preceding sections necessitates the execution of hundreds of flow simulations to systematically explore the waveform parameter space. To render this approach computationally tractable whilst maintaining physical fidelity, a multi-fidelity numerical strategy is employed, combining computationally efficient implicit LES for optimisation with high-fidelity DNSs for validation.

All simulations are performed using the Xcompact3d framework (Bartholomew et al. Reference Bartholomew, Deskos, Frantz, Schuch, Lamballais and Laizet2020), wherein the incompressible Navier–Stokes equations are solved in a turbulent channel flow at friction Reynolds number ${\textit{Re}}_\tau =200$ . For DHT investigations, the momentum equations are coupled with a passive scalar transport equation for temperature at unity Prandtl number ( ${\textit{Pr}}=1$ ). The numerical methodology employs sixth-order compact finite-difference schemes for spatial derivatives combined with a third-order Runge–Kutta scheme for temporal advancement, ensuring high accuracy whilst maintaining computational efficiency. For further details on the numerical configurations and governing continuous equations, refer to Guérin et al. (Reference Guérin, Flageul, Cordier, Grieu and Agostini2024) for the DNS set-up and Guérin et al. (Reference Guérin, Cordier, Flageul, Grieu and Agostini2025) for the LES configuration.

The dimensionless governing equations, using the channel half-height $H$ and the bulk velocity $U_b$ as reference scales, are expressed as

(2.13) \begin{equation} \begin{cases} \boldsymbol{\nabla} ^{\star } \boldsymbol{\cdot }\boldsymbol{u}^{\star } = 0, \\[3pt] \displaystyle \frac {\partial \boldsymbol{u}^{\star }}{\partial t^{\star }} + \left (\boldsymbol{u}^{\star } \boldsymbol{\cdot }\boldsymbol{\nabla} ^{\star }\right )\boldsymbol{u}^{\star } = -\boldsymbol{\nabla} ^{\star } p^{\star } + \frac {1}{\textit{Re}}\,\Delta ^{\star }\boldsymbol{u}^{\star } + \boldsymbol{f}^{\star }, \\[8pt] \displaystyle \frac {\partial \varTheta ^{\star }}{\partial t^{\star }} + \boldsymbol{u}^{\star } \boldsymbol{\cdot }\boldsymbol{\nabla} ^{\star }\varTheta ^{\star } + A^{\star }u^{\star } = \frac {1}{\textit{Pe}}\,\Delta ^{\star }\varTheta ^{\star }, \end{cases} \end{equation}

where $\boldsymbol{u}^{\star } = \boldsymbol{u}/U_b$ , $p^{\star } = p/U_b^2$ , $t^{\star } = tU_b/H$ , $\textit{Re} = U_b H / \nu$ is the bulk Reynolds number, $\textit{Pe} = \textit{Re}\,\textit{Pr}$ the Péclet number and $\boldsymbol{f}^{\star }$ the body force maintaining CFR. The temperature $\varTheta$ is non-dimensionalised following the MBC formulation of Kasagi et al. (Reference Kasagi, Tomita and Kuroda1992), with $\varTheta = (T - \langle T_{\textit{w}} \rangle _t)/(T_b - \langle T_{\textit{w}} \rangle _t)$ , and $A^{\star }$ represents the streamwise temperature gradient arising from the constant wall-heat-flux condition. Periodic boundary conditions are imposed in the streamwise ( $x$ ) and spanwise ( $z$ ) directions, with no-slip conditions for velocity and $\varTheta _{\textit{w}} = 0$ at both channel walls. The flow is driven at CFR through an adaptive streamwise body force. The divergence-free condition is enforced through a fractional step method, wherein a Poisson equation for the pressure is solved using three-dimensional fast Fourier transforms (Laizet & Lamballais Reference Laizet and Lamballais2009). Hereafter, the star superscripts are dropped for brevity.

The investigation employs two distinct grid configurations, each designed to fulfil specific computational requirements within the optimisation-validation workflow. Table 2 summarises the complete specifications for all configurations, with resolution requirements motivated by the physical scales that must be adequately captured for each problem category.

Table 2.

Grid specifications for the numerical configurations employed in the present investigation, at ${\textit{Re}}_\tau \approx 200$ and ${\textit{Pr}}=1$ . Here, $S_1$ is an implicit LES configuration used for all optimisation workflows (drag reduction and DHT), whilst $S_2$ provides DNS-level resolution for final validation.

The $S_1$ configuration is employed for all optimisation campaigns, encompassing drag reduction, net energy gain and DHT investigations. This configuration comprises a computational domain of $(L_x, L_y, L_z) = (8, 2, 6)$ discretised with $128 \times 65 \times 96$ grid points, leading to streamwise and spanwise resolutions of $\Delta x^+ = 11.0$ and $\Delta z^+ = 11.0$ , respectively. The wall-normal resolution varies from $\Delta y^+_{\textit{min}} = 2.0$ at the wall to $\Delta y^+_{\textit{max}} = 15.7$ at the channel centreline, ensuring adequate representation of near-wall turbulent structures. The relatively compact streamwise domain proves sufficient for parametric studies, as the SWO primarily influence the near-wall region wherein turbulent streaks and quasi-streamwise vortices develop over shorter streamwise distances. This configuration balances computational efficiency with adequate physical fidelity for both momentum- and thermal-transport phenomena, requiring approximately 96 CPU hours per simulation. This set-up employs a fourth-order accurate spectral vanishing viscosity (Lamballais, Fortuné & Laizet Reference Lamballais, Fortuné and Laizet2011) to provide controlled numerical dissipation at the smallest resolved scales, thereby characterising these simulations as implicit LES. This approach avoids explicit subgrid-scale modelling whilst maintaining numerical stability, a strategy that has been rigorously validated against DNS benchmarks in prior work (Guérin et al. Reference Guérin, Cordier, Flageul, Grieu and Agostini2025). For DHT investigations, an additional passive scalar temperature equation is resolved at ${\textit{Pr}}=1$ , with numerical dissipation increased by a factor of 2.25 through the spectral vanishing viscosity scheme to address small instabilities observed in some drag/heat-transfer increase cases. The LES methodology demonstrates accurate prediction of primary performance metrics, with typical discrepancies in friction coefficient and Nusselt number remaining below 3 % and 5 %, respectively, relative to fully resolved DNS. It is noted that the $S_1$ resolution is already close to standard DNS resolution for this Reynolds number; the implicit LES classification arises from the use of spectral vanishing viscosity rather than from substantially coarsened grid spacing. Preliminary tests with coarser grids revealed excessive numerical dissipation during the oscillation-cycle phases in which near-wall structures are fully tilted in the spanwise direction by the Stokes layer, corrupting in particular the thermal-transport metrics. The present resolution thus represents a practical lower bound for reliable performance evaluation within the optimisation framework, as it ensures adequate capture of the spanwise-distorted streaks during maximum-tilt phases without excessive artificial damping.

The $S_2$ configuration provides DNS-level resolution for final validation of optimal configurations identified through the optimisation workflow. This high-fidelity grid comprises $400 \times 221 \times 200$ points distributed over the extended domain $(L_x, L_y, L_z) = (24, 2, 6)$ , leading to substantially refined resolution ( $\Delta x^+ = 10.7$ , $\Delta z^+ = 5.3$ , $\Delta y^+_{\textit{min}} = 0.43$ ) capable of fully resolving all dynamically significant scales. The $S_2$ configuration serves as the authoritative standard against which final performance metrics are established, ensuring that the conclusions drawn from the optimisation-based exploration are substantiated through rigorous high-fidelity validation. Each DNS simulation requires approximately 18 000 CPU hours, rendering systematic parameter space exploration infeasible without the multi-fidelity strategy.

This hierarchical approach enables a computationally pragmatic workflow: extensive parameter exploration is performed through hundreds of LES at reduced computational cost, whilst final reported performance metrics are verified against DNS standards. The computational advantage is substantial, with each LES evaluation requiring approximately 0.5 % of the resources necessary for an equivalent DNS run, thereby permitting comprehensive exploration of the continuous waveform parameter space that would otherwise be prohibitively expensive. The subsequent validation of optimal configurations through DNS ensures scientific rigour whilst maintaining practical feasibility for discovering novel control strategies.

2.5. Computational strategy and statistical convergence

The PBO-LES framework employs a fundamentally different computational strategy compared with traditional parametric studies. Whilst Cimarelli et al. (Reference Cimarelli, Frohnapfel, Hasegawa, De Angelis and Quadrio2013) performed 90 high-fidelity DNS runs, each requiring integration times of approximately 10 000 viscous time units with 95 wash-out times to ensure statistical convergence, the present investigation executes hundreds of shorter LES per optimisation case (typically 400–720 simulations depending on problem complexity and convergence behaviour).

Each PBO-LES comprises two distinct phases. First, a warm-up period of 3200 wall units is enforced to ensure the flow adapts to the imposed actuation, thereby reaching a statistically steady state prior to data collection. This extended warm-up duration, substantially exceeding typical transient time scales, guarantees that the modified flow regime is fully established before statistical sampling commences. Second, statistical collection proceeds over multiple complete oscillation cycles, with the duration determined adaptively according to the period through the relation $t^+_{\textit{collect}} = T^+ \lfloor 1630/T^+ \rfloor$ , where $\lfloor \boldsymbol{\cdot }\rfloor$ denotes the floor function. This formulation ensures that data are collected over an integer number of complete cycles, typically spanning 4–15 cycles depending on the period (corresponding to approximately 1630 wall units of statistical data per simulation). For the drag-reduction case at $T^+ = 111$ , this leads to approximately 14–15 complete cycles; for the DHT waveform case at $T^+ = 350$ , this corresponds to 4–5 complete cycles. This duration, whilst substantially shorter than traditional DNS validation runs, proves sufficient for the optimisation framework owing to a fundamental distinction between single-simulation convergence and population-based convergence. Individual simulations of limited duration inevitably retain residual statistical noise in their performance metrics ( $P_{\textit{sav}}$ , $\overline {A_n}$ ), and 4–5 cycles would be insufficient for fully converged second-order statistics in a traditional validation study. However, the PBO algorithm does not rely on the absolute accuracy of individual evaluations; rather, it exploits the aggregate signal emerging from simultaneous evaluation of multiple individuals per generation ( $n_{\textit{ind}}=10$ ) within the parameter neighbourhood. The algorithm updates its policy based on the collective reward distribution across these individuals, such that statistical noise from any single short-duration simulation is effectively filtered through population averaging. Convergence is evidenced by the consistency of rewards across final generations, wherein both the best individual and the generation average plateau at comparable performance levels. This population-based approach contrasts fundamentally with single-point evaluation strategies, wherein each simulation must achieve high statistical fidelity independently. In the PBO framework, individual simulations serve primarily for relative comparison within and across generations, with absolute performance subsequently validated through high-fidelity DNS. The optimal configurations identified through this efficient exploration methodology undergo rigorous validation using high-fidelity DNS. By integrating computationally efficient LES for exploration with DNS for final validation, this two-stage approach ensures both scientific rigour and practical feasibility. The computational advantage is significant: each LES requires approximately 96 CPU hours, compared with 18 000 CPU hours for an equivalent DNS run. This efficiency enables the execution of hundreds of simulations, facilitating a comprehensive exploration of the continuous parameter space that would otherwise be infeasible with DNS alone.

2.6. Optimisation strategy and problem formulation

Having established the computational infrastructure, i.e. waveform parameterisation, optimisation algorithm, numerical framework, performance metrics and statistical convergence strategy, it is pertinent to delineate the structured experimental programme through which these methodological components are deployed. The investigation is organised around four distinct optimisation problems, which must be addressed in a specific sequential order dictated by methodological necessity. The experimental design follows a systematic progression from validation on partially established benchmarks (wherein optimal-waveform topologies are known from discrete sampling but continuous parameter space exploration remains incomplete) to discovery on unexplored problems, with each case building upon the insights and confidence established by its predecessors. This hierarchical approach ensures that the framework’s capabilities are rigorously established before application to novel research questions wherein ground-truth solutions remain unknown. Table 3 summarises the four cases, differentiated by objective (or reward) function, signal parameterisation and control parameter space.

Table 3.

Summary of optimisation problems. Cases 1–2 validate the waveform framework against established results from drag-reduction studies. Cases 3–4 address the primary research objective: discovering optimal actuation signals for DHT. An auxiliary sinusoidal optimisation at $W^+=30$ was additionally performed to provide an iso-amplitude comparison baseline for case 4 analysis (see § 3.4).

2.6.2. Discovery strategy: cases 3–4

Following successful validation on drag-reduction benchmarks, the framework is applied to the primary research objective: discovering optimal actuation signals for dissimilar-heat-transfer enhancement. Unlike drag reduction, wherein optimal parameters are well established through extensive prior investigations, the thermal flow-control problem lacks such foundational parametric knowledge. The discovery phase therefore requires a systematic two-stage approach to establish rigorous performance baselines before exploring waveform effects.

The investigation proceeds through two sequential cases. Case 3 first identifies the optimal sinusoidal parameters through comprehensive exploration of the ( $T^+$ , $W^+$ ) space, establishing the best achievable performance using conventional sinusoidal forcing. This ‘sinusoidal champion’ then serves as the performance benchmark for case 4, wherein waveform shape is optimised whilst holding amplitude below the sinusoidal optimum. This sequential structure is essential for isolating waveform topology contributions: without first identifying optimal sinusoidal parameters, any performance improvements in case 4 could not be unambiguously attributed to beneficial temporal structures, as they might instead arise merely from improved parameter selection relative to an arbitrary baseline. The two-stage progression thus ensures that any superiority of non-sinusoidal waveforms is demonstrably attributable to shape effects rather than confounded by comparison against a suboptimal sinusoidal reference.

Case 3 optimises period and amplitude ( $T^+$ , $W^+$ ) to maximise the analogy factor $\overline {A_n}$ over an extensive parameter space ( $10 \leqslant T^+ \leqslant 750$ , $10 \leqslant W^+ \leqslant 40$ ). This broad exploration ensures comprehensive coverage and reliable identification of the global sinusoidal optimum, avoiding premature convergence to local maxima that might arise from restricted search domains. The resulting configuration, designated the ‘sinusoidal champion’, provides substantial improvement over the non-optimised reference of $\overline {A_n}=1.064$ reported by Guérin et al. (Reference Guérin, Flageul, Cordier, Grieu and Agostini2024), thereby establishing a rigorous performance target for waveform-based strategies. Importantly, case 3 identifies the optimal amplitude $W^+_{\textit{opt}}$ , subsequently employed as a reference for parameter selection in case 4. The comprehensive nature of this exploration ensures that the sinusoidal baseline represents the true optimum rather than a convenient but potentially suboptimal comparison point.

Case 4 investigates the central research hypothesis: whether non-sinusoidal temporal structures can achieve superior DHT performance relative to the optimal sinusoidal baseline, ideally at reduced energetic cost. The amplitude is fixed at $W^+=30$ , selected to operate within the moderate-amplitude regime whilst remaining below the upper bound explored in case 3. This strategic amplitude selection serves dual purposes: it enables assessment of whether waveform optimisation can compensate for reduced forcing intensity, potentially providing a more energy-efficient pathway to DHT enhancement, whilst simultaneously reducing optimisation dimensionality to concentrate computational effort on waveform parameter space exploration. Success is defined as identification of a waveform topology leading to $\overline {A_n}$ (see (2.11)) approaching or exceeding the sinusoidal champion performance despite operating at moderate amplitude. Such an outcome would demonstrate that optimisation of the temporal waveform offers a viable alternative to simple amplitude escalation for enhancing DHT, with potential implications for energy-efficient thermal management applications.

The complete parameter ranges and reward functions employed for each optimisation case, with ranges selected to enable comprehensive exploration whilst maintaining physical realisability are conveyed by table 4. Waveform parameters are constrained to $(a_1, a_2) \in [0,1]^2$ for all waveform cases, ensuring the generated signals remain physically realisable and bounded within the feasible domain defined by the Lagrangian parameterisation. Period and amplitude ranges are adjusted according to the specific physical regime under investigation and the particular objective function being optimised, with bounds informed by prior parametric studies and physical considerations regarding Stokes-layer penetration and control authority.

Table 4.

Systematic optimisation programme: parameter ranges, signal types and reward functions for the four sequential cases. Cases 1–2 (blue, drag reduction) validate the framework’s capability to identify both complex waveforms (case 1) and simple solutions (case 2) without bias. Cases 3–4 (red, DHT) constitute the discovery phase, establishing optimal sinusoidal baseline (case 3) before exploring waveform effects (case 4). Parameter ranges $(a_1, a_2) \in [0,1]^2$ apply to waveform cases; sinusoidal optimisation fixes $a_1= 1/({\sqrt {2}(1+\sqrt {2})}), a_2= {2}/({1+\sqrt {2}})$ to generate pure sinusoidal signals.