Jet mixing enhancement with Bayesian optimization, deep learning, and persistent data topology

The configuration is a jet flow exiting a circular nozzle of diameter $D$ in figure 1. The flow is described in a Cartesian coordinate system $(x,y,z)$ where $x$ represents the streamwise direction and the origin coincides with the center of the nozzle. The computational domain starts from the exit and covers a rectangular region with size $12D\times 16D\times 12D$. The actuation $\bm{u}_{b}(r,\theta,t)$ is imposed with the mean streamwise velocity $u_{m}(r)$ as the inlet velocity profile $u(r,t)$,

where $r$ measures the radial distance from the centerline, and $\theta$ is the azimuthal angle. The mean streamwise component has a hyperbolic-tangent profile,

where $U_{\rm j}$ is the jet centerline velocity, $u_{c}=0.03U_{\rm j}$ is the co-flow velocity to mimic a natural suction process, and $\delta_{2}=R/20$ is the momentum boundary layer thickness of the initial shear layer. At the side boundaries, we impose the vertical velocity equals $u_{c}$, and the remaining velocity components equal zero. The pressure at the side boundaries is computed from the Neumann condition ${\mathbf{n}}\cdot\nabla p=0$ with ${\mathbf{n}}$ as the vector normal to the boundary. At the outlet plane, the velocity is computed from a convective boundary condition $\partial\bm{u}/\partial t+\tilde{V}_{C}\partial\bm{u}/\partial n=0$, where $\tilde{V}_{C}$ is the instantaneous convection velocity $V_{C}$ limited to positive values: $\tilde{V}_{C}=\max(V_{C},0)$. $V_{C}$ is the velocity averaged over the outlet plane. The pressure at the outflow equals zero. Such defined outflow boundary condition ensures stable simulations and has negligible impact on the turbulent flow structures leaving the computational domain (Tyliszczak & Geurts, 2014; Tyliszczak, 2018).

As introduced in § 1, the jet control techniques for mixing enhancement are usually designed according to the instability mode, described by their azimuthal wavenumber at order $0$ (axisymmetric mode) or $1$ (helical mode). The perturbation is either axial or radial. We combine both axial and radial perturbation and define the actuation $\bm{u}_{b}$ with a general expression of $\theta$ and $t$, without assumption on the forcing mode. Therefore, the term $\bm{u}_{b}$ includes an axisymmetric streamwise bulk forcing $u^{\alpha}(r,t)$, and a peripheral forcing with the streamwise component $u^{\beta}(r,\theta,t)$, and the radial $u^{\gamma}(r,\theta,t)$:

The forcing components are the product of a perturbation $f^{m}(\theta,t)$ and a radial profile $g^{m}(r)$: $u^{m}(r,\theta,t)=f^{m}(\theta,t)g^{m}(r)$, $m=\alpha,\beta,\gamma$ with $g^{\alpha}(r)=1$ for $r\leq R$ and $0$ for $r>R$, and $g^{\beta}(r)=g^{\gamma}(r)=\exp(-1000(R-r)^{2.5})$. The profiles of the three forcing components are depicted in figure 1. The perturbation terms $f^{m}(\theta,t)$ are defined as the sum and product of space- and time-harmonic functions:

where $\alpha_{i}$, $\beta_{ij}$, $\gamma_{ij}$ and $\omega_{\alpha}$, $\omega_{\beta}$, $\omega_{\gamma}$ are the actuation amplitudes and angular frequencies, respectively. $\Theta_{i}(\phi)$ is the harmonic function basis defined as: $\Theta_{i}(\phi)=\sin(i\phi)$ for $i>0$, $\Theta_{i}(\phi)=1$ for $i=0$, and $\Theta_{i}(\phi)=\cos(i\phi)$ for $i<0$. The forcing ansatz can approximate any periodic function of $\theta$ and $t$ as the expansion order increases. In this study, focus is placed on the first order expansion ($L=1$) of (2-4). Thus, the control law is parameterized by a $22$-dimensional vector $\mathbf{b}$,

