Controlled mixing enhancement in turbulent rectangular jets responding to periodically forced inflow conditions

On: 13 April 2015, At: 13:20 Publisher: Taylor & Francis

Informa Ltd Registered in England and Wales Registered Number: 1072954 Registered office: Mortimer House, 37-41 Mortimer Street, London W1T 3JH, UK

Click for updates

Journal of Turbulence

Publication details, including instructions for authors and subscription information:

http://www.tandfonline.com/loi/tjot20

Controlled mixing enhancement in

turbulent rectangular jets responding

to periodically forced inflow conditions

a

Department of Mechanical Engineering and Computer Science, Czestochowa University of Technology, Czestochowa, Poland b

Multiscale Modelling and Simulation, University of Twente, Enschede, The Netherlands

c

Anisotropic Turbulence, Faculty for Applied Physics, Eindhoven University of Technology, Eindhoven, The Netherlands

Published online: 08 Apr 2015.

To cite this article: Artur Tyliszczak & Bernard J. Geurts (2015) Controlled mixing enhancement

in turbulent rectangular jets responding to periodically forced inflow conditions, Journal of Turbulence, 16:8, 742-771, DOI: 10.1080/14685248.2015.1027345

To link to this article: http://dx.doi.org/10.1080/14685248.2015.1027345

PLEASE SCROLL DOWN FOR ARTICLE

Taylor & Francis makes every effort to ensure the accuracy of all the information (the “Content”) contained in the publications on our platform. However, Taylor & Francis, our agents, and our licensors make no representations or warranties whatsoever as to the accuracy, completeness, or suitability for any purpose of the Content. Any opinions and views expressed in this publication are the opinions and views of the authors, and are not the views of or endorsed by Taylor & Francis. The accuracy of the Content should not be relied upon and should be independently verified with primary sources of information. Taylor and Francis shall not be liable for any losses, actions, claims, proceedings, demands, costs, expenses, damages, and other liabilities whatsoever or howsoever caused arising directly or indirectly in connection with, in relation to or arising out of the use of the Content.

This article may be used for research, teaching, and private study purposes. Any substantial or systematic reproduction, redistribution, reselling, loan, sub-licensing, systematic supply, or distribution in any form to anyone is expressly forbidden. Terms &

and-conditions

Controlled mixing enhancement in turbulent rectangular jets

responding to periodically forced inflow conditions

Artur Tyliszczaka∗_{and Bernard J. Geurts}b,c

a_{Department of Mechanical Engineering and Computer Science, Czestochowa University of}

Technology, Czestochowa, Poland;b_{Multiscale Modelling and Simulation, University of Twente,}

Enschede, The Netherlands;c_{Anisotropic Turbulence, Faculty for Applied Physics, Eindhoven}

University of Technology, Eindhoven, The Netherlands

(Received 8 October 2014; accepted 5 March 2015)

We present numerical studies of active flow control applied to jet flow. We focus on rectangular jets, which are more unstable than their circular counterparts. The higher level of instability is expressed mainly by an increased intensity of mixing of the main flow with its surroundings. We analyse jets with aspect ratio Ar= 1, Ar= 2 and Ar= 3 at Re= 10,000. It is shown that the application of control with a suitable excitation (forcing) at the jet nozzle can amplify the mixing and qualitatively alter the character of the flow. This can result in an increased spreading rate of the jet or even splitting into nearly separate streams. The excitations studied are obtained from a superposition of axial and flapping forcing terms. We consider the effect of varying parameters such as the frequency of the excitations and phase shift between forcing components. The amplitude of the forcing is 10% of the inlet centreline jet velocity and the forcing frequencies correspond to Strouhal numbers in a range St= 0.3–0.7. It is shown that qualitatively different flow regimes and a rich variety of possible flow behaviours can be achieved simply by changing aspect ratio and forcing parameters. The numerical results are obtained applying large eddy simulation in combination with a high-order compact difference code for incompressible flows. The solutions are validated based on experimental data from literature for non-excited jets for Ar= 1 at Re = 1.84 × 105 and Ar= 2 at Re = 1.28 × 105. Both the mean velocities as well as their fluctuations are predicted with good accuracy.

Keywords: active and passive controls; mixing enhancement; bifurcation phenomenon; LES

1. Introduction

Interest in flow control techniques is driven by the potential to gain considerable improve-ment of performance, safety and efficiency of various technical devices with limited energy input. Flow control may be divided into two categories: passive and active.[1,2] The former most often relies on optimisation procedures which are based on geometric shaping or adding fixed elements (obstacles, swirlers, etc.) to the flow domain. Active methods require the input of energy to the flow whose type and level may be constant (predetermined control methods) or varying in response to the instantaneous flow behaviour (interactive methods). An evident advantage of passive flow control is its low cost and simplicity. However, modifications to the flow domain are usually not optimal simultaneously for different flow

∗_{Corresponding author. Email: atyl@imc.pcz.czest.pl}

C

2015 Taylor & Francis

conditions. From this point of view, active flow control, while being more intricate to realise, is much more flexible and could result in good overall response under a variety of different flow conditions.

Fundamental research is conducted presently into both active and passive flow control techniques. A prominent example of successful alteration of flow dynamics by flow control is found for the class of jet flows. Research on passive control techniques concentrates on geometrical modifications of the jet nozzle. It turns out that non-symmetric jets emanating from rectangular or elliptical nozzles enhance mixing between the jet and the surrounding flow. Particularly, the large-scale mixing of rectangular jets was found substantially larger than in the case of classical circular jets. This aspect is important in many applications (e.g., combustion processes), and therefore, rectangular jets were extensively studied for many years, both experimentally [3–5] and numerically using large eddy simulation (LES) [6–9] and direct numerical simulation (DNS).[10,11] A comprehensive discussion on theoretical issues related to rectangular jets and a list of their possible applications may be found in review papers.[12,13]

Concerning active flow control methods applied to jet type flows, most of the research is devoted to circular jets. The work of Crow and Champagne [14] was probably the first in this category. It was reported that for properly chosen forcing (excitation), the jet behaviour may change qualitatively. An enhanced mixing and the existence of two maxima in the turbulence intensity as function of forcing frequency were found, not seen previously in natural jets. This seminal work initiated many experimental studies which revealed the large potential of active control techniques.[15–23]

In this paper, we combine passive and active flow control techniques to rectangular turbulent jets, and analyse the resulting flow field. The former is obtained by using a rectangular shape with various aspect ratios of the jet and the latter is obtained by imposing at the jet inlet an excitation with a specific amplitude and frequency. In [18], it was shown that for plane jets the influence of a symmetric excitation on the mean and fluctuating velocity fields is much weaker than that in the circular jet. Recently, DNS studies [24] of a planar jet with a steady modulation characterised in terms of a Beltrami flow have been performed. It was shown that the dynamics and size of the mixing layers can be considerably changed by varying the wave number of the applied modulation. The large-scale mixing could be stimulated when the length scale of the imposed modulation had a size comparable to the width of the jet. In the present work, unlike in [18,24], the excitation is a combination of the symmetric and flapping forcing modes and has an unsteady character. We focus on excited jets with aspect ratio Ar = H/W = 1, Ar = 2 and Ar = 3, where H and W

are the inlet dimensions along its major and minor axes, respectively. Emphasis is put on determining inflow perturbation procedures that give rise to a global alteration of the flow pattern and modification or even creation of large-scale flow structures. It is found that such a combined excitation significantly changes the spreading rate of evolving jets, amplify the vortices and for some forcing parameters even initiates a division of the main jet stream into separate streams. This last phenomenon is called the bifurcation phenomenon [25] and is a spectacular example of a modified flow pattern arising from an upstream forcing. The existence of a bifurcation was revealed also for square [11] and rectangular jets.[9] However, they were found to bifurcate with greater difficulty as compared to circular jets. This was attributed to the presence of rib vortices at the corners. In the present work, we extend these studies and quantify the effect of the frequency of the excitation on the bifurcation phenomenon.

The paper is organised as follows. The next section gives details of the governing equations, LES modelling and applied numerical method. In Section 3, the simulation

parameters and forcing method are described. The results of simulation are presented in Sections4and5. First, the validation of solutions is performed based on comparison with literature data for non-excited jets at aspect ratio Ar= 1 and Ar= 2. Next, the simulation

