Modal and non-modal stabilities of ﬂow around a stack of plates

Modal and non-modal stabilities of ﬂow around a stack of ﬂat plates are investigated by means of asymptotic stability and transient growth analyses respectively. It is observed that over the parameters considered, both the base ﬂow and the stabilities vary as a function of ReW 2 / ( W − 1) 2 , i.e. the product of the Reynolds number and the square of the expansion ratio of the stack. The most unstable modes are found to be located downstream of the recirculation bubble while the global optimal initial perturbations (resulting in maximum energy growth over the entire domain) and the weighted optimal initial perturbations (resulting in maximum energy growth in the close downstream region of the stack) concentrate around the stack end owing to the Orr mechanism. In direct numerical simulations (DNS) of the base ﬂow initially perturbed by the modes, it is noticed that the weighted optimal initial perturbation induces periodic vortex shedding downstream of the stack much faster than the most unstable mode. This observation suggests that the widely reported vortex shedding in ﬂow around a stack of plates, e.g. in thermoacoustic devices, is associated with perturbations around the stack end.


Introduction
Flow around a stack of plates has been widely encountered, e.g. in thermoacoustic devices applied in electricity generators, refrigeration or heat recovery systems for conversion between thermal and acoustic energies. Parameter dependence of the vortex shedding in flow past a stack of parallel plates has been extensively studied, e.g. effects of stack positioning, blockage 10 ratio, plate aspect ratio and Reynolds number [7,8,9,10]. To reduce the number of independent parameters, the flow past a stack can be simplified as flow past in confined planar wake flow [11] and flow around arrays of rectangular cylinders 15 [12] or circular cylinders [10], where the critical Reynolds number above which instabilities or vortex shedding occur was calculated. For flow at Reynolds number well above the critical value, three-dimensional stability of the wake flow around a thin plate has been studied [13]. In the rest of this work, the methodology to calculate perturbations relevant to vortex shedding, e.g. modal stability (asymptotic stability) theory to calculate the modal modes (most unstable modes) and nonmodal stability 30 (non-normality, transient growth) theory to calculate the nonmodal modes (optimal initial perturbations) are presented, followed by the numerical setup, and then the modal and nonmodal modes as well as their nonlinear developments to vortex shedding are discussed.

35
Both modal and nonmodal studies performed in this work involve the linearisation of the incompressible Navier-Stokes (NS) equations: where u is the velocity vector, p is the modified pressure and Re is the Reynolds number. The flow field can be decomposed into the summation of a steady base flow and a perturbation flow as (u, p) = (U , P ) + (u , p ).
(2) Substituting eq. (2) into eq. (1) and linearizing the convective term, the linearized NS equations are obtained: In the Cartesian coordinates with x, y and z denoting the streamwise, vertical and spanwise directions respectively as illustrated in fig. 1, the base flow 3 around a stack of plates is assumed to be homogenous in the spanwise direction z. Therefore any perturbation at time t = T can be further decomposed as the sum of modal modes u (x, y, z, T ) = û ij (x, y)e σiT e βj z whereû ij denotes the modal mode with spanwise wavenumber β j , growth rate 40 Re(σ i ) and frequency Im(σ i ). At each given spanwise wavenumber, the most unstable mode can be calculated as the modal mode with the largest growth rate. Clearly if all the growth rates are negative, all the modal modes decay in time and the flow is stable to perturbations, while if at least one growth rate is positive, the flow is asymptotically unstable.

45
To calculate the most unstable mode and its growth rate at a given β, define an operator A, which evolves a perturbation from t = 0 to t = T by integrating the linearized NS equations: The growth rate of the most unstable mode can be obtained as σ max = T −1 lnλ, where λ is the eigenvalue of A. This (dominant) eigenvalue can be calculated by implementing an Arnoldi method to a Krylov sequence built by iterative calls of the linearized NS equations [14]. In a convergence test, it is observed that as expected the growth rate is independent on the value of T . However for larger 50 values of T , the growth rate converges over less iterative calls of the governing equations but each iteration costs more CPU hours. In this work T = 5 is adopted.
While the modal stability analyses focus on asymptotic growth of perturbations over infinite time horizons, dynamics of perturbations over finite time, or transient energy growth, can be investigated by nonmodal analyses. Transient growth is defined with respect to the energy growth of the perturbation over a given time interval [15] and can be quantitatively measured as the maximum ratio of the final perturbation energy and the initial energy across all possible initial perturbations [14]: where τ is a final time and the scale product is defined as (u , u ) ≡ Ω u · u dv, with Ω denoting the computational domain.
To evaluate the energy growth in a particular area of the computational 65 domain, e.g., the region around the exit of the stack where the perturbation has significant impacts on the efficiency of the thermoacoustic system, a weighted energy growth is defined: where F is a non-negative spatial weight function to filter the energy growth in the "uninterested region", e.g. region far downstream of the stack. This 70 weighted maximum transient growth and the corresponding optimal initial perturbation are the largest eigenvalue and the corresponding eigenvector of the operator A * (τ )F 2 A(τ ). In the numerical calculation, F can be considered as a diagonal matrix with entries between 0 and 1. In the following G and G F will be referred to as the global and weighted transient energy growth, respectively.

