2.1. Governing equations

Following Kumar & Pothérat (Reference Kumar and Pothérat2020), we consider a cavity of height $H$ with an upper free surface where hot fluid is fed in, and a cold, porous lower boundary representing a solidification front, as shown in figure 1. The cavity is also assumed infinitely extended in the third direction ( $\mathbf e_z$ ). The lower boundary is semicircular and connected to the flat upper boundary by two solid, adiabatic side walls of height $0.05H$ , with a temperature difference $\Delta T$ between the top and bottom boundaries. To model a liquid metal of density $\rho$ , kinematic viscosity $\nu$ , thermal diffusivity $\alpha$ , thermal expansion coefficient $\beta$ , the flow is assumed Newtonian and since temperature gradients remain moderate, its motion is assumed well described by the Oberbeck–Boussinesq approximation (Oberbeck Reference Oberbeck1879; Boussinesq Reference Boussinesq1903; Chandrasekhar Reference Chandrasekhar1961). This leads to the following non-dimensional governing equations:

(2.1) \begin{align} \frac {\partial \textbf{u}}{\partial t}+ (\textbf{u} \cdot \nabla) \textbf{u} + \nabla p & = RaPr T \textbf{e}_y+Pr \nabla ^2 \textbf{u}, \\[-12pt] \nonumber \end{align}

(2.2) \begin{align} \frac {\partial T}{\partial t}+ (\textbf{u} \cdot \nabla) {T} & = \nabla ^2 T, \\[-12pt] \nonumber \end{align}

(2.3) \begin{align} \nabla \cdot \textbf{u} & = 0, \\[6pt] \nonumber \end{align}

where $\textbf{u}=(u,v,w)$ is the velocity vector field, $T$ is the temperature field, $t$ is the time and $\mathbf g=-g\mathbf e_y$ is the gravitational acceleration. The modified pressure $p$ includes the buoyancy term that accounts for the reference temperature $T_0$ at density $\rho _0$ (Chandrasekhar Reference Chandrasekhar1961; Tritton Reference Tritton1988). The above set of equations are non-dimensionalised using length $H$ , velocity $\alpha /H$ , time $H^2/\alpha$ , pressure $\rho _0(\alpha /H)^2$ and temperature $\Delta T$ . The equations feature two governing non-dimensional groups: the Prandtl number

(2.4) \begin{equation} Pr = \frac {\nu }{\alpha }, \end{equation}

which we fixed to $0.02$ , a typical value of aluminium alloys, and the Rayleigh number $Ra$ , defined as

(2.5) \begin{equation} Ra = \frac {\beta g \Delta T H^3}{\nu \alpha }, \end{equation}

which controls the intensity of buoyancy forces relative to viscous forces. A free-slip boundary condition is applied to the upper boundary at $y=1$ with fluid being poured homogeneously across the boundary at an imposed temperature $\Delta T$ . These boundary conditions are expressed as

(2.6) \begin{eqnarray} \frac {\partial }{\partial y} \mathbf u\times \mathbf e_y=0, \qquad \mathbf u\cdot \mathbf e_y=RePr, \qquad T(y=1)=1, \end{eqnarray}

where $Re$ is the mass flux Reynolds number based on the dimensional feeding velocity $u_0$ :

(2.7) \begin{equation} Re = \frac {u_0H}{\nu }. \end{equation}

At the lower boundary $\mathcal S$ , the solidification front is represented by solid, porous boundary conditions,

(2.8) \begin{eqnarray} \mathbf u_{\mathcal{S}}\times \mathbf e_y=0, \qquad \mathbf u_{\mathcal{S}}\cdot \mathbf e_y=RePr, \qquad T_{\mathcal{S}}=0, \end{eqnarray}

such that the flux of fluid pulled at $\mathcal{S}$ exactly cancels the flux of mass coming from the inlet. Impermeable, no-slip boundary conditions for the velocity field and insulating boundary conditions for the temperature field are imposed at the side walls (see figure 1). These boundary conditions for the temperature field at these side walls ensure that the temperature field remains continuous along the entire periphery of the domain. The pressure field in our calculation is obtained by solving the Poisson equation derived by taking the divergence of (2.1). A consistent Neumann boundary condition for pressure (Karniadakis & Israeli Reference Karniadakis, Israeli and Orszag1991) is applied at $y = 1$ , on the lower boundary $\mathcal{S}$ , and on the side walls. In the third direction $\textbf{e}_z$ , the domain’s infinite extension is represented by periodic boundary conditions for all flow fields.

2.2. Direct and adjoint perturbation equations

The purpose of this work is to suppress linear instabilities with an actuation specifically designed to stifle the growth of linear perturbations. Since this actuation will be based on the linear dynamics of this perturbation, we first need the establish the set of equations that govern its linear dynamics. From Kumar & Pothérat (Reference Kumar and Pothérat2020), linear perturbations grow from a class of equilibrium solutions of (2.1)–(2.3) and boundary conditions (2.6) and (2.8) that are planar and invariant along $\textbf{e}_z$ , hence of the form $\textbf{U}(x,y)$ , $\bar {T}(x,y)$ . Baroclinicity due to the isothermal condition along the solidification front precludes any purely diffusive thermal equilibrium, so $\textbf{U}(x,y)$ is never homogeneously zero and both $\textbf{U}(x,y)$ and $\bar {T}(x,y)$ must be found as a fully nonlinear two-dimensional (2-D) solution of the equations. It follows that perturbations to this equilibrium $\textbf{q}^\prime (x,y,z,t)=(\textbf{u},T)^\top - (\textbf{U}, \bar {T})^\top =(\textbf{u}^\prime (x,y,z,t), T^\prime (x,y,z,t))^\top$ , are governed by the linear system of equations governing the evolution of infinitesimal perturbations,

(2.9) \begin{equation} \frac {\partial \textbf{q}^\prime }{\partial t} = \mathcal{L} \textbf{q}^\prime, \end{equation}

where

(2.10) \begin{equation} \mathcal{L}\textbf{q}^\prime = \begin{pmatrix} - (\textbf {U} \cdot \nabla) \textbf {u}^\prime -( \nabla \textbf{U}) \cdot \textbf {u}^\prime -\nabla p^\prime + RaPrT^\prime \textbf{e}_y+Pr{\nabla }^2\textbf{u}^\prime \\ -\textbf{U} \cdot \nabla T^\prime - (\nabla \bar {T})\cdot \textbf{u}^\prime + \nabla ^2 T^\prime \\ \end{pmatrix}. \end{equation}