where the Strouhal numbers $St_{m}=\omega_{m}D/2\pi U_{\rm j},m=\alpha,\beta,\gamma$. Note that $\alpha_{0}$ is set to 0 as a constant bulk flow can be incorporated into the steady profile. In addition, $\alpha_{-1}=0$ can be assumed by a translation in time. The range of $St_{m}$ is set as $[0.1,1]$ to include the Strouhal number of the preferred mode at $St_{p}=0.3-0.64$ (Crow & Champagne, 1971; Gutmark & Ho, 1983; Sadeghi & Pollard, 2012), and the range of axisymmetric mode $St_{\alpha}\in[0.15,0.8]$ where bifurcating and blooming jets are observed (Lee & Reynolds, 1985; Parekh, 1989; Tyliszczak, 2018). The actuation amplitudes are limited to $-0.1\leq\alpha_{1},\beta_{ij},\gamma_{ij}\leq 0.1$, lower than $0.15$ used by Danaila & Boersma (2000); Gohil et al. (2015); Tyliszczak (2018), and $0.5$ by Koumoutsakos et al. (2001).

This high-dimensional search space allows the actuation to emulate various forcing modes, such as axisymmetric, helical, flapping, bifurcating, dual-mode, and harmonic waves discussed in §  1. In short, the forcing can be axisymmetric, statistically axisymmetric, and non-axisymmetric.

The mixing process in turbulent jets is typically characterized by quantities such as the decay of centerline velocity and its fluctuations, or the entrainment (Nathan et al., 2006). However, these quantities only measure the local statistics and cannot reflect the mixing process of non-axisymmetric jets, such as bifurcating jets and asymmetric jets. For example, the velocity in these jets may locally drop to zero due to a jet-splitting phenomenon, rather than as a result of enhanced mixing. The entrainment is rather a measure of the amount of surrounding fluid entrained into the jet vicinity, without guaranteeing that it mixes with the jet. In asymmetric jets, the amount of fluid flowing towards the jet may characterize significant radial non-uniformity not revealed by the entrainment. Considering the non-axisymmetric forcing in search spaces, we define a new metric, the equivalent mixing radius $R_{\rm eq}$, to estimate the spatial uniformity of the streamwise velocity. $R_{\rm eq}$ is defined as the normalized streamwise velocity variance computed at a given $y-z$ cross-section:

with $\varrho(y,z)=(\langle{u(x=X_{0},y,z)}\rangle_{t}-u_{c})/(U_{\rm j}-u_{c})$, $X_{0}=8D$, and $(y_{c},z_{c})$ the jet center, as an analog to the center of mass: $y_{c}=\int\!\!\!\int y\varrho(y,z)dydz/\int\!\!\!\int\varrho(y,z)dydz$ and $z_{c}=\int\!\!\!\int z\varrho(y,z)dydz/\int\!\!\!\int\varrho(y,z)dydz$. The $\sqrt{2}$ coefficient is chosen so that the equivalent mixing radius of a top flat jet flow of radius $R$ is $R$.

The amplitude and mass flow rate have been adopted to evaluate the control input. Inspired by Parekh (1989), we define the momentum flux $P$ of the actuation to estimate the energy input from a practical perspective. The momentum flux is time-averaged and normalized by the jet axial momentum flux at the inlet:

The optimization problem to maximize the mixing as a response to the actuation input parameterized by $\mathbf{b}$ is formulated as

where $\mathcal{B}=[0.1,1]^{3}\times[-0.1,0.1]^{19}\subset\mathbb{R}^{22}$. Cost function $J$ is defined as the inverse of the equivalent mixing radius and normalized by the unforced case, $J=R_{\rm eq,0}/R_{\rm eq}$. Better mixing with large $R_{\rm eq}$ is related to the decrease of $J$. The assumed optimization goal leads to the spatial uniformity of $\langle{u(x=X_{0},y,z)}\rangle_{t}$.