Numerical setup
In this work, the stack is isolated and a constant inflow condition is implemented on the inflow boundary (see fig. 1). The inflow velocity and the width of the plate are used to define the Reynolds number and therefore the nondimensionalized plate width is D = 1. The length of the plate is fixed at L = 10 80 to generate a thin plate and the domain width W is a free parameter. Therefore the flow expands with a ratio W/(W − 1) when exiting the channel.
In the numerical setup, the centre of the plate is located at the origin, the to be zero and for nonmodal stability studies, both the inflow and outflow velocity conditions are set to Dirichlet zero [14]. When the velocity condition is of Dirichlet type, a computed Neumann pressure condition is adopted [16].

100
In this section, the steady base flow is firstly calculated and then the modal and nonmodal studies are carried out to calculate the most unstable modes and optimal initial perturbations. Finally, the nonlinear developments of the modes are studied to identify their role in activating vortex shedding through DNS of the base flow initially perturbed by the modes.  flow around a stack of plates in thermoacoustic devices [18], as well as in other 115 bluff body flows, e.g. flow around a cylinder cascade [19] and flow past square cylinders [20].

Modal Stability Analyses
The geometry with W = 5.6 matches the experimental setup used in [2] and 135 is therefore adopted as an example to illustrate the instabilities at various spanwise wavenumber and Reynolds number, as presented in fig. 5. It is observed that the boundary of instability, featured by the contour line Re(σ max ) = 0 expands to higher values of spanwise wavenumber as the Reynolds number increases. It is also noted that the two-dimensional mode (β = 0) is always the 140 most unstable one for a given Reynolds number and therefore the following studies focus on two-dimensional instabilities. To quantitatively evaluate effects of Reynolds number and the expansion ratio, the critical Reynolds number, above which the flow becomes unstable, is plotted in fig. 6(b). It is observed that the critical Reynolds number increases linearly with (W −1) 2 /W 2 and Re C W 2 /(W −1) 2 ≈ 0.00907 over the parameters

Nonmodal Stability Analyses
The distributions of all the calculated unstable modes are concentrated far downstream of the stack, and would require a large time interval to reach the stack and activate vortex shedding. To reveal the most effective route from ini- The global transient energy growth G as well as the weighted transient growth G F at W = 5.6 is shown in fig. 8. It is noticed that there is significant transient growth even in the asymptotically stable condition (Re < 74.1 as shown in fig. 6). For the weakly unstable cases, e.g. Re = 100, the transient   fig. 7, the weighted optimal initial perturbations are located around the stack end, and therefore can be convected downstream and activate vortex shedding around the stack faster than the unstable modes. Physically these perturbations could stem from wall roughness around the stack end. As τ increases, the 195 perturbation moves upstream slightly so as to keep the bubble region disturbed over the range 0 ≤ t ≤ τ . Compared with the weighted optimal initial perturbation, the global optimal initial perturbation is located more downstream and therefore relies less on the amplification in the bubble region.
It is worth noting that the structures of the optimal initial perturbations 200 tilt backwards so as to take advantage of the Orr mechanism when they are fig. 10(a) and observed in several other nonmodal stability analyses [21,14].
For further developments of the perturbation until t = 50 (see fig. 10  For all the three points considered, it is seen that the optimal initial perturbation leads to periodic oscillations of the vertical velocity much faster than the most unstable mode. As the most unstable mode travels upstream, it perturbs the region downstream of the bubble first ( fig. 11b, c) and takes over 100 time units to reach the trailing edge region of the stack ( fig. 11a). In contrast, the  It is observed that the base flow varies as a function of ReW 2 /(W − 1) 2 , which is the product of the Reynolds number and the square of the expansion ratio of the stack flow, over the parameters considered.

250
Then the steady base flow is used in modal and nonmodal stability analyses to reveal the asymptotic and transient energy growth of perturbations, respectively. In modal stability studies, the growth rate of the most unstable mode is found to increase at increasing Reynolds number or the stack expansion ratio, indicating a destabilisation effect of the flow expansion or confinement. The crit-255 ical Reynolds number Re C is observed to be a linear function of (W − 1) 2 /W 2 and satisfies Re C W 2 /(W − 1) 2 ≈ 0.00907. Therefore the two free parameters, i.e. the Reynolds number and expansion ratio, can be merged into one parameter, i.e. ReW 2 /(W − 1) 2 , in studies of both base flow and instabilities. The unstable modes, in the form of odd-odd modes as discussed in [10], are located 260 in the region far downstream of the plates and are associated with the shear layers downstream of the recirculation bubble.
In nonmodal stability analyses, it is seen that both the global optimal initial perturbations (resulting in maximum energy growth across all the domain) and the weighted optimal initial perturbation (resulting in maximum energy growth The route from initial perturbations to vortex shedding is investigated through 270 DNS of the base flow initially perturbed by the modes calculated in stability analyses. It is observed that the weighted optimal initial perturbations are convected downstream and trigger vortex shedding much faster than the most unstable modes, which travel upstream to activate vortex shedding. Overall the linear and nonlinear dynamics of the perturbations investigated in this work 275 indicate that even when the steady flow is asymptotically unstable, the most unstable modes are ineffective in activating the bifurcation to vortex shedding since the upstream noise can be much faster in perturbing the flow.