Since $p^\prime$ is determined by the constraint $\nabla \cdot \textbf{u}^\prime = 0$ , it is therefore not included in the state vector $\textbf{q}^\prime$ . We will refer to (2.9) as the direct perturbation equation, and to $\mathcal{L}$ as the direct linear operator. The boundary conditions for the base flow are the same as those for the main variables. As a result, the perturbation variables satisfy the homogeneous counterparts of the boundary conditions associated with the base flow. Since the base flow is invariant along $\textbf{e}_z$ and $\mathcal{L}$ does not explicitly depend upon time, the perturbation can be written as a linear combination of normal modes of the form,

(2.11) \begin{equation} \textbf {q}^\prime (x,y,z,t) = \hat {\textbf {q}}(x,y) e^{ikz + \lambda t}, \end{equation}

where $k$ is the wavenumber along the homogenous direction $\textbf{e}_z$ and $\lambda = \sigma \pm i\omega$ contains the growth rate, $\sigma$ and frequency $\omega$ . The growth rate, frequency and the wavenumber of the most dominant mode $\hat {\textbf {q}}(x,y)=(\hat {u}, \hat {v}, \hat {w}, \hat {T})^\top$ (also referred to as the direct mode) are found by solving the eigenvalue problem for $\lambda$ that result from (2.9) and (2.11). This is done numerically by means of a time-stepper method (Barkley, Blackburn & Sherwin Reference Barkley, Blackburn and Sherwin2008). Here the eigenmodes are normalised such that $\|\hat {\textbf{q}}\|_2=1$ , where $\|\cdot \|_2$ denotes the standard $l^2$ vector norm of $4 \times N_e \times (N+1)^2$ values that make up $\hat {\textbf {q}}$ . In this context, $N_e$ refers to the number of quadrilateral elements. For the details of the eigenvalue solver we refer the reader to § 2.2 of our previous work (Kumar & Pothérat Reference Kumar and Pothérat2020).

Next, we need to work out the form of the actuation that is best suited to suppress the growth of individual normal modes. The idea we pursue relies on the ideas of Giannetti & Luchini (Reference Giannetti and Luchini2007), who showed that the receptive regions of the direct linear modes are mapped by the adjoint eigenmodes of the same linearised equations. We, therefore, need to construct the adjoint operator $\mathcal{L}^*$ with respect to the time-averaged inner product relevant to the problem (Mao, Blackburn & Sherwin Reference Mao, Blackburn and Sherwin2015),

(2.12) \begin{equation} (\textbf{a}, \textbf{b}) = \int _\Omega \textbf{a}\cdot \textbf{b} \, \text{d}V, \qquad \langle \textbf{c},\textbf{d} \rangle = \int _0^\tau (\textbf{c}, \textbf{d}) \, \textrm{d}t, \end{equation}

where $\textbf{a}$ and $\textbf{b}$ are time-averaged vector fields defined on the fluid domain $\Omega$ , while $\textbf{c}$ and $\textbf{d}$ are time-dependent vector fields defined on $\Omega$ and time domain $[0, \tau]$ . The adjoint operator is then defined by the relation

(2.13) \begin{equation} \langle \textbf{q}^*, (-\partial _t +\mathcal{L})\textbf{q}^\prime \rangle - \langle (\partial _t+ \mathcal{L}^*)\textbf{q}^*, \textbf{q}^\prime \rangle =0, \end{equation}

and thus satisfies

(2.14) \begin{equation} -\frac {\partial \textbf{q}^*}{\partial t} = \mathcal{L}^* \textbf{q}^*, \end{equation}

where $\textbf{q}^*(x,y,z,t) = (u^*, v^*, w^*, T^*)^\top$ represents adjoint variables. The expression of $\mathcal L^*$ is readily derived by applying integration by parts to the first term in (2.13),

(2.15) \begin{equation} \mathcal{L}^*\textbf{q}^* = \begin{pmatrix} (\textbf{U} \cdot \nabla) \mathbf{ u^*} - (\nabla \textbf{U})^\top \cdot \mathbf{ u^*} -(\nabla \bar {T})T^* -\nabla p^* + Pr{\nabla }^2\textbf{u}^* \\ (\textbf{U} \cdot \nabla) {T^*}+ RaPr\, v^*+\nabla ^2{T}^* \\ \end{pmatrix}, \end{equation}

with $\nabla \cdot \textbf{u}^* = 0$ . Similarly, the adjoint variables satisfy adjoint boundary conditions imposed by (2.13). These are identical to the boundary conditions satisfied by the direct variables, except for the $\textbf{e}_x$ and $\textbf{e}_z$ components of the velocity field at the upper free surface, that must satisfy Robin conditions (Barkley et al. Reference Barkley, Blackburn and Sherwin2008)

(2.16) \begin{equation} \textbf{e}_y \cdot \nabla u^* - Re \, u^*=0, \qquad \textbf{e}_y \cdot \nabla w^* -Re \, w^*=0. \end{equation}

Like the direct variables, the adjoint variables are decomposed into the normal modes but obtained as the solution of the adjoint eigenvalue problem, instead of the direct one so this time, the same eigenvalue solver yields the adjoint mode $\hat {\textbf {q}}^*(x,y)=(\hat {u}^*, \hat {v}^*, \hat {w}^*, \hat {T}^*)^\top$ . It follows from the biorthogonality property that the eigenvalue of the adjoint mode is the complex conjugate of that of the direct mode (Salwen & Grosch Reference Salwen and Grosch1981). Therefore, the magnitude of the growth rate and frequency of the adjoint and direct modes are identical when the real and imaginary parts of the complex eigenvalue are equated.

2.3. Receptivity and sensitivity