To optimize this $22$-dimensional search space, we employ techniques inspired by Bayesian optimization (BO) (Williams & Rasmussen, 2006). BO has shown to be advantageous in optimizing expensive black-box functions by systematically reducing uncertainty in the black-box mapping and incorporating prior assumptions of the cost function (Shahriari et al., 2015). Through a sequential approach, BO identifies the next actuation to evaluate, or “data point” to acquire, for the purpose of finding the optimum. This is generally achieved via a surrogate model trained on all the queried data and an acquisition function (Williams & Rasmussen, 2006). A sketch of the method used in this study is shown in figure 1. The algorithm is initialized with the evaluation of a set $\mathcal{D}_{0}$ of $N_{0}$ actuation vectors in $\mathcal{B}$ generated by Latin hypercube sampling. $N_{0}$ is equal to $N_{D}+1$ with $N_{D}$ the dimension of the search space $\mathcal{B}$. We recall that for this study, $N_{D}=22$. $\mathcal{D}_{0}$ includes all the evaluated parameter vectors and their cost $\{\mathbf{b}_{i},J_{i}\}_{i=1}^{N_{0}}$. A surrogate model $\tilde{J}$ is trained on the available data to approximate the latent objective function $J$. After the initialization, the algorithm explores the search space $\mathcal{B}$ one new query at a time. At each iteration, BO determines the optimal actuation to implement next by minimizing an acquisition function $a(\mathbf{b})$. The acquisition function leverages the surrogate model $\tilde{J}$ and available data $\mathcal{D}_{n-1}$ to guide the data selection in the search space. After each query, the data set is enriched by the new data point $\{\mathbf{b}_{n},J_{n}\}$ into $\mathcal{D}_{n}$ to further refine the surrogate model. When the query budget is met, the algorithm ends with the best design vector $\mathbf{b}^{*}$ recorded during the optimization.

The two key elements in BO are the choice of the surrogate model and the sequential strategy (Blanchard & Sapsis, 2021). We focus on the former for a better estimation of the high-dimensional flow control system. Gaussian processes (GP) serve as a successful surrogate model in moderate dimensions and can provide closed-form solutions with the posterior distribution. However, the computation of the posterior costs $\mathcal{O}(n^{3})$ where $n$ is the number of observations due to the inverse of the covariance matrix. The number of evaluations required to effectively cover the domain grows exponentially with the dimensionality. This makes GP difficult to scale to large training sets for high-dimensional problems. Recently, the deep operator network (DeepONet) has shown small generalization error for systems where functions are acted upon by an operator (Lu et al., 2021). Different from GP that parameterizes the input, DeepONet can map input functions, which are then transformed by an operator, to an output function or scalar with improved accuracy. This means that DeepONet does not fall victim to the scaling difficulties of GP when training. Therefore, DeepONet is capable of learning from infinite-dimensional functions. Empirically, DeepONet’s utility as an operator surrogate model for Bayesian-inspired experimental design has been shown to significantly outperform GP in several infinite-dimensional systems that exhibit extreme events in Pickering et al. (2022), ranging from stochastic pandemic spikes, to catastrophic structural failure, to rogue wave identification. Through a study of the Bayesian optimizer based on GP (BO) and DeepONet (BO-DeepONet) for the defined 22-dimensional problem (see §4.2), we propose a new algorithm, deep learning-enhanced Bayesian optimization (BO-DL). By incorporating both GP and DeepONet as the surrogate model, BO-DL presents a better explorative capability and faster convergence. We alternate between DeepONet and GP every $10$ iterations to query the next sample. The $10$ is chosen empirically to balance the characteristics of the two models and combine their advantages. If the interval is too long, such as $100$, the models will be more independent than interacting with each other. On the other hand, if the interval is too short, such as $1$, the exploitation may be interrupted by uncertainty due to the exchange of models. In GP implementations, the parameter space of a stochastic process is used for both regression and searching. Instead, DeepONet performs regression in the functional space, leveraging the typically disregarded basis functions associated with the parameterization. Here, the function $\hat{\mathbf{u}}_{b}$ is designed as the actuation command at $9$ points located at the jet exit, including the centerline and $8$ equidistant points on the periphery.

The acquisition function employed is the likelihood-weighted lower confidence bound proposed by Blanchard & Sapsis (2021), with superiority in finding rare extreme behaviors,

Here, $\kappa$ balances exploration (large $\kappa$) and exploitation (small $\kappa$), and is chosen as $1$. The mean model $\mu(\mathbf{b})$ and standard variance model $\sigma(\mathbf{b})$ are estimated by the mean and variance over a $2$-ensemble of trained DeepONet. For the case of GP, these can be calculated in closed form using standard expressions from GP regression (Williams & Rasmussen, 2006). The likelihood ratio $w(\mathbf{b})$ measures relevance by weighting the uncertainty of the point (the input density $p_{\mathbf{b}}$) against its expected impact on the cost function (the output density $p_{\mu}$).

We consider an incompressible flow described by the Navier-Stokes equations in the framework of LES:

where ${u}_{i}$ represents the velocity components, $p$ denotes pressure and $\rho$ is density. The overbar denotes spatially filtered variables, $\bar{f}(\mathbf{x},t)=\int_{\Omega}f(\mathbf{x}^{\prime},t)\mathcal{G}(\mathbf{x}-\mathbf{x}^{\prime},\Delta)d\mathbf{x}^{\prime}$ and $\mathcal{G}$ is the filter function that fulfils the condition $\int_{\Omega}\mathcal{G}(\mathbf{x},\Delta)d\mathbf{x}=1$. A local filter width equals the cube root of the computational cell volume, $\Delta=(\Delta x\Delta y\Delta z)^{1/3}$. The stress tensor includes the large-scale term $\tau_{ij}$ and the sub-grid term $\tau_{ij}^{f}$ defined as

where $\nu$ is the kinematic viscosity and $S_{ij}=\frac{1}{2}\left(\frac{\partial\bar{u}_{i}}{\partial x_{j}}+\frac{\partial\bar{u}_{j}}{\partial x_{i}}\right)$ is the rate of strain tensor of the resolved velocity field. In this work, the sub-filter tensor is modeled by an eddy-viscosity model with $\tau_{ij}^{f}=2\nu_{t}S_{ij}+\tau^{f}_{kk}\delta_{ij}/3$. The diagonal terms $\tau^{f}_{kk}$ are added to the pressure, resulting in the so-called modified pressure $\bar{P}=\bar{p}-\rho\tau^{f}_{kk}\delta_{ij}/3$. The Vreman subgrid-scale model is used for the low computational cost and very good accuracy in simulating jet flows (Wawrzak et al., 2015; Boguslawski et al., 2019).

The simulations are conducted with an in-house high-order LES solver SAILOR. The solution algorithm is based on the projection method for the pressure-velocity coupling for half-staggered meshes where the pressure nodes are shifted half a cell size from the velocity nodes (Tyliszczak, 2014, 2015a). A predictor-corrector method (Adams-Bashforth/Adams-Moulton) is applied for the time integration. Derivative approximations and interpolation on staggered nodes are defined using 6^th and 10^th order compact difference formulas. The SAILOR solver was used in the jet studies with similar dynamic scales as the present work, such as jets undergoing laminar/turbulent transition (Boguslawski et al., 2019) and excited jets (Tyliszczak & Geurts, 2014; Tyliszczak, 2018). The applied high-order discretization schemes led to grid-independent results already at relatively coarse meshes.

The Reynolds number $\Rey=U_{\rm j}D/\nu$ is decided as $3000$ for this study. This allows for the use of a relatively coarse computational mesh to obtain reliable and fast simulations as the database. Two meshes are employed in this study. A coarse mesh with $80\times 160\times 80$ nodes is used for the learning process and a refined mesh with $192\times 336\times 192$ nodes for the validation and flow analysis of selected cases. The mesh points are compacted in the axial direction towards the inlet using an exponential function and radially towards the jet axis by a tangent hyperbolic function. In the region $-1.2D<y,z<1.2D$, the mesh spacing is almost uniform and equal to $\Delta y=\Delta z=0.05D$ (46 nodes) on the coarse mesh and $\Delta y=\Delta z=0.02D$ (115 nodes) on the dense one. In the axial direction, the sizes of the cells in the direct inlet vicinity are $\Delta x=0.067D$ and $\Delta x=0.032D$ for the coarse and dense mesh, respectively. The time step varies according to the Courant-Friedrichs-Lewy condition, with the CFL number equal to 0.5. The jet impulsively injects into quiescent flow and becomes fully developed after $100D/U_{\rm j}$ time units. The time-averaging procedure then starts and lasts for $500D/U_{\rm j}$ time units for the statistics to converge. A single simulation on the coarse mesh takes $20$ CPU-hours. The whole optimization process with $1000$ converged simulations lasts around $21$ days, using 40 CPUs of an AMD EPYC 7742 (2.25GHz) processor. On the dense mesh, the single simulation run takes approximately $576$ CPU-hours. The parallel computation is carried out with the Open MPI interface.

As presented in 3.2, the numerical code employed has been well validated against experimental and numerical data for a series of studies of jet dynamics and control. Here, the code with the assumed perturbation design (2-4) is verified to obtain a well-documented flow pattern of an excited jet, the bifurcating jet at $\Rey=4300$ (Lee & Reynolds, 1985). The well-known jet is shown in figure 2, produced by the combined axisymmetric and flapping excitation at the frequencies $0.4<St_{\alpha}<0.6$ and $St_{\beta}=St_{\gamma}=St_{\alpha}/2$. Based on the current control definition, the bifurcating jet is reproduced by