results for the excited jets are presented and discussed, emphasising (1) effect of the forcing frequencies of the axial and flapping excitations, (2) phase shift between these forcing modes and (3) differences between the excitation effects for varying Ar. Finally, concluding

remarks are given in Section6.

2. Governing equations and numerical algorithm

In this section, we introduce the LES model used for the simulations and summarise the main characteristics of the numerical method. The LES approach has been developing for more than 30 years now, and is regarded as a reliable tool in computational fluid dynamics (CFD) simulations.[26,27] In this period, considerable effort has been put on many different issues, including filtering techniques, commutation errors, sub-filter modelling and also mutual interactions between numerically induced errors and modelling errors (for instance,[28– 33]) All this resulted in the current maturity of the LES approach.

In this paper, we consider incompressible flow described by the continuity and the Navier–Stokes equations. In the framework of LES, we have for a commuting filter:

∂uj ∂xj = 0 (1) ∂ui ∂t + ∂uiuj ∂xj = − 1 ρ ∂ ¯p ∂xi + ∂τij ∂xj + ∂τ_ijf ∂xj (2)

where uiare the velocity components, p the pressure andρ the density. The overbar denotes

spatially filtered variables: ¯f = G ∗ f with G being the filter function.[26,27] The stress tensor of the resolved field,τij, and the unresolved sub-filter stress tensorτ_ijf due to filtering

of the non-linear advection terms are

τij = 2νSij, τijf =



uiuj − uiuj (3)

whereν is the kinematic viscosity and S_ij = 1₂  ∂ ¯ui ∂xj + ∂ ¯uj ∂xi 

is the rate of strain tensor of the resolved velocity field. The sub-filter tensor is modelled in this paper by an eddy-viscosity model with τ_ijf = 2ν_tS_ij+ τ_kkfδ_ij/3. The diagonal terms τ_kkf are added to the pressure resulting in the so-called modified pressure ¯P = ¯p − ρτ_kkfδ_ij/3.[26] Hence, we have

∂ui ∂t + ∂uiuj ∂xj = − 1 ρ ∂ ¯P ∂xi + ∂ ∂xj  (ν + νt) _{∂ ¯u} i ∂xj + ∂ ¯uj ∂xi  (4)

In the present implementation, ¯P is calculated with the help of the projection method. The eddy viscosityνtis computed using the model proposed by Vreman [34]:

νt = C



Bβ

αijαij (5)

αij =∂u_∂xj i, βkl = 2_α mkαml (6) Bβ= β11β22− β122 + β11β33− β 2 13+ β22β33− β232 (7)

where the constant in Equation (5) is taken as C= 2.5 × 10−2.[34] The filter width is

computed as = ( x y z)1/3 _{with x, y, z being the mesh spacings. This model}

yields a vanishing eddy viscosity close to a solid wall, in laminar flows or in pure shear regions. This is a very important property of a sub-filter model for jet flows where turbulence develops in the shear layer region. Other sub-filter models, such as Smagorinsky’s model, could be too dissipative and prevent turbulent fluctuations to develop properly in the flow.

The flow solver used in this work is an academic high-order code (SAILOR) which is based on a projection method with time integration performed by a predictor–corrector (Adams-Bashforth, Adams-Moulton) method. The eddy viscosity is computed every time step at the beginning of the projection step. The spatial discretisation is based on a compact difference method developed for half-staggered meshes,[35,36] where the pressure nodes are shifted half a cell size from the velocity nodes. The solution algorithm requires both interpolation procedures and derivative approximations which are performed using 10th-order and 6th-10th-order formulas,[36] respectively.

The SAILOR code was used previously in various studies, including laminar/turbulent transition in free jet flows,[37–39] near-wall flows,[40] multi-phase flows [41] and flames.[42,43] Recently, the SAILOR code has been used in parametric studies of ex-cited circular jets.[44] It was shown that it enables accurate prediction of the bifurcation phenomenon in circular jets. Grid-refinement studies showed that the applied high-order discretisation schemes yield grid-independent solutions already at relatively coarse meshes. This allowed for a comprehensive analysis of various forcing parameters and identification of the most important ones from the point of view of large-scale alteration of jet flows. This knowledge is partially discussed later in relation to the dynamics of rectangular jets.

3. Simulation set-up and forcing procedure

In this section, we first introduce the flow domain and subsequently address the procedure with which the flow is forced at the jet inflow.

A schematic view of the nozzle geometry is shown inFigure 1. In the computations, we do not explicitly resolve the flow inside the nozzle. The computational domain is a simple rectangular box with inflow conditions defined by the mean velocity components and their fluctuations.

InFigure 1, the symbols zhand xhdenote the widths in the lateral (x) and spanwise (z)

directions at which the mean velocity is half the value in the centre, Ucl. These distances are

introduced to enable comparisons with literature data that are presented in non-dimensional form using zhand xhas the reference values.

The sizes of the computational domain are 12W in x-direction, 16W in y-direction and 12W in z-direction, with W the width of the minor side of the jet nozzle. To check the level of influence of the side boundaries, additional test computations were performed using a wider computational domain with 16W × 16W along the x–z directions. Effects of the boundaries on the dynamics of the jets were found to be negligible, as will be quantified momentarily.

Figure 1. Schematic view of the rectangular nozzle.

The inlet boundary conditions are specified in terms of the mean velocity profile superimposed with fluctuating components as

u(x, t) = umean(x)+ uexcit(x, t) + uturb(x, t) (8)

where the mean velocity is a hyperbolic-tangent profile:

umean(x)= U_c+Uj 4 1− tanh  1 4 Rx δθ  |x| Rx − Rx |x|  × 1− tanh  1 4 Rz δθ  |z| Rz − Rz |z|  (9)

in which Ujdenotes the inlet velocity at the jet axis and Uc= 0.05Ujis the velocity of a

co-flow introduced in order to mimic an inflow entrainment observed when a jet issues from a nozzle into an open domain.[45] Without adding the co-flow, a free jet would become a jet in an enclosure with a recirculation zone created close to the inlet. Adding the co-flowing stream is a standard procedure in jet flow simulations.[23,46,47] At the level Uc< 0.1Uj,

the co-flow has only a minor influence on the jet dynamics.[47] The symbols x, z in Equation (9) denote the in-plane coordinates, Rx= W/2 and Rz= H/2 are the nozzle half-width of

the minor and major axes. The parameterδ_θis the momentum thickness of the initial shear layer. In all simulations performed in this work, we takeδ_θ= 0.05Rxwhich is the same as

in [44,47] for natural and excited circular jets. This choice is also motivated by the fact that for thinner shear layers, i.e.,δ_θ ≤ 0.05Rx, and turbulence intensity of the order of 1% the

jets exhibit very similar behaviour.[39]

The velocity component uexcit(x, t) acts as the deterministic forcing term which in physical experiments is usually produced by a membrane (loudspeaker) located upstream of the nozzle or by plasma actuators or by a mechanical forcing obtained by especially designed flap actuators placed at the nozzle lip.[17,48,49] These excitations change the direction and magnitude of the flow leaving the nozzle. In this work, we add the forcing

Figure 2. Schematic view of the forcing terms: axial mode is shown on the left side of the figure and the flapping mode is displayed on the right-hand side.

only to the streamwise velocity as

uexcit(x, t) = Forcing × Mask (10)

where

Forcing= Aasin2πStaUjt/W+ Afsin2πStfUjt/W + sinπx

W

 (11) is the superposition of the axial forcing mode with amplitude Aaand the flapping forcing

mode with amplitude Af. The mask is defined as

Mask= 1 4 1− tanh  1 4 Rx δθ _|x| Rx − Rx |x|  1− tanh  1 4 Rz δθ _|z| Rz − Rz |z|  (12) and provides a smooth transition with the co-flowing stream. The applied forcing locally alters the magnitude of the inlet velocity but the direction of the flow remains unchanged. This method was successful in previous studies [44,47,50] for the circular jets. The Strouhal numbers of the excitations are defined as Sta= faW/Ujand Stf= ffW/Uj, where fa, ffare the

frequencies of the axial and flapping excitations. The symbol is the phase shift between axial and flapping forcing. In this paper, we consider the flapping excitation along the lateral direction only, as shown schematically inFigure 2. Additional studies performed with the excitation applied along the spanwise direction showed that for the cases with Ar≥ 2, the