Giannetti & Luchini (Reference Giannetti and Luchini2007) showed that the adjoint field represented Green’s function for the receptivity problem. Therefore, the adjoint equations can be used to evaluate the effects of any external actuation or a generic initial condition on the leading eigenmode of the direct problem (2.9). This property forms the basis for the analysis in § 4, where we seek to identify the regions where applying an actuation would most effectively affect the growth of unstable modes. Our intention is to control the growth of the perturbation from an initial condition where the unstable mode has already grown measurably, but sufficiently little to remain within the confines of the linear approximation. For this purpose, we determine the topology of the actuation on the adjoint eigenmode. A common alternative approach is to force the system in such a way as to modify the underlying operator such that the eigenvalue associated with the leading eigenmode remains in the stable region, ignoring the initial condition. The optimal actuation for this purpose is given by the sensitivity map. The relative shift in the eigenvalue, and therefore the growth rate associated with the direct and the adjoint mode incurred by acting at any given location of the flow is given by the sensitivity map (Giannetti & Luchini Reference Giannetti and Luchini2007; Giannetti et al. Reference Giannetti, Camarri and Luchini2010; Qadri, Mistry & Juniper Reference Qadri, Mistry and Juniper2013)

(2.17) \begin{equation} S_{ij}(x,y)=\frac {\hat {\textbf {q}}_i \hat {\textbf {q}}^*_j }{\int _\Omega \hat {\textbf {q}}^\top \hat {\textbf {q}}^* \, \text{d}x\text{d}y}.\end{equation}

Note that the sensitivity tensor above is based on feedback localised in space, of the form ${\mathbf C_0} \delta ({\mathbf x} - {\mathbf x_0}) \hat {\textbf {q}}$ . Here, $\mathbf C_0$ denotes a constant coefficient matrix, $\mathbf x_0$ indicates the position where the feedback acts and $\delta ({\mathbf x} - {\mathbf x_0})$ denotes the Dirac delta function. Tensor $S_{ij}$ determines the relative local intensity of the feedback of individual component $\hat {\textbf {q}}^*_j$ onto the individual component $\hat {\textbf {q}}_i$ of the eigenmode. This quantity also locates the regions of the flow acting as wavemakers. In particular, since ${\mathbf q}^\prime$ and $\mathbf q^*$ contain all three components of velocity and the temperature, the knowledge of $S_{ij}$ indicates whether the actuation should be of a thermal or mechanical nature and if mechanical, which component of the velocity (or combination of all four components including temperature) is most efficient at altering the growth of the unstable mode. Both the direct and adjoint problems being linear, amplitudes are relative so we may further normalise the direct and adjoint modes by choosing (Giannetti & Luchini Reference Giannetti and Luchini2007)

(2.18) \begin{equation} \int _\Omega \hat {\textbf {q}}^\top \hat {\textbf {q}}^* \, \text{d}x\text{d}y = 1, \end{equation}

so that the sensitivity tensor is simply expressed as $S_{ij}=\hat {\textbf {q}}_i \hat {\textbf {q}}^*_j$ .

At this point, we reiterate that our study considers the feedback on the perturbed field. This differs from the approach of Marquet et al. (Reference Marquet, Sipp and Jacquin2008), which examines the effect of forcing on the base flow. The objective of our study is to apply a forcing based on the adjoint mode (receptivity) to suppress instability. In general, the forcing influences both the base flow and the perturbation and this effect can be analysed in detail by studying the sensitivity to base flow modification as Marquet et al. (Reference Marquet, Sipp and Jacquin2008) do or by validating the approach using nonlinear direct numerical simulation (DNS) as we do in this paper. To summarise the difference between these two strategies in a nutshell, our strategy consists in applying a forcing that suppresses the instability before either the perturbation or the forcing has sufficiently grown to affect the base flow. The strategy proposed by Marquet et al. (Reference Marquet, Sipp and Jacquin2008), by contrast, is to modify the base flow to prevent the perturbation from growing at all.

Figure 2. Streamlines of the steady 2-D base flow and temperature field for the simulation cases (a) C1, (b) C2, (c) C3 and (d) C4.

2.4. Methodology and choice of parameters

For the reminder of the paper, we shall proceed as follows to find and assess the actuation best suited to damp the growth of linear instabilities. First, the steady 2-D base flow solutions are obtained using DNS of (2.1)–(2.3) together with the associated boundary conditions. Second, the LSA of the 2-D base flow against three-dimensional (3-D) perturbations is carried out by solving the eigenvalue problem for operator $\mathcal L$ . These first two steps were previously carried out over an extensive range of governing parameters and wavenumbers $k$ in Kumar & Pothérat (Reference Kumar and Pothérat2020). Here, we repeat the same approach for flows that are weakly supercritical so as to focus on cases where the instability is driven by a small number of unstable modes. The idea behind this strategy is that even in the fully nonlinear evolution, if the instability remains driven by a small-enough number of modes, it may be enough to prevent the growth of the most unstable of them to stop the growth of instabilities altogether. This approach is not expected to be successful if the base flow is destabilised by a broad spectrum of fast-growing unstable modes, as may be the case in more strongly supercritical cases. On this basis, we select four typical weakly supercritical cases illustrating the different instability mechanisms identified in this previous work.

1. (i) Simulation case C1 ( $Re = 0; Ra=7 \times 10^3$ ). This case corresponds to a purely convective base flow with zero inlet (and outlet) mass flux, i.e. no fluid crosses the boundaries of the flow domain (see figure 2 a). The flow becomes unstable to a travelling wave at $Ra_c=5.975 \times 10^3$ , through a supercritical Hopf bifurcation. The corresponding branch in the complex eigenmode spectrum is labelled type II.

2. (ii) Simulation case C2 ( $Re = 50; Ra=7 \times 10^3$ ). This case features an inflow through the upper boundary as shown in figure 2(b). The flow becomes unstable to a type II travelling wave through a supercritical Hopf bifurcation.

3. (iii) Simulation case C3 ( $Re = 100; Ra=4 \times 10^4$ ). This case is similar to C2 but with more intense inflow. The base flow is presented in figure 2(c). At $Ra_c=3.621 \times 10^4$ , it becomes unstable to a leading mode from a different branch (type I), that is still oscillatory but arises out of a subcritical Hopf bifurcation.

4. (iv) Simulation case C4 ( $Re = 150; Ra=8 \times 10^4$ ). Here the flow is mostly driven by the inflow (see figure 2 d). The instability corresponds to a further branch of the eigenvalue spectrum labelled type III and yields a non-oscillatory mode through a supercritical pitchfork bifurcation.

Details of the above selected parameters are tabulated in table 1.