that produces the axisymmetric excitation and

as the flapping mode, simulating the orbital motion of the nozzle tip in the experiment. We take $\alpha_{1}=0.17$, $St_{\alpha}=0.5$ and assume $St_{\beta}=St_{\gamma}=St_{\alpha}/2$, $\beta_{-1,1}^{2}+\gamma_{-1,1}^{2}=\alpha_{1}^{2}$ with $\beta_{-1,1}=\alpha_{1}\cos(20^{\circ})$. This type of excitation is also used in the previous LES simulations of the bifurcating jet (da Silva & Métais, 2002; Tyliszczak & Geurts, 2014).

Figure 2(b1-b3) shows the time-averaged axial velocity profiles along the radius in the bifurcating plane at the distance $x/D=5.0,\,6.5,\,8.0$ from the inlet. These results were obtained for $\Rey=4300$, as in Lee & Reynolds (1985), and for $\Rey=3000$ assumed in the present study. The effect of the Reynolds number on the velocity profiles is small. We attribute such behavior to a dominating role of the perturbation.The location and level of two peaks, which are associated with the split jet arms, are well predicted by the numerical solutions. The impact of the mesh density on the solution is also negligible, owing to the employed high-order numerical method. This also holds for the optimized jet, see figure 6 of § 4.4.

We first study the capability of the surrogate model, GP, and DeepONet, to predict the cost function $J$ as a response to the excitation input $\mathbf{b}$. Then, the Bayesian optimizers based on each of the two surrogate models (BO and BO-DeepONet) are tested on our plant. Finally, the performance of the proposed method deep learning-enhanced Bayesian optimization (BO-DL) in 3.1 is illustrated.

A $k$-fold cross-validation training ($k=5$) of GP and DeepONet model is performed over $1000$ data points with $80/20$ train/test split. The data are extracted randomly out of the database from realizations of Bayesian optimizers. Figure 3(a) shows the prediction $\tilde{J}(\mathbf{b})$ versus the truth $J(\mathbf{b})$ obtained. The distribution of points along the diagonal shows that DeepONet achieves a lower prediction error than GP. This is further explained by the correlation coefficient $R=0.89$ for DeepONet, and $0.71$ for GP. The average error of the $k$ tests is measured by the mean squared error (MSE). The MSE for GP model is $0.01$, $1\%$ of the range of $J$ value.The prediction of DeepONet model is superior, with a MSE equal to $0.005$. The learning process of Bayesian optimizer with GP (BO) and DeepONet (BO-DeepONet) is shown in figure 3(b1-2). In figure 3(b1), the learning curve of BO displays a plateau after the initial samples (triangles). After $160$ samples, new optima are found and followed by continuous exploitation of the samples near the learning curve. The final solution is reached with $J=0.274$ after $745$ evaluations. When DeepONet is employed (figure 3b2), a better solution $J=0.256$ is found quickly within $300$ samples. This may be attributed to DeepONet’s capibility to generalize better for previously unseen data than GP as the cross-validation indicates (Lu et al., 2021). After $m=300$, the newly tested parameters cover the entire range of $J$, but no further improvement is observed in the learning curve. This suggests that the optimizer focuses on exploration of the search space rather than exploitation like BO. Based on the above observations, a joint surrogate model is proposed for this study to combine the advantages of GP in local exploitation and DeepONet in exploring new minima. The Bayesian optimizer based on this new model is described in § 3.1 and referred to as BO-DL. The learning process of BO-DL is given in figure 3(b3) with the samples queried by GP and DeepONet. As indicated by the data points on the learning curve, the queries made by DeepONet (red dots) discovers a new minimum with significant reduction of $J$, and GP (blue dots) continues to descend. The best solution is obtained at $J=0.237$ within $600$ evaluations.