only effects of such forcing were an increase of the turbulence intensity in the flow field, but not a global change of the flow pattern.

Finally, the turbulent fluctuationsuturb(x, t) are computed by applying a digital filtering method to random inflow perturbations as proposed by Klein et al. [51]. This method guarantees spatially correlated velocity fields which can be tuned to reflect real turbulent flow conditions.[51–55]

At the side boundaries, the streamwise velocity is assumed equal to Ucand the remaining

velocity components are equal to zero. Hence, there is no flow through the side boundaries. The pressure is computed from the Neumann condition n· ∇p = 0 with n the outward normal vector. At the outlet, all velocity components are computed from a convective boundary condition∂ui/∂t + VC∂ui/∂n = 0. The convection velocity VCis computed every

time step as the mean velocity in the outlet plane, limited such that VC = max (VC, 0).

Table 1. Parameters of the computational cases for non-excited jets, specified by aspect ratio Ar, Reynolds number Re and the computational grids. The cell sizes are expressed as x∗ = x/W,

y∗_{= y/W, z}∗_{= z/W with x, y and z taken from the non-uniform grid at x = y = z = 0.}

Case Ar Re Nx × Ny × Nz ( x∗, y∗, z∗) A1-1 1 1× 104 _{192 × 256 × 192 (0.034, 0.042, 0.034)} A1-2 1 1.84× 105 _{192 × 256 × 192 (0.034, 0.042, 0.034)} A1-3 1 1.84× 105 _{320 × 384 × 320 (0.020, 0.028, 0.020)} A2-1 2 1× 104 _{192 × 256 × 192 (0.034, 0.042, 0.054)} A2-2 2 1.28× 105 _{192 × 256 × 192 (0.034, 0.042, 0.054)} A2-3 2 1.28× 105 _{320 × 384 × 320 (0.020, 0.028, 0.032)} A3-1 3 1× 104 _{192 × 256 × 192 (0.034, 0.042, 0.058)}

The term ∂ui/∂n is discretised with a second-order tri-points upwind scheme.[56] The

pressure at the outflow is assumed constant and equal to zero. This formulation of the outflow boundary conditions is stable and allow to pass turbulent flow structures without any visible distortion.[42–44]

4. Natural, non-excited jets – validation test cases

The first set of simulations is performed for non-excited jets. These solutions will be used as point of reference for comparison with the excited jets analysed in the next section. Additionally, the results for the jets with aspect ratios Ar = 1 and Ar = 2 illustrate the

accuracy of the applied numerical method and assess the dependence of the solution on inlet parameters and on the density of the computational grid.

We refer to experimental data [3,4] obtained for Reynolds numbers Re= 1.84 × 105

for Ar= 1 and Re = 1.28 × 105for Ar= 2. The Reynolds number for rectangular jets

is defined as Re= UjDe/ν, where D_e= 2

√

WH/π is the so-called equivalent diameter.

In order to check the dependence of the solutions on Re and also to enable comparison between the solutions for different Ar, we additionally consider flow at Re= 1 × 104.

The computational domain is 12W× 16W × 12W. The numerical grids are compacted towards the jet centre and along the axial direction. In the x- and z-directions, the mesh nodes are distributed using a tangent hyperbolic function with parameters chosen such that almost uniform cell sizes are generated across the jet inlet. In the axial direction, the grid nodes are distributed using an exponential function. Two computational meshes are used. The first one consists of 192 × 256 × 192 nodes and will be further referred to as the basic mesh. We also include a refined mesh with 320 × 384 × 320 nodes to assess the level of dependence of the results on the mesh density.Table 1reports the analysed test cases and gives details of the applied numerical grids. The cell sizes x, z and y correspond to the location x= z = y = 0. The symbols ‘A1’, ‘A2’ and ‘A3’ indicate the aspect ratio and the different flow cases are numbered 1, 2 and 3. Ratios of the maximum (near the boundaries) to minimum (in the inlet jet region) cell sizes are (1) 6.41 and 6.47 in the x-direction, (2) 2.081 and 2.083 in the y-direction, (3) 2.06 and 2.07 in the z-direction for Ar= 2 and 1.56

for Ar= 3, where the first number refers to the coarse mesh and the second one to the ratio

on the dense mesh.

LES computations for the reference cases analysed in [3,4] were carried out previously in [7] for three different velocity profiles. Zero and a non-zero spanwise and lateral velocity components were analysed and compared with mean values taken from experiments. It was

found that the best agreement with the reference data was obtained when the streamwise velocity profile was assumed as a sharp rectangular profile and the remaining velocity components were set to zero. Here, we also adopt this approach and take the mean spanwise and lateral velocities zero, whereas the streamwise velocity profile is approximated by Equation (9). This profile closely resembles both the experimental profile and the profile used in [7]. The mean profile (9) with δ_θ = 0.05Rxyields a shear layer thicknessδ99 =

0.48Rxwhich is defined here as the region whereU_c< umean(x)≤ 0.99U_j. In the lateral

direction, the number of grid nodes located acrossδ99is equal to 7 and 12 on the basic and dense mesh, respectively. In the spanwise direction for Ar= 2, there are, respectively, five

and eight nodes depending on the mesh density, and for Ar= 3, there are four nodes on the

basic mesh.

Concerning the distribution of turbulence intensities and turbulent length/time scales, no detailed data were documented in the experiments. Only the minimum value of the turbulence intensity on the jet axis was reported equal to Ti= 0.005Uj. In LES reported in

[7], the turbulent fluctuations were neglected and the inflow was treated as laminar arguing that a precise reproduction of the time-dependent inflow conditions is very difficult. In the present simulations, we generate approximate inflow perturbations such that the fluctuating velocity components are calculated from an assumed magnitude of Ti and assumed turbulent length and time scales ll, lt needed to generate the turbulent inflow conditions using the

digital filtering method. Such fluctuations are then added to the mean velocity profile computed from Equation (9). For simplicity, we assume a uniform turbulence intensity across the nozzle at Ti = 0.01Uj and will investigate the influence of variations in ll, lt

later.

In the simulations, initially, the velocity field is set to Uceverywhere in the computational

domain except the inflow plane where the velocity is computed from Equation (8). During an initial transient phase of about 100T0 (T0 = W/Uj), the jet develops in the domain.

Subsequently, the solution is averaged over 800T0which was found to be sufficient to obtain

steady statistics. The convergence of the time averaging of the solutions was monitored by comparing the averaged velocity profiles between 600T0 and 800T0, observing less than

3% maximal deviation.

4.1. Influence of the turbulent scales ll, lt

In several simulations, we analyse the influence of ll, lton results obtained for the cases A1-2

and A2-2. The computations are performed with ll= 0.2W, 0.1W, 0.05W and lt= 0.5W/Uj, W/Uj, 2W/Uj, 5W/Uj. The differences between the results obtained are readily apparent.

Figure 3shows the profiles of the time-averaged streamwise velocity and its fluctuations

along the centreline. The comparison is presented for four cases with ll= 0.1W, ll= 0.2W

and lt= W/Uj, lt= 2W/Uj. Results with fluctuations generated using random noise are also

included inFigure 3, denoted as ll= lt= 0. Typical for the randomly generated perturbations

is their rapid dissipation immediately downstream of the inlet (for instance, see Figure 1 in [52]). Indeed, in the zooms showing the profiles of fluctuations close to the inlet, it is seen that the fluctuations corresponding to ll = lt = 0 first decrease to zero and recover only

after some distance. Contrary, results obtained with the turbulence generated using digital filtering show that the initial level of Ti is approximately preserved. Further downstream, the results are dependent on ll, ltand the best agreement with the experimental data seems

to be obtained for the case with ll= 0.2W and lt= 2W/Uj. We adopt these parameters in

the sequel of this paper.

y/D_e Ucl /U ma x (u u ) 1/ 2 /U cl 0 2 4 6 8 10 0 0.2 0.4 0.6 0.8 1 1.2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 ll=0W; lt=0W/Uj ll=0.1W; lt=1W/Uj ll=0.1W; lt=2W/Uj ll=0.2W; lt=1W/Uj ll=0.2W; lt=2W/Uj

Quinn & Militzer(1988)