Table 1. Parameters of the weakly supercritical cases, C1, C2, C3 and C4, were selected for the analysis. Here $(Ra/Ra_c)-1$ represents the level of criticality, where $Ra_c$ is the critical Rayleigh number; $N_U$ represents the number of unstable modes; $k$ represents the most unstable mode. Note that $Ra_c$ for the each case is obtained from table 2 of Kumar & Pothérat (Reference Kumar and Pothérat2020).

Additionally to this previous analysis, we now need to calculate adjoint modes to perform the RSA. These are obtained by solving the eigenvalue problem for the operator $\mathcal{L}^*$ using the same method as for the LSA. Third, the evolution of the normal mode targeted for suppression is calculated using the linearised Navier–Stokes equations (2.9), to which a forcing term based on the adjoint eigenmode is added on the left-hand side. We use two types of actuation to that effect. As the unstable mode is suppressed by the forcing, its amplitude varies, and so does the amplitude of the actuation. Otherwise, once the unstable mode is suppressed, the actuation would act as an external force taking the flow away from its equilibrium. Hence, to assess the ability of an actuation based on the receptivity map to suppress the unstable mode, we use a forcing of the form

(2.19) \begin{equation} \textbf{f}(x,y,z,t) = \alpha (t)\Re \{A\hat {\textbf {q}}^*(x,y)\exp {[ \sigma t + i(kz + \omega t + \phi)]}\}, \end{equation}

where $\alpha (t)$ represents the amplitude of the unstable mode normalised by its initial value. In practice, adapting the forcing to the amplitude in real time would require sensing and processing capable of extracting the evolution of the unstable mode instantaneously, which, in the harsh environment of continuous casting, is impractical. Instead, it is much easier to set a threshold on the sensor output above which an actuation of constant amplitude is applied. For this purpose, we use a simpler form of actuation,

(2.20) \begin{equation} \textbf{f}(x,y,z,t) = \Re \{A\hat {\textbf {q}}^*(x,y)\exp {[ \sigma t + i(kz + \omega t + \phi)]}\}. \end{equation}

In both cases, $A$ and $\phi$ represent the initial real amplitude and real phase of the forcing, respectively, relative to the unstable mode. Note that our strategy is not to choose an actuation intended to manipulate the eigenvalue associated with the leading eigenmode. This well-established technique entails applying a force whose linear dependence on the unstable eigenmode shifts the leading eigenvalue of the linear evolution operator to the stable region. The topology of the optimal forcing map is, in this case, provided by the structural sensitivity, or base-flow sensitivity maps (Giannetti & Luchini Reference Giannetti and Luchini2007; Marquet et al. Reference Marquet, Sipp and Jacquin2008). Instead, we chose to apply a predetermined force, i.e. that does not depend linearly on variable $\hat {\textbf {q}}$ . That actuation is designed to suppress the amplitude of the unstable mode from a small but finite initial value. Since the actuation does not linearly depend on $\hat {\textbf {q}}$ , the evolution equation is not strictly linear, its linear part does not change and neither does the leading eigenvalue. The actuation f has four components. The first three components represent mechanical forces, and the last component represents thermal forces. This step serves two purposes. First, it acts as a validation of the strategy. Since the forcing is derived from linearised adjoint equations, that ignore nonlinear interactions, it must at least be effective at damping the direct linear mode, to carry any hope of preventing the full nonlinear growth of the instability. Second, building the forcing on the topology of the adjoint mode involves a choice of amplitude and phase relative to the direct modes. To find out the optimal values of both, we carry out a series of linear simulations where they are varied.

Finally, with the knowledge of the optimal amplitude and phase of the actuation, we test the suppression of instability with its full nonlinear dynamics through 3-D DNS (the detail of individual simulations is given in § 5).

2.5. Numerical set-up

The methodology outlined in the previous section involves four types of numerical computations: 2-D (nonlinear) DNS; direct and adjoint eigenvalue problems; 3-D evolution of individual eigenmodes through the direct linearised equations; 3-D (nonlinear) DNS. The solution of the direct and adjoint eigenvalue problems are found by means of the time-stepping method implemented and tested in detail in Kumar & Pothérat (Reference Kumar and Pothérat2020). The novelty compared with this previous work is the solution of the adjoint eigenvalue problem, which was validated by making sure the eigenvalues obtained from both the LSA and RSA yielded the same results down to machine precision.

The nonlinear governing equations (2.1)–(2.3), the direct perturbation equation (2.9) and the adjoint equation (2.14) are solved using the spectral-element code Nektar++ (Cantwell et al. Reference Cantwell2015; Moxey et al. Reference Moxey2020). We adopted a spectral-element discretisation in the $x{-}y$ plane with a mesh consisting of 348 quadrilateral elements. For the 3-D simulations, we used a Fourier-based spectral method (Bolis et al. Reference Bolis, Cantwell, Moxey, Serson and Sherwin2016) for discretisation in the $\textbf{e}_z$ direction. The computational domain extends along $\textbf{e}_z$ by $2\pi$ . A third-order implicit–explicit method (Vos et al. Reference Vos, Eskilsson, Bolis, Chun, Kirby and Sherwin2011) is used for time stepping. For all four types of numerical calculations, the time step $\Delta t$ was kept constant so that the maximum local Courant number $C_{max }$ remained below unity everywhere in the domain at all times. The numerical implementation is described and tested in detail in Kumar & Pothérat (Reference Kumar and Pothérat2020). Figure 1 shows the details of a 2-D $x{-}y$ mesh generated using the Gmsh package (Geuzaine & Remacle Reference Geuzaine and Remacle2009) with polynomial order $N = 3$ as an example of spatial–spectral discretisation used in the $(x,y)$ plane. Elements at the edges are more densely packed than in the bulk, with a ratio of four between the edge sizes of the largest and the smallest elements. On each element, the flow variables are projected onto the polynomial basis represented at Gauss–Lobatto–Legendre points. As in our previous work, we perform a convergence test on the polynomial order for the leading eigenvalue for each case to ensure the solution is independent of the spectral order. For example, the leading eigenvalue computed for the simulation case C4 with the polynomial order $N=8$ differs by less than 0.04 % from the calculation performed with $N=9$ . Convergence of the leading eigenvalue with the polynomial order is presented in table 2. Similar convergence tests have been conducted for other simulation cases. For all 3-D DNS calculations, we have retained 32 Fourier modes along $\textbf{e}_z$ , which are deemed adequate based on our previous convergence test conducted in Kumar & Pothérat (Reference Kumar and Pothérat2020).