The average performance of three Bayesian optimizers above is further studied. Each optimizer is employed for three realizations with a fixed budget of $1000$ evaluations. Figure 3(c) reports the average value of the current optimum $J_{min}$ from each optimizer with the standard deviation (shaded region) of three runs. The learning curve starts from $J_{min}=0.45$, the lowest cost value after initialization of $23$ samples, including the unforced case and the other $22$ controlled cases from Latin hypercube sampling in the search spaces. The unforced case ($J=1$) is omitted for better visibility of the data. The maximum cost of the controlled flow is around $0.7$. With around $750$, $580$ and $570$ queries, the average lowest costs $J_{min}$ achieved by BO, BO-DeepONet, and BO-DL are $J=0.274$ (diamond), $J=0.268$ (square), and $J=0.237$ (star), respectively. On average, BO-DeepONet shows the fastest learning speed (dotted line) but with the largest variation. This is owing to DeepONet’s capability of predicting the potential minima with a small generalization error. Although the descent of BO is the slowest, the optimal results of the three realizations are consistent. This indicates that GP provides better interpolation around the minima than DeepONet due to its deterministic nature. By combining GP and DeepONet, BO-DL not only demonstrates a comparable learning speed as BO-DeepONet but also inherits the small variance of the final solution from BO. Finally, among the three optimizers, BO-DL derives the best solution. In addition, the warm-up phase during queries $0$ to $100$ appears to be significantly shortened, denoted by the rectangle in figure 3(c).

The computational cost of the BO loop is also noteworthy. With BO, the computation of the posterior costs $\mathcal{O}(N^{3})$, where $N$ is the number of observations (Williams & Rasmussen, 2006). This makes the algorithm quite slow, even after only a few hundred observations. The experience of this study shows that the computation of BO increases from $10$ CPU-seconds to $600$ CPU-seconds after $1000$ iterations on an AMD EPYC 7742 (2.25GHz) processor. The BO-DeepONet procedure scales much more favorably. Initially, the first iteration takes $120$ CPU-seconds, increasing only to $180$ CPU-seconds after $1000$ iterations. The combination of the two models in BO-DL compromises the cost to an average level.

For the high-dimensional physical problem, we show that Bayesian optimization can benefit not only from more accurate surrogate models but also from combining the advantages of parametric and non-parametric predictors. The proposed BO-DL holds a fast convergence and efficient exploration with a GP-DeepONet surrogate model. Compared with GP, the proposed surrogate model can provide more accurate predictions by leveraging the hidden functional input with DeepONet and scales better as both data size and dimensionality increase. In addition, a comparison between the Bayesian methods and bio-inspired approaches is given in Appendix A. It is shown that the optimizer with a surrogate model employed, particularly a deep-learning model, shows more advantage in the current problem.

In this section, we explore the learning processes of the BO and BO-DL in the $22$-dimensional space with persistent data topology (Wang et al., 2023a, b). This data analysis identifies the cost function minima and their depth, i.e. their persistence to noise, and was inspired by Edelsbrunner & Harer (2008). The analysis includes the identification of local minima in the high-dimensional actuation space, a dimension reduction to a two-dimensional proximity map and corresponding data visualization, and as shown in figure 4. The $22$-dimensional data obtained by both BO and BO-DL are projected on a 2D proximity map by classical multidimensional scaling. The feature coordinates $\gamma_{ij}$ are chosen to optimally preserve the dissimilarity between control parameters defined by the Euclidean distance $D_{ij}=|\mathbf{b}_{i}-\mathbf{b}_{j}|$. The map features two large basins of attraction with low values of $J$, as well as small basins distributed around the border. A point $\mathbf{b}_{0}$ is supposed as a local minimum $\mathbf{m}$, if there exists a neighborhood $\mathcal{B}$ of $\mathbf{b}_{0}$ that satisfies $J(\mathbf{b}_{0})\leq\min_{\mathbf{b}\in\mathcal{B}}J(\mathbf{b})$. $\mathcal{B}$ is an open set which should include $K$ nearest neighbours of $\mathbf{b}_{0}$ measured by Euclidean distance, $K\geq N_{D}+1$. Note that the local minima are assumed based on the obtained discrete data and may change with additional data. A total of $57$ local minima are extracted from the data, with $36$ found by BO-DL and $21$ by BO. In the proximity map, the unforced case is represented by a black square where both algorithms begin. The other symbols denote the derived minima $\mathbf{m}$ found by BO (blue) and BO-DL (red). The final BO and BO-DL solutions highlighted by the filled diamond ($J=0.27$) and the filled star ($J=0.24$) are located in the large basins of attractions. Most of the minima queried by BO are located in the center of the map, whereas BO-DL also explores outward regions. Forced by the control commands corresponding to these minima, different jet patterns are observed, corroborated with the control modes. The axial puff (circles) are close to the unforced case. The bifurcating type (cross) distributes widely in the cost range. The lower the $J$ value is, the closer to the helix (filled diamond) basin. The jets bifurcating to one side are away from the center, surrounded by the other unidentified patterns (triangles). Helix (diamond) and blooming (star) jets feature the most substantial performance, but the latter is only detected by BO-DL. Among the $20$ minima explored by BO, the identified patterns include $6$ helix, $5$ flapping, and $4$ axial puff. Among the $37$ minima explored by BO-DL, the identified patterns include $6$ flapping, $5$ asymmetric flapping, $2$ helix, $1$ blooming, and $1$ axial puff. In addition, most ($22$) of the $27$ unidentified patterns are detected by BO-DL.