0 0.5 1

0 0.01 0.02

(a) Case: A1-2; Ar= 1, Re=1.84 × 105

y/D_e Ucl /U ma x (u u ) 1/ 2 /U cl 0 2 4 6 8 10 0 0.2 0.4 0.6 0.8 1 1.2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 ll=0W; lt=0W/Uj ll=0.1W; lt=1W/Uj ll=0.1W; lt=2W/Uj ll=0.2W; lt=1W/Uj ll=0.2W; lt=2W/Uj Quinn (1992) 0 0.5 1 0 0.01 0.02

(b) Case: A2-2; Ar= 2, Re=1.28 × 105

Figure 3. Mean stremwise velocity profile and its fluctuating component for the cases A1-2 and A2-2 for different turbulent length and time scales used to generate turbulent inflow conditions based on the digital filtering method.[51] The mean centreline velocity is normalised using the maximum velocity in the flow domain, Umax. The fluctuations are normalised by the centreline velocity, Ucl.

Figure 4. Solution for the non-excited jet with aspect ratio Ar= 1 for the case A1-1. Left figure: isosurface of the vorticity magnitude (| | = 3(Uj/W)). Right figure: contours of the time-averaged axial velocity in the x–y plane.

4.2. Jet with Ar= 1

The instantaneous and the time-averaged results obtained for the non-excited jet with Ar=

1 for the case A1-1 are presented inFigure 4. The instantaneous solution is visualised by the isosurface of the vorticity magnitude equal to | | = 3(Uj/W). A highly turbulent flow

behaviour emerges at the end of the potential core that extends to≈4De(De= 1.128W for Ar= 1). Here, the end of the potential core is taken as the distance from the inlet, where

the time-averaged axial velocity at the centreline drops to 98% of Uj.Figure 4shows also

the contours of the time-averaged axial velocity in the lateral cross-section plane ‘x–y’. Note that because of the symmetry of the flow for Ar= 1, the results in the streamwise

plane ‘z–y’ are the same which confirms that the averaging time is sufficiently long. The profiles of the time-averaged axial velocity obtained from the simulation A1-1, A1-2 and A1-3 are shown inFigures 5and6. These figures present the mean values and fluctuating component along the lateral direction. The solutions are compared with experimental data [3] at the axial locations y/De= 0.28, 2.658, 4.484, 7.088, where the measurements were

performed. The results are normalised by the centreline velocity Ucland the half-velocity

width xhcomputed at the respective y/Delocations.

In Figure 5(a), it is seen that the experimental data exhibit behaviour typical for jets

issuing from sharp-edged slots or from short converging nozzles, as the mean axial velocity at the shear layer is higher than the centreline value. As shown in [3], the overshoot at the nozzle edges was about 20%. This is not taken into account in the present simulations, and therefore in the flow region adjacent to the inflow plane, some discrepancies may be expected. Indeed, concerning the mean velocity profiles, small differences between the experimental and numerical results are observed for the solutions at y/De= 0.28 and y/De=

2.658 as shown inFigure 5(a) and 5(b). Further downstream at y/De= 4.484 and y/De=

7.088, the numerical results match the experimental data very accurately. It is apparent that the influence of the Reynolds number is very small and the profiles for different Re coincide to a high degree.

InFigure 6, the velocity fluctuations are seen to be very close to each other regardless of

Re. Additionally, it is worth noting that the velocity fluctuations of solution A1-3 obtained on the dense mesh are virtually the same as the fluctuations of the solution found on the basic mesh. This means that the mesh density 192 × 256 × 192 is sufficient to provide high-fidelity results for mean and fluctuating quantities. The profiles of the velocity fluctuations agree reasonably well with the experimental data at y/De= 0.28 (Figure 6(a)) and y/De=

x/x_h U/ Ucl -3 -2 -1 0 1 2 3 0.0 0.2 0.4 0.6 0.8 1.0 1.2 (a) y/De= 0.28 x/x_h U/ Ucl -3 -2 -1 0 1 2 3 0.0 0.2 0.4 0.6 0.8 1.0 1.2 (b) y/De= 2.658 x/x_h U/ Ucl -3 -2 -1 0 1 2 3 0.0 0.2 0.4 0.6 0.8 1.0 1.2 (c) y/De= 4.484 x/x_h U/ Ucl -3 -2 -1 0 1 2 3 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Quinn (’88) A1-1 A1-2 A1-3 (d) y/De= 7.088

Figure 5. Profiles of the mean axial velocity along the lateral direction for the jet with Ar= 1 for the cases A1-1, A1-2 and A1-3.

7.088 (Figure 6(d)), whereas in the transient region where turbulent flow develops, i.e., at y/De= 2.658 and y/De= 4.484 (Figure 6(b) and 6(c)), differences are clearly visible.

At these locations, the shapes of the profiles reflect the experimental data correctly, but the maxima of the fluctuations significantly overpredict the measured values, particularly at x/xh = ±1. This may be caused by the fact that in the simulation the level of Ti and

the turbulent scales were taken uniform across the jet radius which may lead to a larger perturbation growth downstream.

4.3. Jet with Ar= 2

An isosurface of instantaneous vorticity magnitude and the contours of the time-averaged axial velocity obtained for the non-excited jet with Ar= 2 for the case A2-1 are presented in

Figure 7. The results are presented both in the lateral and the spanwise planes. The potential

x/x_h (uu) 1/ 2 /U cl -3 -2 -1 0 1 2 3 0.00 0.05 0.10 0.15 0.20 0.25 0.30 Quinn (’88) A1-1 A1-2 A1-3 (a) y/De= 0.28 x/x_h (uu) 1/ 2 /U cl -3 -2 -1 0 1 2 3 0.00 0.05 0.10 0.15 0.20 0.25 0.30 _{Quinn (’88)} A1-1 A1-2 A1-3 (b) y/De= 2.658 x/x_h (uu) 1/ 2 /U cl -3 -2 -1 0 1 2 3 0.00 0.05 0.10 0.15 0.20 0.25 0.30 Quinn (’88) A1-1 A1-2 A1-3 (c) y/De= 4.484 x/x_h (uu) 1/ 2 /U cl -3 -2 -1 0 1 2 3 0.00 0.05 0.10 0.15 0.20 0.25 0.30 Quinn (’88) A1-1 A1-2 A1-3 (d) y/De= 7.088

Figure 6. Profiles of the axial velocity fluctuations along the lateral direction for the jet with Ar= 1 for the cases A1-1, A1-2 and A1-3.

core is shorter than for the jet with Ar= 1 (seeFigure 4) and extends only to 2.2De(De=

1.595W for Ar= 2 ), which in terms of Deis about 40% less than found at Ar= 1. Analysis

of the time-averaged results allows to find a location of the so-called crossover point, i.e., the distance from the inlet where the jet dimensions along the lateral and spanwise axes become equal. In the present case, this is found at y/De≈ 4.25 which fits within the range of

crossover points, i.e., 3< y/De< 7 as put forward in [12]. The profiles of the axial velocity

obtained for the cases A2-1, A2-2 and A2-3 are compared with the experimental data from

[4] inFigure 8. The solutions are presented along the lateral and spanwise directions at the

axial distances y/De= 2, y/De= 5 and y/De= 10 as selected in [4].

The numerical results for the lateral direction agree very well with the experimental data, as seen inFigure 8. The solutions are practically independent of the Reynolds number and the mesh density. The same applies to the results in the spanwise direction. Small discrepancies are found only at y/De= 2 in the shear layer region, where the numerical

Figure 7. Solution for the non-excited jet with aspect ratio Ar= 2 for the case A2-1. Left figures: isosurface of the vorticity magnitude (| | = 3(Uj/W)) in the x–y (upper figures) and z–y plane (bottom figures). Right figures: contours of the time-averaged axial velocity in the x–y and z–y plane.

results show a slightly slower velocity decay. We note that a similar effect was observed in lattice Boltzmann LES reported in [7], where the same experimental data were used as point of reference. Hence, the faster velocity decay in the shear layer region in the spanwise direction may be characteristic for the experimental set-up used in [4]. Finally, the comparison of the turbulent kinetic energy (TKE) defined as q2= 0.5( < u2> + < v2> + < w2_{> ) is shown in}

Figure 9. The shape of the TKE distribution is captured correctly,