Table 2. We examine the relationship between the leading eigenvalues and the polynomial order $N$ . The leading eigenvalues are computed on the mesh at $Re=150$ , $Ra=8 \times 10^4$ and $k=6$ . The relative error is calculated with respect to the case of the highest polynomial order ( $N=9$ ).

4.1. Methodology

In the first instance, we seek the linear response of the system. This step acts as a validation, as the sensitivity analysis already provided us with an indication of the linear response we should expect. Two very important aspects still remain to be clarified by calculating the linear evolution of the leading eigenmode under the effect of the actuation: first, the time scale of the response and the duration of the effect are not considered in the receptivity analysis; second, the receptivity does not specify how the relative amplitude and phase of the receptive mode affect the linear response. Indeed, in (2.19) and (2.20), the amplitude and phase of the actuation are both relative to the leading eigenmode. Since the receptivity mode is the solution of a homogeneous linear problem, none of them is specified. On the other hand the relative amplitude and phase of the actuation can be expected to greatly influence the system’s response to it. For these two reasons, and to identify the combination of relative phase and amplitude that optimises the suppression of the instability, we run a series of linear simulations by adding a source term representing the receptivity-based actuation on the right-hand side of (2.9):

(4.1) \begin{equation} \frac {\partial \textbf{q}^\prime }{\partial t} = \mathcal{L} \textbf{q}^\prime + \textbf{f}. \end{equation}

We conduct a parametric study for the actuation of constant amplitude (2.20), varying the phase $\phi$ from $0$ to $2\pi$ in increments of $\pi /10$ and the amplitude between 0.01 and 0.5 for most of the cases. Then, we conduct a single linear simulation with adaptive actuation amplitude (2.19), using the most effective combination of $A$ and $\phi$ found in the parametric analysis, to assess whether the receptivity-based actuation indeed fully damps the unstable mode asymptotically. Each linear simulation is initiated using the unstable mode obtained from the LSA as the initial condition with amplitude such that the normalisation condition (2.18) is satisfied, from which both $A$ and $\phi$ are fixed. In the linear simulation, normal modes are decoupled from each other so the time-dependent solution obtained by initialising the solution with a single normal mode reflects the evolution of that particular mode only. As such, linear simulations provide a direct measure of the ability of the actuation based on the receptivity mode to affect the evolution of the unstable mode. To quantitatively asses this effect, we monitor the time-dependent energy of the mode, $(\textbf{q}^\prime, \textbf{q}^\prime)$ . The lowest value $(\textbf {q}^\prime, \textbf{q}^\prime)_{min }$ reached by this quantity gives an indication of the damping achieved by the actuation, and the time of occurrence of this minimum measures the time scale over which this damping is achieved, denoted as $t_f$ . Going back to the example of the continuous casting of alloys as one of the motivations for this work, the unstable mode may not need to be suppressed indefinitely. Instead, maintaining the unstable mode to a low level for the entire duration of the casting operation, which is finite, would be sufficient to ensure it does not impact the final quality of the alloy. Hence the importance of $t_f$ is both fundamental and practical.

Figure 4. For the simulation case C1 ( $Re=0$ , $Ra=7 \times 10^3$ and $k=6$ ) (a) Minimum value of the strength of perturbation $(\textbf{q}^\prime, \textbf{q}^\prime)_{min }$ as a function of the amplitude $A$ and the phase $\phi$ of the linear receptivity-based actuation. (b) Dependence of the time $t_f$ at which $(\textbf{q}^\prime, \textbf{q}^\prime)_{min }$ occurs on the amplitude $A$ and phase $\phi$ . (c) The strength of perturbation $(\textbf{q}^\prime, \textbf{q}^\prime)$ as a function of time $t$ for both unforced and forced cases. The inset represents the evolution of both constant amplitude actuation and adaptive actuation from $t=0$ to $t \simeq t_f$ .

4.2. Analysis of cases C1–C4

We shall now examine the outcome of this approach in each of the C1–C4 cases, defined in § 2.4. Figure 4(c) illustrates four examples of the evolution of $(\textbf{q}^\prime, \textbf{q}^\prime)$ for case C1. The blue curve represents a simulation with no actuation, i.e. one with $A=0$ amplitude. As predicted by LSA, the instability grows exponentially with oscillations of frequency determined by the imaginary part of the mode’s eigenvalue. The red curve represents the case with receptivity-based actuation of amplitude $A=0.5$ and phase $\phi =\pi$ . In this case, $(\textbf{q}^\prime, \textbf{q}^\prime)$ decreases over time and reaches its minimum value at $t_f=3.65$ . For $t\gt t_f$ , $(\textbf{q}^\prime, \textbf{q}^\prime)$ increases, and crosses the blue curve at $t=7.16$ . Hence, in this particular case, the actuation results in an effective suppression of the instability until $t=t_f$ . Now, to determine the optimal value of $A$ and $\phi$ , i.e. values that achieve the smallest value of $(\textbf{q}^\prime, \textbf{q}^\prime)_{min }$ , we vary the phase in the range of $0$ to $2\pi$ , and for amplitudes between 0.07 and 1.0. The grid map of the corresponding values of $(\textbf{q}^\prime, \textbf{q}^\prime)_{min }$ is shown in figure 4(a) for case C1. It turns out that $A=0.5$ and $\phi =\pi$ are the optimal values of amplitude and phase with $(\textbf{q}^\prime, \textbf{q}^\prime)_{min } = 0.06$ . Simulation conducted with adaptive actuation (2.19) shows that the unstable mode decays very similarly to the case of constant amplitude actuation until $t\simeq t_f$ , but continues to decay asymptotically to zero as $t \to \infty$ . This shows that the receptivity-based actuation indeed suppresses the unstable mode completely.

Figure 5. For the simulation case C2 ( $Re=50$ , $Ra=7 \times 10^3$ and $k=6$ ). (a) Minimum value of the strength of perturbation $(\textbf{q}^\prime, \textbf{q}^\prime)_{min }$ as a function of the amplitude $A$ and the phase $\phi$ of the linear receptivity-based actuation. (b) Dependence of the time $t_f$ at which $(\textbf{q}^\prime, \textbf{q}^\prime)_{min }$ occurs on the amplitude $A$ and phase $\phi$ . (c) The strength of perturbation $(\textbf{q}^\prime, \textbf{q}^\prime)$ as a function of time $t$ for both unforced and forced cases. The inset represents the evolution of both constant amplitude actuation and adaptive actuation from $t=0$ to $t \simeq t_f$ .