BO-DL explores not only more minima than BO but also more diverse flow patterns beneficial to the mixing. This is probably owing to DeepONet’s capability to extrapolate the mapping from the high-dimensional actuation to the mixing response more accurately. Two solutions with large basins of attractions in the search space are revealed — the optimal solution with a $7$-armed blooming jet generated, and the suboptimal with a double-helix shape.

Here, we include three solutions for the discussion: an ad hoc forcing with the best mixing in Tyliszczak (2018), BO optimized solution, and the optimal solution of BO-DL. The forcing command, the instantaneous snapshots, and the mean flow fields are presented in figure 5. The forcing commands are expressed by the operators in an order of constant, spatial-periodic, temporal-periodic, and traveling waves in figure 5(a1), 5(b1) and 5(c1). The axial excitation combining axisymmetric and helical modes has been widely employed to study the bifurcating and blooming jets since Lee & Reynolds (1985). A parametric study of the blooming jets with this type of excitation was performed in Tyliszczak (2018), under the same Reynolds number as this study. Among various multi-armed jets, the one with 11 arms led to the best mixing performance. The excitation was imposed on the axial velocity and combined the axisymmetric mode with Strouhal number $St_{a}=0.45$ and the helical mode with $St_{h}=0.164$ at the same amplitude, $15\%$ the bulk jet velocity (figure 5a1). The BO solution contains mainly the axisymmetric mode at an amplitude of $8\%$ with $St_{\alpha}=0.497$ for the bulk, a helical mode at an amplitude of $7\%$, and a flapping mode at an amplitude of $4\%$ with $St_{\gamma}=0.232$ for radial components in the periphery. After removing the expansions with negligible amplitudes, less than $1\%$ of the bulk jet, the control law approximately reads

Because the removed terms hold an amplitude lower than the turbulence intensity at the jet outlet, the approximation hardly changes the flow patterns, with the relative cost difference being less than $1\%$. The BO-DL solution contains mainly the axisymmetric mode at an amplitude of $6\%$ with $St_{\alpha}=0.519$ for the bulk, a helical mode at an amplitude of $8\%$ with $St_{\gamma}=0.223$ for radial components in the periphery. The simplified control law reads

Two significant factors to be noted are the axisymmetric forcing Strouhal number $St_{\alpha}$, and the frequency ratio between the axial and helical modes, $\alpha=St_{\alpha}/St_{\gamma}$. For both BO and BO-DL solutions, the axisymmetric forcing Strouhal number falls into the range $0.4\lesssim St_{\alpha}\lesssim 0.6$ to observe bifurcating and blooming jets, and coincides with around $St_{\alpha}=0.5$ where the peak spreading occurs (Lee & Reynolds, 1985; Shaabani-Ardali et al., 2020; Gohil et al., 2015). The BO-DL actuation takes a frequency ratio of $2.34$, which very well agrees with a theoretically derived value $\alpha=7/3$ Tyliszczak (2015b); Gohil et al. (2015). Interestingly, the ratio of the BO solution which produces a helix jet ($\alpha=2.14$) also falls into this range. Moreover, different from the ad hoc excitation using only the axial forcing, the radial component in the periphery plays an important role in solutions optimized by both BO and BO-DL. Shaabani-Ardali et al. (2020) also concludes radial forcing is the dominant component of helical modes to maximize the spreading angle of a bifurcating jet. We extend the importance of radial forcing to the jet spreading globally. From an estimate of the momentum flux, the solutions in this study take only $2.8\%$ (BO-DL) and $4.8\%$ (BO) of the main jet, one order lower than the ad hoc excitation ($25\%$). One reason is the low amplitudes, and another is the forcing applied into the local boundary region (see § 2.1) rather than the whole jet, which leads to a more efficient control. This represents the physical reality of small actuators installed on the wall of the inlet nozzle, like flap arrays in Suzuki et al. (1999), only affecting the boundary layers.