although the maxima of the TKE differ from the experiment by 5%–30%, depending on the simulation (A2-1,2,3) and the y/Delocation. At the distance y/De= 10, not only the level

of TKE is underestimated, but also the profiles are flattened in the central part of the jet (−1 ≤ x/xh≤ 1; −1 ≤ z/zh≤ 1). This may well be related to the proximity of the outflow

boundary. Simulations with a slightly longer domain of 12.5Dealso displayed the tendency

towards flattened profiles.

Figure 10shows amplitude spectra of the axial velocity for the case A2-1 calculated

at the location of the shear layer x/W= 0.5. The Strouhal number used on the horizontal axis inFigure 10is defined as StW = fW/Uj, where f is the frequency of the harmonics

of the velocity time signal. The presented spectra were calculated from time series taken from mesh points located at the axial distances y/De= 0.0, 0.5, 1.0, 3.0. FromFigure 10,

it can be seen that there is no distinct dominant frequency which could be interpreted as the preferred mode frequency, instead a broad region of more intense fluctuations is observed. The results for the jets with Ar= 1 and Ar= 3 also do not show the existence of

a well-defined peak at a specific isolated frequency, as usually seen in circular jets.[57,58]

x/x_h U/ Ucl -3 -2 -1 0 1 2 3 0.0 0.2 0.4 0.6 0.8 1.0 1.2 y/De=2 x/x_h U/ Ucl -3 -2 -1 0 1 2 3 0.0 0.2 0.4 0.6 0.8 1.0 1.2 y/D_e=5 x/x_h U/ Ucl -3 -2 -1 0 1 2 3 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Quinn (’92) A2-1 A2-2 A2-3 y/D_e=10 z/z_h U/ Ucl -3 -2 -1 0 1 2 3 0.0 0.2 0.4 0.6 0.8 1.0 1.2 y/De=2 z/z_h U/ Ucl -3 -2 -1 0 1 2 3 0.0 0.2 0.4 0.6 0.8 1.0 1.2 y/De=5 z/z_h U/ Ucl -3 -2 -1 0 1 2 3 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Quinn (’92) A2-1 A2-2 A2-3 y/D_e=10

Figure 8. Profiles of the mean axial velocity along the lateral direction (upper figures) and the spanwise direction at y/De= 2, 5, 10 for the jet with Ar= 2, cases A2-1, A2-2 and A2-3.

x/x_h q 2/U cl 2 -3 -2 -1 0 1 2 3 0.00 0.02 0.04 0.06 Quinn (’92) A2-1 A2-2 A2-3 y/De=2 x/x_h q 2/U cl 2 -3 -2 -1 0 1 2 3 0.00 0.02 0.04 0.06 0.08 Quinn (’92) A2-1 A2-2 A2-3 y/De=5 x/x_h q 2/U cl 2 -3 -2 -1 0 1 2 3 0.00 0.02 0.04 0.06 0.08 0.10 Quinn (’92) A2-1 A2-2 A2-3 y/De=10 z/z_h q 2/U cl 2 -3 -2 -1 0 1 2 3 0.00 0.02 0.04 0.06 Quinn (’92) A2-1 A2-2 A2-3 y/De=2 z/z_h q 2/U cl 2 -3 -2 -1 0 1 2 3 0.00 0.02 0.04 0.06 0.08 Quinn (’92) A2-1 A2-2 A2-3 y/D_e=5 z/z_h q 2/U cl 2 -3 -2 -1 0 1 2 3 0.00 0.02 0.04 0.06 0.08 0.10 Quinn (’92) A2-1 A2-2 A2-3 y/De=10

Figure 9. Profiles of the turbulent kinetic energy q2_{along the lateral direction (upper figures) and} the spanwise direction at y/De= 2, 5, 10 for the jet with Ar= 2, cases A2-1, A2-2 and A2-3.

StW Amp lit u d e 10-2 ₁₀-1 ₁₀0 ₁₀1 0.00 0.05 0.10 0.15 y/De=0.0 StW Amp lit u d e 10-2 10-1 100 101 0.00 0.05 0.10 0.15 y/De=0.5 StW A m plit ude 10-2 ₁₀-1 ₁₀0 ₁₀1 0.00 0.05 0.10 0.15 y/De=1.0 StW Amp lit u d e 10-2 10-1 100 101 0.00 0.05 0.10 0.15 y/De=3.0

Figure 10. Axial velocity spectrum for the non-excited jet for the case A2-1 at the location of the shear layer x/W= 0.5 at four distances from the inlet y/De= 0.0, 0.5, 1.0, 3.0.

In laboratory experiments [59,60], the preferred mode frequency for a jet with Ar = 2

was found in the range StW, pref= 0.25–0.35. In the present results, the peak at StW= 0.3

is also visible; however, it is not distinctly larger than the fluctuation level in the broad region corresponding to the range of frequencies StW= 0.3–0.7. Hence, from the current

simulations it cannot be said univocally that the peak at StW = 0.3 defines the preferred

frequency. The discrepancies may result from different characteristics of the inlet conditions in the experiments compared to the present studies. It would be very meaningful if a more complete identification of the experimental inlet conditions could be determined to allow further scrutiny of the remaining differences.

4.4. Jet with Ar= 3.

The results obtained for the non-excited jet with Ar = 3 in case A3-1 are presented in

Figure 11. As before, the vorticity modulus and the time-averaged axial velocity in two

cross-section planes are presented. Compared to the results obtained for the jet with Ar=

2, larger vortical structures are formed closer to the inlet, at y/De ≈ 1.75 (De= 1.954W

for Ar= 3). The length of the potential core becomes shorter y/De= 1.5, and, on the other

hand, the crossover location is longer and shifts downstream to y/De≈ 7.

5. Excited jets

The simulations for the non-excited jets established the level of accuracy of the applied LES and correctness of the numerical settings (grid density, spatial disretisation, time integration). It was demonstrated that the inlet boundary conditions defined by Equation (9) allow for a physical representation of the dominant features of jets issuing from square and rectangular nozzles without actual inclusion of the flow upstream of the nozzles in the computational domain. The mean velocity profiles agreed very well with the measurements both for Ar= 1 and Ar= 2. Also, turbulent quantities, e.g., the axial velocity fluctuations

and TKE, were found to be in acceptable agreement with the experimental data. This is

Figure 11. Solution for the non-excited jet with the aspect ratio Ar = 3 for the case A3-1. Left figures: isosurface of the vorticity magnitude (| | = 3(Uj/W)) in the x− y (upper figures) and z − y plane (bottom figures). Right figures: contours of the time averaged axial velocity in the x− y and

z− y plane.

expected to carry over to the simulation results for the excited jets presented in this section, since the dynamic range of scales is not essentially larger compared to the unforced jets.

We concentrate on identifying forcing conditions causing an increased spreading rate of the jet or even a splitting of the jet into two separate branches (bifurcation phenomenon). In round jets, it has been observed that this phenomenon occurs in a wide range of Reynolds numbers. Experiments as well as DNS and LES [25,47,50] studies of round jets confirmed this for Reynolds numbers 1.5 × 103 _{< Re < 10}5_{. Moreover, it was found that this}

phenomenon is only weakly dependent on Re. Necessary conditions leading to a bifurcation were formulated in terms of forcing frequencies yielding a Strouhal criterion Sta/Stf= 2

with 0.35< Sta< 0.7,[25] with Sta, Stfas defined in Equation (11). In parametric studies of

excited circular jets [44], it was shown that the strength of the bifurcation directly depends on Sta. The effect of the forcing was observed to be largest when Stais close to the preferred

mode frequency.

As shown in the previous subsection in the rectangular jets, a specific preferred mode frequency could not be clearly identified, but rather a broader region of more prominent frequencies was observed. For that reason, the simulations are performed for a wide set of frequencies chosen to cover the particular broad region observed inFigure 10, i.e., we consider Sta = 0.3, 0.4, 0.5, 0.6, 0.7, and use a fixed ratio Sta/Stf= 2. In all cases, the

Reynolds number is Re= 1 × 104for which the non-excited jets closely resemble high-Re findings in the near -field (y/De≤ 10). The influence of Re on properties of turbulent square