Figure 6. For the simulation case C3 ( $Re=100$ , $Ra=4 \times 10^4$ and $k=4$ ). (a) Minimum value of the strength of perturbation $(\textbf{q}^\prime, \textbf{q}^\prime)_{min }$ as a function of the amplitude $A$ and the phase $\phi$ of the linear receptivity-based actuation. (b) Dependence of the time $t_f$ at which $(\textbf{q}^\prime, \textbf{q}^\prime)_{min }$ occurs on the amplitude $A$ and phase $\phi$ . (c) The strength of perturbation $(\textbf{q}^\prime, \textbf{q}^\prime)$ as a function of time $t$ for both unforced and forced cases. The inset represents the evolution of both constant amplitude actuation and adaptive actuation from $t=0$ to $t \simeq t_f$ .

Further calculations were performed for cases C2 and C3, which, despite corresponding to different instability branches, return similar results, presented in figures 5 and 6. As compared with case C1, case C2 has an optimal amplitude that is twice as large, with a phase shift of $1.7\pi$ radians. Therefore, the optimal values are $A=1.0$ and $\phi =1.7\pi$ . The value of $t_f$ is slightly reduced to 2.57. On the other hand, for case C3, the optimal amplitude is reduced to 0.1, and the optimal phase is shifted to $1.5\pi$ . The value of $t_f$ , however, is slightly increased to 3.93, which implies that the instability can be attenuated for a slightly longer period of time. As in case C1, the adaptive actuation (2.19) leads to a decay similar to the constant amplitude actuation until $t\simeq t_f$ , but it continues to decay asymptotically to zero as $t \to \infty$ , so in cases C2 and C3 too, the receptivity-based actuation fully suppresses the unstable mode.

Figure 7. For the simulation case C4 ( $Re=150$ , $Ra=8 \times 10^4$ , and $k=6$ ). (a) Minimum value of the strength of perturbation $(\textbf{q}^\prime, \textbf{q}^\prime)_{min }$ as a function of the amplitude $A$ of the linear receptivity-based actuation. Inset represents the amplitude $A=-0.03$ which corresponds to the lowest value of $(\textbf{q}^\prime, \textbf{q}^\prime)_{min }$ . (b) Dependence of the time $t_f$ at which $(\textbf{q}^\prime, \textbf{q}^\prime)_{min }$ occurs on the amplitude $A$ . (c) The strength of perturbation $(\textbf{q}^\prime, \textbf{q}^\prime)$ as a function of time $t$ for both unforced and forced cases. The inset represents the evolution of both constant amplitude actuation and adaptive actuation from $t=0$ to $t \simeq t_f$ .

Case C4 differs from the other cases in that the leading eigenmode and its adjoint, which represents receptivity, are non-oscillatory. Consequently, the actuation is also non-oscillatory and we do not need to determine the optimal phase, only the optimal amplitude. In the absence of a phase, the possibility still remains to reverse the actuation, which we incorporated into the sign of $A$ . The variations of $(\textbf{q}^\prime, \textbf{q}^\prime)_{min }$ with $A$ are shown in figure 7(a). For $A\geq 0$ , $(\textbf{q}^\prime, \textbf{q}^\prime)_{min }$ remains very close to the initial value of $(\textbf{q}^\prime, \textbf{q}^\prime)$ so no suppression is achieved. For $A\lt 0$ , by contrast, we observe a sudden drop in the value of $(\textbf{q}^\prime, \textbf{q}^\prime)_{min }$ at very small amplitudes, with a minimum value at $A=-0.03$ (see inset of figure 7 a).

Figure 7(c) shows the evolution of $(\textbf{q}^\prime, \textbf{q}^\prime)$ over time for both the unforced and forced cases with $A=-0.03$ . It appears that the actuation achieves a suppression of the non-oscillatory mode up to time $5.45$ . Here again, the adaptive actuation (2.19) suppresses the unstable mode completely and acts over the same time scale $t_f$ as the constant amplitude actuation.

4.3. The salient features of the linear response to an actuation of constant amplitude

In all cases, the receptivity-based actuation with adaptive amplitude suppresses the unstable mode asymptotically. This result shows that the theoretical foundations for the receptivity-based instability control are sound. However, the actuation with constant amplitude is more useful in practical situations, but induces a more complex response, which we now discuss in more detail. All four cases show strong similarities, despite the differences in the nature of the unstable modes. First, the linear actuation has a stabilising effect up to a finite time $t_f$ , effectively suppressing the instability. Beyond this point, the unstable mode grows again and at a faster pace than when unforced. This means that for $t\gt t_f$ , the linear actuation becomes destabilising. The reason is that, since the base flow is fixed in the linear equations we simulate, the actuation’s only effect is to suppress the unstable mode. Once the mode is suppressed, the system is back to the equilibrium point defined by the steady base flow. The actuation then acts as an additional force, moving the system away from this equilibrium. Nonlinear simulations are required to verify whether the base flow is indeed affected by the actuation and whether the phenomenology observed in the linear framework is indeed valid. Table 3 shows that $t_f$ remains of the order of the inverse of the growth rate of the unstable mode, as the ratio $\sigma t_f$ remains of the order of unity for all suppression-optimising cases. The term ‘suppression-optimising’ refers to seeking the amplitude and phase of the forcing that achieve the maximum reduction in mode amplitude. The reduction in mode amplitude achieved at $t=t_f$ is of approximatively two orders of magnitudes (see table 3) for the suppression-optimising cases. This suggests that the actuation is more effective when applied for a finite time $t_f$ , then stopped until the unstable mode regrows to its initial amplitude (i.e. during a time $\sim \sigma ^{-1}$ ), at which point the procedure can be repeated. Whether this strategy is realistic depends on how this pattern is affected by nonlinearities.