The flow structures are presented by the bottom view of the instantaneous Q-parameter isosurfaces (figure 5a2, b2, and c2). The arms of the ad hoc blooming jet are not explicitly observed due to the interaction between the closely aligned vortex rings. A double-helix jet is formulated by the BO solution. The jet bifurcates into two branches, which rotate with a specific frequency (see figure 8) and then experiences continuous bifurcation along the azimuth until the vortex rings break. This type of jet has not been reported in the literature so far. We reserve it for future investigation. The BO-DL optimized jet produces a $7$-armed blooming jet, with the vortex rings eventually propagating along $7$ different trajectories. The contour slices of the time-averaged streamwise velocity also confirm the spreading observed from the vortex rings. The 11 branches generated by the ad hoc forcing can be traced in $x=8D$. The BO-optimized jet shows more continuous distribution along the circumference due to the azimuthal bifurcation of two helix-shaped arms. The blooming jet features the earliest and furthest spreading. This leads to the largest effective mixing radius, $4.22R_{\rm eq,0}$ at $x=8D$, followed by BO optimized jet with $R_{\rm eq}=3.65R_{\rm eq,0}$, and the ad hoc forced jet with $R_{\rm eq}=2.77R_{\rm eq,0}$.

Figure 6 shows the axial profiles of the time-averaged centerline velocity and its fluctuation for the unforced and forced jets. The results obtained on the coarse and dense meshes agree well, except for slight discrepancies in the region $5D<x<8D$. Compared with the unforced case in figure 6(a), the length of the potential core shortens significantly from $7.5D$ to $2D$ for both helical and blooming jets. Beyond the potential core, the velocity drops steeply and even reaches small negative values in the blooming jet. As a similar behavior observed by Tyliszczak (2018), this is the effect of jet splitting resulting in a local pressure drop. As a result, the reversal flow amplifies the local turbulence intensity to the peak at $x\approx 5D$, as indicated in figure 6(b). The initial fluctuation level in both jets corresponds to the imposed forcing amplitudes of the bulk forcing term. A small decrease at $x<1D$ is caused by the lack of energy in a range of low-wave numbers (Kempf et al., 2005). For the controlled jets, the waves of the fluctuation profiles around the peak are attributed to the interactions between the Kelvin-Helmholtz instability and the forcing disturbance. Further downstream, the fluctuations drop nearly to zero along with a low mean velocity. The fluctuation profile for the unforced jet shows a very low turbulence level until $x\approx 7D$, and then slowly increases to the maximum around $x=10D$.

Figure 7 shows the amplitude spectra of the centerline velocity at four localizations along the axis, $x=0D$, $2D$, $4D$, and $6D$. These results are presented versus the Strouhal number $St_{D}=\omega D/2\pi U_{\rm j}$. The spectrum of the unforced jet is nearly flat at the inlet as the imposed turbulent signal does not contain any characteristic frequency. The high-frequency components ($St_{D}>1$) are dampened downstream, and a broadband peak emerges around $St_{D}=0.52$. This falls within the range of the preferred mode frequency $St_{p}=0.3-0.64$ (Crow & Champagne, 1971; Gutmark & Ho, 1983; Sadeghi & Pollard, 2012). Note that the optimal $St_{\alpha}$ predicted by BO-DL for the blooming jet perfectly matches the current preferred mode $St_{p}=0.52$. This finding is consistent with previous studies (Tyliszczak & Geurts, 2014; Gohil et al., 2015; Tyliszczak, 2018) which concluded that the jet splitting phenomenon is most pronounced when $St_{\alpha}$ is equal to $St_{p}$. The initial spectra of the helix and blooming jets characterize a distinct peak at $St_{\alpha}$. The peaks related to the helical forcing $St_{\gamma}$ can be observed from $x=2D$. The high-frequency harmonics also appear due to the interactions between generated toroidal vortices. In the case of the helical jet, the peak at $St_{D}\approx 0.032$ is also noteworthy. We find that this frequency coincides with the azimuthal motion of the helical arms, with a period $T_{r}$ equal to $62.5D/U_{\rm j}$. Figure 8 shows the snapshots of the helical jet, depicting the positions of its arms during the period of $31.25D/U_{\rm j}$, which corresponds to the detected $St_{D}\approx 0.032$. The relationship between the frequency of the rotation and the one associated with forcing terms is left for future study.