Table 2. Parameters of the computational cases for the excited jets simulations. The superscript ‘fine’ denotes simulations on denser mesh and ‘ext’ simulations in a wider computational domain. The symbols ‘St0.3-’, ‘St0.4-’ etc. denote the axial forcing frequency and ‘0’ and ‘π₄’ denote the phase shift. Case Ar Sta= faW/Uj  A1-St0.3− 0.7-0,π₄ 1 0.3, 0.4, 0.5, 0.6, 0.7 0,π4 A2-St0.3− 0.7-0,π₄ 2 0.3, 0.4, 0.5, 0.6, 0.7 0,π4 A2fine_-St 0.3− 0.7-0 2 0.3, 0.5, 0.7 0 A2exp_-St 0.3− 0.7-0 2 0.3, 0.5, 0.7 0 A3-St0.3− 0.7-0,π4 3 0.3, 0.4, 0.5, 0.6, 0.7 0, π 4

jets in a range 8 × 103_{≤ Re ≤ 5 × 10}4_{was analysed in recent experimental work by Xu}

et al.[5] It was shown that the near-field region (y/De≤ 10) is not significantly affected by

Re. Results in the far field were found quite independent of Re only beyond Recrit≥ 3 ×

104. In the current simulations, the analysis is limited to the near-field region, where the excitation plays the most important role.

Table 2gives details of the analysed test cases. The notation inTable 2is constructed

as follows: ‘A1-’, ‘A2-’, ‘A3-’ denote the aspect ratio; ‘St0.3-’, ‘St0.4-’ etc. denote the axial

forcing frequency which determines the flapping mode forcing frequency from Sta/Stf= 2;

‘0’ and ‘π

4’ denote the phase shift (see Equation (11)). The computations for the phase

shift =π₄ were motivated by the results obtained in [44] for circular jets. There, the phase shift between the axial and flapping forcing caused qualitative differences. In all cases, the turbulence intensity is equal to Ti= 0.01Ujand the excitation amplitudes are taken as Aa= Af= 0.1Uj. Such a relatively high level of excitation is chosen following [44], where it was

shown that the excitation should have an amplitude larger or at least comparable to the Ti level.

Most of the computations are carried out using the basic mesh. The test computations aiming to asses the influence of the mesh density are performed for the jet with Ar = 2

using the mesh with 320 × 384 × 320 nodes. InTable 2, these computations are denoted by the superscript ‘fine’. Moreover, to check to which extent the side boundaries influence the inner part of the flow domain, an additional set of simulations is carried out using an expanded domain 16W× 16W (in the lateral and spanwise directions). In this case, denoted by the superscript ‘exp’, the dense mesh is used with stretching parameters chosen such that the cell sizes in the vicinity of the jet are smaller than in the basic domain. They are equal to x∗ = 0.017 and z∗ = 0.036 (cf.Table 1). Hence, this case allows to check further the influence of the mesh density in the ‘x–z’ plane. Sample results obtained from these validation tests are shown inFigure 12presenting the radial profiles of the time-averaged axial velocity (mean and fluctuations) at two locations downstream of the inlet. It is seen that the results obtained on the basic mesh and using the basic computational domain (A2-St0.5-0) are only slightly different from those obtained on the denser meshes on the

basic computational domain (A2fine_-St

0.5-0) and on the extended domain (A2exp-St0.5 -0). Further analyses (not presented here) show that small differences are observed only for turbulent quantities (u, v, uv, etc.) and only in the region of high turbulence intensity

y/De= 2.0–4.0. Similar behaviour is found in simulations with Sta= 0.3 and Sta = 0.7.

Influence of the side boundaries of the basic domain on dynamics of the flow in the central part is also very small. Later, in contour plots, it will be shown that for a few cases the forcing leads to very large spreading rate, which causes that, in the region far from the

x/D_e U/ Uj (uu) 1/ 2 /U j 0 1 2 3 0 0.2 0.4 0.6 0.8 1 1.2 0 0.1 0.2 0.3 0.4 0.5 A2-St_0.5- ₀ A2fine -St_0.5- ₀ A2exp -St_0.5- ₀ y/De= 3.0 x/D_e U/ Uj (uu) 1/ 2 /U j 0 1 2 3 0 0.2 0.4 0.6 0.8 1 1.2 0 0.1 0.2 0.3 0.4 0.5 A2-St0.5- 0 A2fine -St0.5- 0 A2exp -St0.5- 0 y/D_e= 7.0

Figure 12. Comparison of the axial velocity profiles at axial distances y/De= 3 and y/De= 7 in the lateral cross-section plane at z= 0. Cases for A2-St0.5-0, A2fine-St0.5-0and A2exp-St0.5-0.

inlet (y/De≈ 10), the jets approach the boundaries of the computational domain. In these

situations, the influence of the boundaries is expected to be significant; however, the far-field region is not of primary importance for the present studies. Hence, one may conclude that the basic mesh and the basic domain are sufficient to obtain accurate numerical results. The dynamics of the excited jets are illustrated in Figures 13and 14 showing the isosurfaces of instantaneous vorticity modulus and the contours of the time-averaged axial velocity for the cases A1-St0.3-0and A2-St0.3-0. It can be seen that almost in the whole

flow domain, the excited jets in the lateral planes are significantly wider compared to the non-excited jets (cf.Figures 4and7). In the non-excited jet for Ar= 2, the crossover point

was found at y/De≈ 4.25, while in the present case, the crossover point does not exist in the

entire computational domain, and the jet in the lateral plane is everywhere wider than in the spanwise plane. The isosurfaces of the vorticity clearly show that the potential core region of the excited jets is very short and strong vortical structures are seen already at y/De≈ 1,

resembling deformed toroidal rings due to the axial forcing. In fact, careful inspection of the zoomed part ofFigures 13and14allows identifying the so-called rib vortices (pointed by the arrows) with elongated shapes originating from the corners of the nozzle. These structures are also present in the non-excited jets but much shorter and less pronounced.

Figure 13. Solution for the excited jet with aspect ratio Ar= 1, the case A1-St0.3-0. Left figure: isosurface of the vorticity magnitude (| | = 3(Uj/W)) in the x–y plane. Right figure: contours of the time-averaged axial velocity in the x–y plane.

Figure 14. Solution for the excited jet with aspect ratio Ar= 2, the case A2-St0.3-0. Left figures: isosurface of the vorticity magnitude (| | = 3(Uj/W)) in the x–y (upper figures) and z–y planes (bottom figures). Right figures: contours of the time-averaged axial velocity in the x–y and z–y planes.

For Ar= 1, hairpin vortices are created in the corners, whereas for the jet with Ar≥ 2, only

single vortical tubes are observed, which is consistent with [12].

The velocity spectra for the excited jets for the cases A2-St0.3-0, A2-St0.5-0 and

A2-St0.7-0are shown inFigure 15. At y/De= 0, two peaks corresponding to the axial and

flapping excitation are very pronounced. Further downstream, a series of higher harmonics

10-2 ₁₀-1 ₁₀0 ₁₀1 10-6 10-4 10-2 100 St_a=0.7 Stf=Sta/2 StW Amp lit u d e 10-2 10-1 100 101 10-6 10-4 10-2 100 _St a=0.3 St_f=St_a/2 10-2 10-1 100 101 10-6 10-4 10-2 100 Sta=0.5 Stf=Sta/2 (a) y/De= 0.0 10-2 ₁₀-1 ₁₀0 ₁₀1 10-6 10-4 10-2 100 St_a=0.7 Stf=Sta/2 10-2 10-1 100 101 10-6 10-4 10-2 100 Sta=0.5 Stf=Sta/2 Amp lit u d e StW 10-2 10-1 100 101 10-6 10-4 10-2 100 _St a=0.3 St_f=St_a/2 (b) y/De= 0.5 10-2 10-1 100 101 10-6 10-4 10-2 100 Sta=0.7 Stf=Sta/2 Amp lit u d e StW 10-2 ₁₀-1 ₁₀0 ₁₀1 10-6 10-4 10-2 100 _St a=0.3 St_f=St_a/2 10-2 ₁₀-1 ₁₀0 ₁₀1 10-6 10-4 10-2 100 Sta=0.5 Stf=Sta/2 (c) y/De= 1.0 10-2 ₁₀-1 ₁₀0 ₁₀1 10-6 10-4 10-2 100 _St a=0.3 Stf=Sta/2 10-2 10-1 100 101 10-6 10-4 10-2 100 Sta=0.7 Stf=Sta/2 Amp lit u d e StW 10-2 ₁₀-1 ₁₀0 ₁₀1 10-6 10-4 10-2 100 _St a=0.5 Stf=Sta/2 (d) y/De= 3.0