Table 3. The normalised value of $t_f$ with respect to the growth rate $\sigma$ for various suppression-optimising and time-optimising simulation cases. Comparison of $(\textbf{q}^\prime, \textbf{q}^\prime)_{min }$ relative to the initial value of $(\textbf{q}^\prime, \textbf{q}^\prime)$ is shown for each suppression-optimising and time-optimising simulation cases.

This pattern also raises the question of whether instead of applying suppression-optimising, it may be beneficial to seek configurations that achieve the largest value of $t_f$ so as to maximise the time during which the suppression is effective. The term ‘time-optimising’ is introduced to denote the search for amplitude and phase settings that maximise the value of $t_f$ . The corresponding optimal values are shown on figures 4(b), 5(b), 6(b) and 7(b). Unsurprisingly, actuation parameters optimising the suppression and $t_f$ differ. We shall compare how both optimals perform on the nonlinear response in the next section.

In addition to identifying optimal suppression parameters, our results highlight that certain combinations of amplitude and phase can instead enhance instability. For instance, as shown in figure 4(c), the case $A=1$ and $\phi =0$ leads to the most destabilising effect for $Re=0$ , and similarly, as shown in figure 6(c), for $ Re=100$ . In contrast, for $Re=50$ (figure 5 c), the configuration $A=1$ and $\phi =\pi$ is more destabilising than $ A=1$ and $ \phi =0$ . This observation underscores the importance of carefully selecting both the amplitude and phase in actuation design, as an inappropriate choice can inadvertently amplify perturbations rather than suppressing them.

5.1. Methodology

The linear simulations confirmed that an actuation designed from the topology of the receptivity, with a frequency matched to the leading eigenmode was indeed capable of stifling the linear mechanism responsible for the growth of the leading eigenmode. We also identified the relative amplitude and phase that maximised its suppression. In reality, the finite amplitude of the perturbation, whether subject to the actuation or not, activates a nonlinear response of the system. Although nonlinearities may not measurably deflect the perturbation’s evolution in its early stages, when its amplitude is still small, it is likely to govern its dynamics at longer time scales, especially at the point where the actuation ceases to be efficient. For this reason, the nonlinear effects are essential to determine how effective the actuation may be at suppressing instabilities, and how long it may remain so.

Including the nonlinear dynamics, however, raises several technical difficulties. First, the optimal amplitude and phase were determined without any consideration of nonlinearity, so we may question whether these remain optimal for the nonlinear evolution. The answer lies partly in how the suppression is assessed. Following the approach we adopted in the linear study, we sought the point of lowest amplitude for the leading eigenmode and the time of this minimum. Since the evolution from an infinitesimal perturbation to that point involves only very small amplitudes, it is legitimate to assume that nonlinearities would play little role there and that, consequently, the linear and the nonlinear evolutions would remain very similar during this phase (Drazin & Reid Reference Drazin and Reid2004). One should, however, keep in mind that if phase and amplitude were sought so as to optimise the saturated state of the evolution, the full nonlinear equations would need to be included in the optimisation process (Pringle, Willis & Kerswell Reference Pringle, Willis and Kerswell2012). However, our specific purpose of understanding the influence of the nonlinear effects is better served by keeping a similar approach to the linear study. Additionally, since the linear study showed that the growth of the leading eigenmode was only effectively suppressed until $t=t_f$ , we shall only apply the actuation up to that time, and let the flow evolve freely after that time.

Second, the definition of the initial conditions for the nonlinear problem differs from the linear one: the linear equations indeed return the same evolution regardless of the amplitude of the initial condition (up to a multiplicative factor), and only the relative amplitude and phase of the actuation mattered. To transpose our approach into the nonlinear framework, we shall also specify the amplitude and phase of the actuation relative to the perturbation. In the nonlinear equations, however, if the unstable mode is left to grow ‘naturally’ from noise, its amplitude and phase are not specified in the initial condition, so actuation cannot be fixed a priori. A workaround would be to let the unstable mode develop up to a trigger amplitude, where both amplitude and phase can be measured, and then apply the optimal actuation based on these. Indeed, such an approach would be required in a real process where the onset of instability would need to be detected for the actuation to be activated. For instance, electromagnetic sensors (Thomas et al. Reference Thomas, Yuan, Sivaramakrishnan, Shi, Vanka and Assar2001; Cho & Thomas Reference Cho and Thomas2019) can detect the onset of such instabilities. Since our previous study of the free nonlinear evolution of the leading eigenmode confirmed that it emerged naturally from noise (Kumar & Pothérat Reference Kumar and Pothérat2020), we shall directly initiate the nonlinear simulation with the leading eigenmode set to an amplitude such its ratio to the base flow,

(5.1) \begin{equation} \mathcal{R}_A = \frac {(\textbf{q}^\prime, \textbf{q}^\prime)}{(\textbf{Q},\textbf{Q})}, \end{equation}

remains much smaller than unity (see table 4). Here, $\textbf{Q}=(\textbf{U}, \bar {T})^\top$ represents the steady base flow. Note that for the nonlinear simulations, $\textbf{q}^\prime$ is obtained by subtracting the base flow from the entire flow field.

Table 4. Parameters of our numerical calculations: Reynolds number $Re$ , Rayleigh number $Ra$ , order of polynomial $N$ , time step $\Delta t$ . Here amplitude ratio $\mathcal{R}_A$ , the forcing cutoff time $t_f$ , nonlinear efficiency time $t_{{ eff}}=t_2-t_1$ and actuation efficiency $\eta _a=t_R/t_f$ based on recovery time $t_R$ such that $(\textbf{q}^\prime, \textbf{q}^\prime)_{t=t_R}=(\textbf{q}^\prime, \textbf{q}^\prime)_{t=0}$ for the 3-D DNS.

Figure 8. As a function of time $t$ , the nonlinear evolution of the perturbation energy $(\textbf{q}^\prime, \textbf{q}^\prime)$ , in the simulation cases (a) C1, (b) C2, (c) C3 and (d) C4, without actuation (blue), with actuation maximising $t_f$ (orange) and with actuation achieving maximum reduction of energy (red).

5.2. Analysis of cases C1–C4

The evolution of the perturbation with and without actuation for case C1 is shown in figure 8(a). The blue curve illustrates the free nonlinear evolution with the flow initialised with the leading eigenmode for case C1 and no actuation applied. Initially, the energy of the perturbation $(\textbf{q}^\prime, \textbf{q}^\prime)$ grows exponentially, as predicted by the linear theory, and this confirms that nonlinear effects do not affect the early stages. These indeed induce a saturation at $t=t_1=25.4$ . The envelope of the curve provides an estimate of the saturation time, denoted as $t_1$ . We consider the saturation time to be reached when the derivative of the envelope with respect to time falls below $10^{-3}$ . The red curve represents the evolution of the perturbation for the same initial conditions, but this time, with the suppression-optimising actuation applied as we just described. Until $t \approx 10$ , the actuation results in a decrease of the perturbation’s amplitude, after which the perturbation begins to grow slowly. Interestingly, the decay of the perturbation continues well beyond the time of minimum amplitude found in the linear simulation. Since, however, nonlinear effects are not expected to play a role for such low amplitudes, as suggested by the exponential form of the actuation-free evolution, the difference more likely results from switching off the actuation at $t_f$ in the nonlinear simulations. Given the very low amplitude of the perturbation around $t=t_f$ , the system remains governed by linear dynamics so the reason for the extended period of suppression is that exponential growth from such a low level of perturbation simply extends over a longer period of time. This suggests that applying the actuation beyond the time of maximum growth suppression in fact promotes the growth of the eigenmode. As such, switching off the actuation at $t_f$ is optimal. The underlying suppression strategy, is then to suppress the leading eigenmode down to the lowest possible amplitude so that the linear mechanism takes as long as possible to act. The nonlinear dynamics only provide a measure of when this strategy ceases to be effective.

Indeed, after $t \approx 20$ , the exponential growth becomes sufficient to activate nonlinearities and the perturbation amplitude reaches nonlinear saturation at $t=t_2$ ( $t_2$ is calculated using the same method as $t_1$ by analysing the envelope of the curve). In order to estimate the effectiveness of the actuation at suppressing the instability, we define the nonlinear efficiency time, $t_{{ eff}}=t_2 - t_1$ , which for case C1 is $19.1 \simeq 5.23t_f$ . In other words, the actuation prevents nonlinear saturation during $5.23$ times the duration of its application, i.e. it approximately doubles the time to saturation. Even more interestingly, the unstable mode returns to its initial amplitude at a time $t_R$ , $\eta _a=5.21$ times longer than $t_f$ , so that it could potentially be indefinitely maintained at that level by applying the actuation during $t_f$ every $\eta _at_f$ , which costs $\eta _a$ less energy than applying the actuation continuously. These results are very encouraging since they validate the linear approach, and show that the receptivity forcing obtained from the linear analysis is indeed capable of suppressing the growth of the instability in the full nonlinear system for a long time, not only compared with the time scale of growth of the instability and to the duration of the actuation $t_f$ .

In addition to applying actuation based on the minimum value of $(\textbf{q}^\prime, \textbf{q}^\prime)$ , we also tested the actuation strategy based on the forcing maximising $t_f$ . For case C1, the highest $t_f$ value of 21.8 is achieved for $A=0.08$ and $\phi =\pi$ . The simulation corresponding to these parameters is depicted by the orange curve in figure 8(a). Interestingly, the perturbation energy remains below that of the unforced case until around $t \approx 20$ , yet it consistently stays higher than that of the mode with suppression-optimising actuation ( $A=0.5$ ), before reaching nonlinear saturation. Since the perturbation energy for $A=0.08,\phi =\pi$ exceeds that of $A=0.5, \phi =\pi$ , its nonlinear saturation occurs prior to the case of $A=0.5$ . This means that suppression-optimising actuation outperforms time-optimising actuation.

Case C2 returned exactly the same phenomenology for the suppression-optimising actuation, despite the difference in the nature of the leading eigenmodes (see figure 8 b). The values of $t_{{ eff}}/t_f$ and $\eta _a$ are not significantly different than for case C1, but the time-optimising actuation incurs a different behaviour: although the unstable mode does not initially decay, it decays briefly at a much later time than for the actuation optimising suppression, down to a similar level. Shortly after, however, both curves practically coincide and simultaneously evolve to saturation. In this sense both actuations perform equally well: this is the only case where the suppression-optimising actuation does not significantly outperform the time-optimising actuation. Nevertheless, since the perturbation remains at a much lower level in the initial stages of evolution for the suppression-optimising actuation, it is still overall better at suppressing the unstable mode.

Case C3 also produced the same phenomenology shown on figure 8(c), but with increased values of $t_{{ eff}}/t_f=8.29$ and $\eta _a=8.27$ , i.e. with a much higher efficiency of the suppression-optimising actuation. By contrast with case C2, however, the time-optimising actuation is practically ineffective at suppressing the unstable mode as its amplitude grows to saturation barely later than without any actuation.

Lastly, we applied constant actuations in the nonlinear simulation of case C4, since the unstable mode is non-oscillatory. Once again, we used the unstable mode as the initial condition. The non-oscillatory nature of the mode does not seem to play a role as the evolution of the mode’s energy in all three simulations, without actuation, with suppression-optimised actuation and with time-optimising actuation show very similar features to cases C1 and C3, with $t_{{ eff}}/t_f \approx \eta _a=3.93$ .

5.3. The salient features of the nonlinear evolution

The overall outcome of the nonlinear simulations is that the phenomenology seen in the linear simulation remains valid until the nonlinear effects become important and crucially, this happens when the saturation starts. If an actuation based on the adjoint mode, applied up to the time of maximum suppression $t_f$ the unstable mode does not regain its initial amplitude until a time $\eta _a t_f$ , roughly an order of magnitude longer than $t_f$ ( $\eta _a=8.27$ in the most favourable case C3), so as long as the amplitude and phase of the actuation are optimised for suppression (and not for $t_f$ ). Until that point, the evolution follows mostly the linear dynamics, as the energy of the unstable mode remains small compared with that of the base flow. The nonlinear effects act shortly after this time, and when they do so, they incur growth up to the same point of saturation as when no actuation is applied.

In light of this behaviour, an on–off control strategy could offer a viable control strategy. Repeating a sequence where the suppression-optimising actuation is applied until $t_f$ , then switched off to let the flow evolve until $\eta _a t_f$ , may indefinitely keep the flow evolution in the linear regime, where the actuation remains effective. Thus, this strategy may confine oscillations to very small amplitudes for an arbitrary length of time. Though, this strategy needs to be further explored to assess how feasible this approach is in practical applications.