Figure 15. Axial velocity spectrum for the excited jet with Ar = 2 at the location of the shear layer x/W= 0.5 at four distances from the inlet y/De= 0.0, 0.5, 1.0, 3.0 for the cases A2-St0.3-0, A2-St0.5-0and A2-St0.7-0corresponding to different Staexplicitly indicated in the figures.

appear with amplitudes that are quite comparable with the forcing amplitudes. The pres-ence of higher frequency harmonics reveals itself already at y/De≈ 0.1 and indicates the

development of smaller flow scales which, at the end of the potential core region, intensify and trigger a sudden transition to turbulent flow, as could also be seen inFigure 14. The spectra obtained for the jets with Ar= 1 and Ar= 3 show qualitatively the same behaviour.

The forcing frequencies are characteristic of the flow even quite far downstream. Higher harmonics are visible mainly in 0.1< y/De< 3.0.

5.1. Time-averaged solutions

The contours of the time-averaged axial velocity in the lateral and spanwise planes obtained from the simulations detailed inTable 2are presented inFigure 16for Ar= 1,Figure 17

for Ar= 2 andFigure 18for Ar= 3. The effects of changing the excitation frequency are

readily noticed for all three aspect ratios. Also, the effect of the phase shift between the

Figure 16. Axial velocity contours in the excited jet with aspect ratio Ar = 1 for the cases A1-St0.3− 0.7-0and A1-St0.3− 0.7-π

4. The results are shown in the lateral plane (first and second rows

of figures) and in the spanwise plane (third and fourth rows).

axial and flapping excitation is pronounced. In the case of the jet with Ar= 1, the results

obtained with = 0 (upper row ofFigure 16) show that for some forcing frequencies, the maxima of the velocity are located in the shear layer regions, while for =π₄, they are always in the centreline, as can be seen in the second row of figures. This effect is even more pronounced in the case with Ar= 2 for the forcing frequency Sta= 0.5 (seeFigure 17). In

the simulation with = 0, the jet almost splits into separate branches, while a phase shift of = π/4 completely prevents this phenomenon. For the jets with Ar= 2, the phase shift

also influences the spreading rate. For instance, in simulations with Sta= 0.7 and  = π₄,

the spread of the jet is smaller than in case = 0. On the other hand, in the simulations with Sta = 0.3, the opposite trend is observed. The results obtained for the jets Ar= 3,

shown inFigure 18, are only slightly dependent on the phase shift and the solutions for

 = 0 and  = π

4 differ only quantitatively. It is surprising that the tendency to splitting

clearly observed for Ar= 2 is practically not seen at all for Ar= 3.

Figure 17. Axial velocity contours in the excited jet with aspect ratio Ar = 2 for the cases A2-St0.3− 0.7-0and A2-St0.3− 0.7-π

4. The results are shown in the lateral plane (first and second rows

of figures) and in the spanwise plane (third and fourth rows).

The simulations show that for all values of Sta the jets are wider in the lateral plane

in which the flapping excitation is applied compared to the spanwise plane. This effect increases with increasing Ar. In case Ar= 1, although the flow structures of the jets in the

lateral and spanwise planes differ, their sizes are quite similar, whereas in the simulations with Ar= 2 and Ar= 3, the jets in the lateral plane are significantly wider. In many cases,

the contours of the velocity even reach the side boundaries, for instance, in the simulation A2-St0.3-π₄ shown inFigure 17in the second row of the first column. In these situations,

the results far from the inlet are somewhat biased by the influence of the nearby boundaries. Such a large spread of the jet in the lateral plane implies a corresponding reduction of the jet spreading in the spanwise direction. It is seen that in some cases, for instance, A2-St0.5-0

or A3-St0.4-0, in the region y/De< 2, the jets first expand and then suddenly narrow to

a size smaller even than its inlet dimension. Analysis of the vector fields (not presented here) shows that in this region the velocity vectors are oriented towards the centreline. This

Figure 18. Axial velocity contours in the excited jet with aspect ratio Ar = 3 for the cases A3-St0.3− 0.7-0 and A3-St0.3− 0.7-π

4. The results are shown in the lateral plane (1st and 2nd row of

figures) and in the spanwise plane (3rd and 4th row).

illustrates the remarkable level of control that can be achieved on the far-field development of the flow.

The axial velocity profiles are shown inFigures 19–21. These profiles correspond to the solutions at distances y/De = 1, y/De = 3 and y/De = 7 along the lateral direction

(sub-figures (a), (c) and (e)) and along the spanwise direction ((b), (d) and (f)). At y/De

= 1, the effect of excitation is small and the results are very similar regardless of Sta.

This means that the differences between the solutions observed further downstream are not directly related to the mechanical forcing, but rather indirectly, resulting mainly from interactions between natural instability modes and the flow structures created by the forcing at different Sta. A clear illustration of the effect of varying Staappears beyond y/De≈ 3,

where the profiles exhibit quite complicated patterns with a flattening in the lateral plane around x/De= 1. This effect is most pronounced in the jets with Ar= 1.

The solutions at y/De = 7 clearly confirm that the excitation leads to a significant

widening of the jet in the lateral direction and narrowing in the spanwise direction, compared

x/D_e U/ Uj 0 1 2 3 0 0.2 0.4 0.6 0.8 1 1.2 no forcing St_a=0.3 Sta=0.4 Sta=0.5 Sta=0.6 Sta=0.7 y/De=1.0 (a) z/D_e U/ Uj 0 1 2 3 0 0.2 0.4 0.6 0.8 1 1.2 no forcing St_a=0.3 Sta=0.4 Sta=0.5 Sta=0.6 Sta=0.7 y/De=1.0 (b) x/D_e U/ Uj 0 1 2 3 0 0.2 0.4 0.6 0.8 1 1.2 no forcing St_a=0.3 St_a=0.4 Sta=0.5 Sta=0.6 Sta=0.7 y/De=3.0 (c) z/D_e U/ Uj 0 1 2 3 0 0.2 0.4 0.6 0.8 1 1.2 no forcing St_a=0.3 St_a=0.4 Sta=0.5 Sta=0.6 Sta=0.7 y/De=3.0 (d) x/D_e U/ Uj 0 1 2 3 0 0.2 0.4 0.6 0.8 1 1.2 no forcing St_a=0.3 St_a=0.4 Sta=0.5 Sta=0.6 Sta=0.7 y/De=7.0 (e) z/D_e U/ Uj 0 1 2 3 0 0.2 0.4 0.6 0.8 1 1.2 no forcing St_a=0.3 St_a=0.4 Sta=0.5 Sta=0.6 Sta=0.7 y/De=7.0 (f)

Figure 19. Mean axial velocity profiles for the jet with Ar= 1 for the non-excited case and for the cases A1-St0.3–0.7-0. Solutions at axial distances y/De= 1, y/De = 3 and y/De = 7 in the lateral cross-section plane at z= 0 (figures in the column on the left side) and in the spanwise cross-section plane at x= 0.

x/D_e U/ Uj 0 1 2 3 0 0.2 0.4 0.6 0.8 1 1.2 no forcing Sta=0.3 Sta=0.4 Sta=0.5 St_a=0.6 St_a=0.7 y/De=1.0 (a) z/D_e U/ Uj 0 1 2 3 0 0.2 0.4 0.6 0.8 1 1.2 no forcing Sta=0.3 Sta=0.4 Sta=0.5 St_a=0.6 St_a=0.7 y/De=1.0 (b) x/D_e U/ Uj 0 1 2 3 0 0.2 0.4 0.6 0.8 1 1.2 no forcing Sta=0.3 Sta=0.4 Sta=0.5 Sta=0.6 Sta=0.7 y/De=3.0 (c) z/D_e U/ Uj 0 1 2 3 0 0.2 0.4 0.6 0.8 1 1.2 no forcing Sta=0.3 Sta=0.4 Sta=0.5 Sta=0.6 Sta=0.7 y/De=3.0 (d) x/D_e U/ Uj 0 1 2 3 0 0.2 0.4 0.6 0.8 1 1.2 no forcing St_a=0.3 Sta=0.4 Sta=0.5 Sta=0.6 Sta=0.7 y/De=7.0 (e) z/D_e U/ Uj 0 1 2 3 0 0.2 0.4 0.6 0.8 1 1.2 no forcing St_a=0.3 Sta=0.4 Sta=0.5 Sta=0.6 Sta=0.7 y/De=7.0 (f)

Figure 20. Mean axial velocity profiles for the jet with Ar= 2 for the non-excited case and for the cases A2-St0.3–0.7-0. Solutions at axial distances y/De= 1, y/De= 3 and y/De= 7 in the lateral cross-section plane at z= 0 (figures in the column on the left side) and in the spanwise cross-section plane at x= 0.

x/D_e U/ Uj 0 1 2 3 0 0.2 0.4 0.6 0.8 1 1.2 no forcing St_a=0.3 Sta=0.4 Sta=0.5 Sta=0.6 Sta=0.7 y/De=1.0 (a) z/D_e U/ Uj 0 1 2 3 0 0.2 0.4 0.6 0.8 1 1.2 no forcing St_a=0.3 Sta=0.4 Sta=0.5 Sta=0.6 Sta=0.7 y/De=1.0 (b) x/D_e U/ Uj 0 1 2 3 0 0.2 0.4 0.6 0.8 1 1.2 no forcing St_a=0.3 St_a=0.4 Sta=0.5 Sta=0.6 Sta=0.7 y/De=3.0 (c) z/D_e U/ Uj 0 1 2 3 0 0.2 0.4 0.6 0.8 1 1.2 no forcing St_a=0.3 St_a=0.4 Sta=0.5 Sta=0.6 Sta=0.7 y/De=3.0 (d) x/D_e U/ Uj 0 1 2 3 0 0.2 0.4 0.6 0.8 1 1.2 no forcing St_a=0.3 St_a=0.4 Sta=0.5 Sta=0.6 Sta=0.7 y/De=7.0 (e) z/D_e U/ Uj 0 1 2 3 0 0.2 0.4 0.6 0.8 1 1.2 no forcing St_a=0.3 St_a=0.4 Sta=0.5 Sta=0.6 Sta=0.7 y/De=7.0 (f)

Figure 21. Mean axial velocity profiles for the jet with aspect ratio Ar= 3 for the non-excited case and for the cases A3-St0.3–0.7-0. Solutions at axial distances y/De= 1, y/De= 3 and y/De= 7 in the lateral cross-section plane at z= 0 (figures in the column on the left side) and in the spanwise cross-section plane at x= 0.

to the non-excited case. For instance, in case Ar= 3, the profiles presented inFigure 21(e)

show that at the lateral distance x/De≈ 1.7, the velocity of the unexcited jet approaches

the co-flow velocity, while in all excited cases (A3-St0.3–0.7-0), it is comparable with the

centreline velocity. On the other hand, in the spanwise direction, the jets are significantly thinner, compared to the non-excited case. The velocities of excited jets reaches zero at

z/De≈ 1.5, where the velocity of non-excited jet still has substantial value of about 20% of Uj.

5.1.1. Bifurcation phenomenon

The results presented inFigures 16–18show that the bifurcation phenomenon as observed for round jets, strictly speaking, does not appear in any of the presented cases. In the round jets [25,44,47,50] and also in the square jets [9,11] at low Reynolds numbers, a characteristic effect of the bifurcation phenomenon was a clear separation of the jets into two distinct branches and almost total vanishing of the flow near the axis. In those simulations, the jets exhibited a ‘Y’ shape when seen as a cross section of time-averaged velocity. In the present simulations, there are only some indications of this behaviour, as may be seen, for instance, in the results for A2-St0.5-0 inFigures 17and20(c), but a distinct splitting of

the main stream does not occur. The main effect appears a strongly anisotropic spreading rate at some forcing conditions, e.g., a much stronger spreading rate in the lateral compared to the spanwise direction. This illustrates the qualitative modification of the flow that is achievable via the inflow forcing.

6. Conclusions

The paper presented the results of LES of non-excited and excited square and rectangular jets with aspect ratios Ar = 1, Ar = 2 and Ar = 3. The simulations performed for

non-excited jets validated the applied numerical approach (subgrid modelling, mesh densities, boundary conditions) and assessed the dependency of the solutions on the parameters, i.e., Reynolds number and inlet turbulence time/length scales. Importance of the inlet turbulence characteristics was demonstrated comparing the axial velocity along the jet centrelines. This allowed to chose the time and length scales to best match experimental data in the region downstream of the inlet. The simulation results were in good agreement with the experimental data taken from the literature. In the downstream region y/De< 10, the profiles

of the mean axial velocity almost exactly matched the measurements. Discrepancies were seen for the fluctuating components, which points toward differences in small-scale features of the inflow conditions between the experiments and the simulations. It turned out that the Reynolds number in the range 1 × 104_{to 1.8 × 10}5 _{has only a minor influence on the}

near-field results.

The computations showed that active and passive flow control methods (excitation and nozzle shaping) may be successfully combined and used to increase the mixing of the jets or to alter their behaviour even quite far downstream of the inflow. It was shown that a rich variety of mean flow patterns can be induced simply by changing aspect ratio and forcing parameters. Modifications of the forcing frequencies allowed alteration of the large-scale flow in the region downstream of the inlet boundary. Analysis of the solutions in the spectral space showed that at a distance up to y/De≤ 3, the induced axial and flapping forcing result

in smaller scales which were identified by high-frequency harmonics. In this region, the

contour plots and also one-dimensional plots revealed a very complex flow pattern. Further downstream, the observed effects of excitation were as follows:

• Increased spreading of the jet in the lateral plane—significantly larger than in the case of the unexcited jet; the widening of the jet in the lateral direction caused its narrowing in the spanwise plane and decreasing of the velocity on the centreline. • For some values of the excitation frequency, the jets exhibited a tendency to the

separation of the main flow direction; the occurrence of this phenomenon for the jets with Ar = 1 and Ar = 2 was conditioned by the phase shift between the axial

and flapping forcing; when the forcing terms were in phase, the maxima of the axial velocity could be found away from the centreline of the jet which is characteristic for the bifurcation phenomenon; when forcing terms were shifted byπ/4, this effect disappeared. For the jet with Ar = 3, the phase shift was found to lead to only

quantitative changes in the flow; large qualitative alterations in the flow structure were not observed.

• The full splitting of the jet, in the way as for the round jets, was not observed in the present simulations; earlier works [9,11] showed that the jets splitted into completely separate branches; however, these results were reported for low Reynolds number flows, Re=1000 and Re=3200, and appear not to be robust enough to survive higher Reynolds numbers.

In simulations of round jets [44], the excitation amplitude had to be large compared to the level of turbulence intensity in order to induce a strong bifurcation. Taking into account that the rectangular jets are more unstable than the round jets, it is likely that the amplitude of excitation will play a key role on the bifurcation in two or more jets at high Reynolds numbers. Further research in this direction is planned.

Disclosure statement

No potential conflict of interest was reported by the authors.

Funding

This work has been partially supported by the Polish National Science Center [grant number DEC-2011/03/B/ST8/06401]. Computations have been carried out at SARA Computing Centre (Amster-dam) [grant number SH-061]; Cyfronet Computing Centre (Krakow) within the PL-Grid infrastruc-ture.

References

[1] Kral LD. Active flow control technology. New York (NY): ASME Fluids Engineering Division; 1998. (Technical brief).

[2] Jahanmiri M. Active flow control: a review. G¨oteborg (Sweden) Chalmers University of Tech-nology; 2010. (Technical report 2010:12).

[3] Quinn WR, Militzer J. Experimental and numerical study of a turbulent free square jet. Phys Fluids. 1988;31:1017–1025.

[4] Quinn WR. Turbulent free jet flows issuing from sharp-edged rectangular slots: the influence of slot aspect ratio. Exp Thermal Fluid Sci. 1992;5:203–215.

[5] Xu M, Pollard A, Mi J, Secretain F, Sadeghi H. Effects of Reynolds number on some properties of a turbulent jet from a long square pipe. Phys Fluids. 2013;25:1–19.