Computational Characterization of Turbulent Flow in a Microﬂuidic Actuator

: In this contribution, an unsteady numerical simulation of the ﬂow in a microﬂuidic oscillator has been performed. The transient turbulent ﬂow inside the device is described by the Unsteady Reynolds Averaged Navier–Stokes equations (URANS) coupled with proper turbulence models. The main characteristics of the complex ﬂuid ﬂow inside the device along one oscillation cycle was analyzed in detail, including not only velocity contours but also the pressure and turbulent kinetic energy ﬁelds. As a result, two-dimensional simulations provided good estimations of the operating frequency of the ﬂuidic actuator when compared with experimental measurements in a range of Reynolds numbers. Moreover, with the objective of altering the operating frequency of the apparatus and, in order to adapt it to different applications, geometrical modiﬁcations of the feedback channels were proposed and evaluated. Finally, a fully three-dimensional simulation was carried out, which allowed for the identiﬁcation of intricate coherent structures revealing the complexity of the turbulent ﬂow dynamics inside the ﬂuidic oscillator.


Introduction
Fluidic oscillators are devices that generate periodic jet pulses of constant frequency from a steady supply of gas or liquid. These devices can be monolithic in their construction, without moving parts, capable of withstanding high temperatures and pressuresdepending on their construction material, long runs without maintenance, etc. These qualities make fluidic oscillators simpler and more reliable than traditional flow control devices. The fluidic oscillation happens naturally in a self-induced and self-sustained manner, which means that no input from the user is required to make the device function and run indefinitely, as long as the inlet flow is on.
Analogous to electronic oscillators, fluidic oscillators must embody four essential elements: a continuous supply of energy-which is provided by the inlet flow-delivered at the right intensity to compensate for the hydraulic losses as well as turbulent dissipation, a feedback mechanism-which is provided by the feedback loops-a frequency determining mechanism-which is provided by the length of the feedback loops-and an amplifier, which constitutes the mixing chamber. The frequency-determining mechanism delivers flow to the amplifier, along with the main flow, which in turn must produce enough gain so part of its flow is fed back to the feedback loops to sustain the oscillations and of course produce enough output flow. As with electronic oscillators, the purpose of a fluidic oscillator should be to deliver a periodic steady flow of fluid at a constant frequency and constant amplitude, although it depends on the specific application.
Fluidic oscillators come in a variety of geometrical forms, working principles, materials and applications. Tesar [1] published an extensive review of fluidic oscillator topologies and Ghanami and Farhadi [2] reviewed fluidic oscillator geometries, fabrication techniques, pologies and Ghanami and Farhadi [2] reviewed fluidic oscillator geometries, fabrication techniques, materials, applications and mechanisms. Guzmán et al. [3] have made an elaborated description of the functioning of a fluidic oscillator device.
In recent years, there has been a flurry of activity both in the design of new fluidic oscillators and new applications. Fluidic oscillators have found novel applications in the flow control of swirl flame dynamics [4], fuel combustion in gas turbines [5], control of fluid flow for aerodynamics applications [6], ultrasonic testing [7], jet actuators for flow control [8], fuel injection in hydrogen engines [9], even in controlling the flow in fume cupboards in chemistry laboratories [10]. One very interesting application is that presented in [11], where a fluidic oscillator was coupled with a low-temperature plasma dielectric barrier discharge device to produce ozone-rich microbubbles as a cost-effective means to lyse cyanobacterial cells and degrade the cyanotoxins they produce. The use of microbubbles is advantageous as heat and mass transfer are enhanced as the size of the bubbles is decreased; in addition, the use of ozone produced in a microplasma has energy, capital investment and maintenance cost savings that make the use of advanced oxidation technologies more attractive for water treatment. Lozano and Zimmerman [12] established the kinetic conditions under which ozone can be cheaply and efficiently produced in a microplasma, and Tesar [13] made an extensive review of the state of the art in fluidic oscillators for microbubble generation.
The fluidic oscillator considered in this work is that studied experimentally in Bobusch et al. [14] and numerically in Krüger et al. [15]. This device consists of an inlet nozzle, also known as a power nozzle, which injects a liquid (water) into a mixing chamber-also known as an interacting region. It has two feedback loops that link the rear part of the mixing chamber with its front part, and an outlet port connected to two expanding branches with outlets exposed to ambient pressure (see Figure 1 below), through which the oscillating flow leaves the device. The experiments employed planar particle image velocimetry (PIV) to obtain time-resolved velocity fields of the fluid inside the actuator. From the numerical point of view, the commercial software ANSYS-CFX was used to perform URANS (Unsteady Reynolds Averaged Navier-Stokes) simulations aimed at describing the turbulent flow time evolution. Comparison of the two-dimensional (2D) computed and measured contour velocity fields was satisfactory; moreover, the smallest difference between the numerical and experimental operational frequency of the fluid oscillator was found when the k − ω Shear Stress Transport (SST) turbulence model was applied. Therefore, the authors conclude that 2D computations are sufficient to capture the principal features of the flow inside the oscillator.  [14] illustrating its parts as well as the boundary conditions. More recently, Baghaei et al. [16] performed three-dimensional simulations of the turbulent flow in the previous oscillator. That study focused on the evaluation of a few geometric modifications aimed at altering the oscillator's operating frequency. They used a hybrid Detached Eddy Simulation (DES) turbulence model in OpenFoam, which is an  [14] illustrating its parts as well as the boundary conditions. More recently, Baghaei et al. [16] performed three-dimensional simulations of the turbulent flow in the previous oscillator. That study focused on the evaluation of a few geometric modifications aimed at altering the oscillator's operating frequency. They used a hybrid Detached Eddy Simulation (DES) turbulence model in OpenFoam, which is an open CFD software, to simulate the flow characteristics. The width of the inlet and exit sections of the mixing chamber as well as the angles of the inner walls of the exit nozzle were modified, and curves of frequency and amplitude of the oscillation versus the corresponding dimension values were produced. It was shown that the frequency can be increased by augmenting the width of the exit section, whereas the opposite effect was found when it was reduced. On the other hand, it was found that the inlet of the mixing chamber needed to be within a certain range of values in order to produce oscillations; too large or too small widths caused the jet to go straight through the mixing chamber without oscillation. Therefore, the conclusion was that, for a fixed Reynolds number, the oscillation frequency can be effectively altered by modifying some internal dimensions of the device; in that respect, the operation of the oscillator can be adapted for different applications. However, such modifications would not be realizable without introducing moving parts in the design, which defeats one of the main advantages of the device, which is that related to being a monolithic piece without moving parts.
Motivated by the previous findings, in this work we focus on modifying the feedback loops in [14] rather than the inlet or outlet widths, since with such modifications we can still have a monolithic device, with all the advantages it entails. Therefore, CFD simulations have been chosen for evaluating the effect of such modifications due to their ability of handling multiple conditions in a cheap and efficient way, reducing the developing times and costs of a multitude of engineering devices and processes [17]. Moreover, some preliminary results obtained by a fully three-dimensional simulation are presented, illustrating the fascinating complexity of the turbulent flow dynamics inside the fluidic oscillator.

Geometrical Configuration and Simulation Set-Up
The geometry of the fluidic oscillator considered in this work is that investigated in [14,15], which is presented in Figure 1, showing its constitutive parts and also the employed boundary conditions.
A uniform velocity profile is imposed at the inlet, while a zero-pressure condition is set at the outlets; finally, at the walls, the no-slip condition is enforced. According to the experimental conditions in [14], the working fluid is water with density ρ = 998.2 kg/m 3 and kinematic viscosity ν = 1.004 × 10 −6 m 2 /s; therefore, its dynamics will be described by the transient incompressible Navier-Stokes equations.
The principle of operation of the fluidic oscillator is summarized as follows. The flow is injected at the inlet with constant velocity; it accelerates in the nozzle (also called power nozzle) due to the reduction in cross section and enters the mixing chamber. In principle, the generated jet goes straight to the end of the mixing chamber, crosses the exit nozzle and encounters the splitter, which divides the flow through the branches towards both outlets. This flow configuration is unstable, and the jet starts to perform oscillations in the cross flow direction. Such oscillations soon increase their amplitude, and the jet collides either with the upper or lower wall at the end of the mixing chamber with the effect that the flow divides, and part of it is directed through one of the feedback channels. The flow through that feedback channel reaches the beginning of the mixing chamber, pushing the main jet towards the opposite wall, generating the self-sustained jet oscillation dynamics.
The flow Reynolds number based on the hydraulic diameter D h is Re ≈ 16, 000 [15]. This means that the flow inside the device is turbulent, which in this study is described using turbulence models under the Unsteady Reynolds Averaged Navier-Stokes (URANS) approximation. In the first instance and for validation purposes, the Shear Stress Transport (SST) k-ω model [18] was used, as in Krüger et al. [15]. This model combines both k-ε and k-ω models through a blending function. The k-ε model is used far away from the wall, whilst the k-ω model is used near the wall up to the edge of the boundary layer. The SST kω model performs better than other two-equation classical turbulence models in the case of flows with strong adverse pressure gradients and flow separation and reattachment, which are situations commonly found in fluidic oscillators. The governing equations of such a model can be found in the original paper of Menter [18] and will not be repeated here.
Since the fluidic oscillator has a planar geometry, two-dimensional simulations were performed with the ANSYS Fluent v. 19 commercial software. Preliminary results of a 3D simulation will also be presented in Section 6. The 2D geometrical domain was discretized by means of a structured grid based on quads. To perform the mesh independence study, three grids have been considered in this work: a coarse one with around 10 5 elements, intermediate 3.43 × 10 5 cells and refined with 11.75 × 10 5 elements. Figure 2a shows a panoramic view of the intermediate mesh, whilst Figure 2b is a zoom-in of the red rectangle enclosing the mixing chamber: there, it is seen that the grid is pretty uniform. In this graph it can also be seen that the grid is substantially refined close to the walls in order to resolve the boundary layer, giving maximum values of y + ∼ 1 (Table 1), which is a requirement of the turbulence model used.
separation and reattachment, which are situations commonly found in fluidic oscillators. The governing equations of such a model can be found in the original paper of Menter [18] and will not be repeated here.
Since the fluidic oscillator has a planar geometry, two-dimensional simulations were performed with the ANSYS Fluent v. 19 commercial software. Preliminary results of a 3D simulation will also be presented in Section 6. The 2D geometrical domain was discretized by means of a structured grid based on quads. To perform the mesh independence study, three grids have been considered in this work: a coarse one with around 10 5 elements, intermediate 3.43 × 10 5 cells and refined with 11.75 × 10 5 elements. Figure 2a shows a panoramic view of the intermediate mesh, whilst Figure 2b is a zoom-in of the red rectangle enclosing the mixing chamber: there, it is seen that the grid is pretty uniform. In this graph it can also be seen that the grid is substantially refined close to the walls in order to resolve the boundary layer, giving maximum values of ~1 (Table 1), which is a requirement of the turbulence model used.  The employed spatial discretization schemes for all the differential equations were the second-order centered scheme for diffusive terms and the second-order upwind for the convective terms. Discretization in time used an implicit second-order scheme with a time step determined after performing a sensitivity study, which is presented below. Pressure-velocity coupling is handled by the transient Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) scheme. At each time step, a maximum of 200 iterations has been fixed with a residual convergence criterion of 10 −7 , which is typically attained at every time step.
The simulation approach starts with a steady-state simulation to obtain an initial estimation of the flow field. Then, the transient simulation is started with first-order schemes to establish the oscillation; in the next step, the spatial discretization schemes are switched to second order and, finally, the temporal scheme is also set to second order to improve accuracy. The simulation is finally run for five complete oscillation periods in order to extract the simulation results.
The grid independency study was carried out taking the mass flow in outlet 1 as the monitored variable, which is normalized with the inlet mass flow. Figure 3 shows the time evolution of this variable for the three considered meshes. It is important to notice that in this figure, positive values of mass flow indicate fluid entrance to the actuator through the corresponding surface, and negative values mean that the fluid is exiting the device. It can be seen that the maximum difference of the mass flow signal between the  The employed spatial discretization schemes for all the differential equations were the second-order centered scheme for diffusive terms and the second-order upwind for the convective terms. Discretization in time used an implicit second-order scheme with a time step determined after performing a sensitivity study, which is presented below. Pressurevelocity coupling is handled by the transient Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) scheme. At each time step, a maximum of 200 iterations has been fixed with a residual convergence criterion of 10 −7 , which is typically attained at every time step.
The simulation approach starts with a steady-state simulation to obtain an initial estimation of the flow field. Then, the transient simulation is started with first-order schemes to establish the oscillation; in the next step, the spatial discretization schemes are switched to second order and, finally, the temporal scheme is also set to second order to improve accuracy. The simulation is finally run for five complete oscillation periods in order to extract the simulation results.
The grid independency study was carried out taking the mass flow in outlet 1 as the monitored variable, which is normalized with the inlet mass flow. Figure 3 shows the time evolution of this variable for the three considered meshes. It is important to notice that in this figure, positive values of mass flow indicate fluid entrance to the actuator through the corresponding surface, and negative values mean that the fluid is exiting the device. It can be seen that the maximum difference of the mass flow signal between the intermediate and refined grids is lower than 1%; therefore, the intermediate grid has been chosen as a compromise between accuracy and total computation time. This will be denoted hereafter as base case. intermediate and refined grids is lower than 1%; therefore, the intermediate grid has been chosen as a compromise between accuracy and total computation time. This will be denoted hereafter as base case. Furthermore, a time-step sensitivity study has been performed for the base case. Two small time steps ∆ = 10 , 10 s have been tested. It can be seen in Figure 4 that there are no significant differences between the two mass flow signals, so a value of 10 s has been set for the time step.

Validation of Simulations
As previously mentioned, the numerical simulation of the flow in the fluidic oscillator was carried out with the commercial software ANSYS-Fluent v. 19. Computations are validated versus experimental data [14] considering the resulting oscillation frequency at four Reynolds numbers; the comparison is shown in Figure 5. In this figure, the single result in [15] for ≈ 16,000, using the SST k-ω turbulence model, is also included, as well as the computed frequencies obtained in [16] using 3D simulations and a hybrid DES turbulence model, which were performed with the software OpenFoam. Furthermore, a time-step sensitivity study has been performed for the base case. Two small time steps ∆t = 10 −5 , 10 −6 s have been tested. It can be seen in Figure 4 that there are no significant differences between the two mass flow signals, so a value of 10 −5 s has been set for the time step.
Appl. Sci. 2022, 12, x FOR PEER REVIEW 5 of 17 intermediate and refined grids is lower than 1%; therefore, the intermediate grid has been chosen as a compromise between accuracy and total computation time. This will be denoted hereafter as base case. Furthermore, a time-step sensitivity study has been performed for the base case. Two small time steps ∆ = 10 , 10 s have been tested. It can be seen in Figure 4 that there are no significant differences between the two mass flow signals, so a value of 10 s has been set for the time step.

Validation of Simulations
As previously mentioned, the numerical simulation of the flow in the fluidic oscillator was carried out with the commercial software ANSYS-Fluent v. 19. Computations are validated versus experimental data [14] considering the resulting oscillation frequency at four Reynolds numbers; the comparison is shown in Figure 5. In this figure, the single result in [15] for ≈ 16,000, using the SST k-ω turbulence model, is also included, as well as the computed frequencies obtained in [16] using 3D simulations and a hybrid DES turbulence model, which were performed with the software OpenFoam.

Validation of Simulations
As previously mentioned, the numerical simulation of the flow in the fluidic oscillator was carried out with the commercial software ANSYS-Fluent v. 19. Computations are validated versus experimental data [14] considering the resulting oscillation frequency at four Reynolds numbers; the comparison is shown in Figure 5. In this figure, the single result in [15] for Re ≈ 16, 000, using the SST k-ω turbulence model, is also included, as well as the computed frequencies obtained in [16] using 3D simulations and a hybrid DES turbulence model, which were performed with the software OpenFoam.
In Figure 5, the error bars, which are superimposed to the experimental points, indicate a 5% deviation for each particular value. As can be seen, the present computations lie within the error bars, although the frequency is always slightly overpredicted. The 3D simulations in [16] show the same trend as the actual two-dimensional simulations and, in fact, the values are very close to each other, with maximum differences between them of around 1%. This result suggests that the 2D simulations suffice to reproduce the global flow features in the oscillator and, in particular, its operation frequency. Finally, in Figure 5, the frequency obtained in Krüger et al. [15] at Re ≈ 16, 000 is also displayed, which in this case is below the experimental measurement but within the 5% error bars. It is necessary to stress that the present simulations and those of [15] are two-dimensional and use the same turbulence model, but were performed with different software (Fluent in this work and CFX in [15]).
In Figure 5, the error bars, which are superimposed to the experimental points, indicate a 5% deviation for each particular value. As can be seen, the present computations lie within the error bars, although the frequency is always slightly overpredicted. The 3D simulations in [16] show the same trend as the actual two-dimensional simulations and, in fact, the values are very close to each other, with maximum differences between them of around 1%. This result suggests that the 2D simulations suffice to reproduce the global flow features in the oscillator and, in particular, its operation frequency. Finally, in Figure  5, the frequency obtained in Krüger et al. [15] at ≈ 16,000 is also displayed, which in this case is below the experimental measurement but within the 5% error bars. It is necessary to stress that the present simulations and those of [15] are two-dimensional and use the same turbulence model, but were performed with different software (Fluent in this work and CFX in [15]).
An interesting exercise is to plot the normalized mass flow at outlet 1 versus a non-dimensional time obtained by dividing the real time over the period of oscillation for the four Reynolds numbers simulated. As can be seen in Figure 6, the four curves overlap each other approximately, which means that the non-dimensional frequency of the fluidic actuator is virtually independent of for the studied range of Reynold numbers, between 8000 and 16,000.   [14] and the CFD numerical simulations in Krüger et al. [15] and Baghaei et al. [16].
An interesting exercise is to plot the normalized mass flow at outlet 1 versus a nondimensional time obtained by dividing the real time over the period of oscillation for the four Reynolds numbers simulated. As can be seen in Figure 6, the four curves overlap each other approximately, which means that the non-dimensional frequency of the fluidic actuator is virtually independent of Re for the studied range of Reynold numbers, between 8000 and 16,000.  [14] and the CFD numerical simulations in Krüger et al. [15] and Baghaei et al. [16].
In Figure 5, the error bars, which are superimposed to the experimental points, indicate a 5% deviation for each particular value. As can be seen, the present computations lie within the error bars, although the frequency is always slightly overpredicted. The 3D simulations in [16] show the same trend as the actual two-dimensional simulations and, in fact, the values are very close to each other, with maximum differences between them of around 1%. This result suggests that the 2D simulations suffice to reproduce the global flow features in the oscillator and, in particular, its operation frequency. Finally, in Figure  5, the frequency obtained in Krüger et al. [15] at ≈ 16,000 is also displayed, which in this case is below the experimental measurement but within the 5% error bars. It is necessary to stress that the present simulations and those of [15] are two-dimensional and use the same turbulence model, but were performed with different software (Fluent in this work and CFX in [15]).
An interesting exercise is to plot the normalized mass flow at outlet 1 versus a non-dimensional time obtained by dividing the real time over the period of oscillation for the four Reynolds numbers simulated. As can be seen in Figure 6, the four curves overlap each other approximately, which means that the non-dimensional frequency of the fluidic actuator is virtually independent of for the studied range of Reynold numbers, between 8000 and 16,000.

Flow Visualization and Analysis
The flow development along the fluidic oscillator through a complete oscillation cycle is shown in Figure 7 through velocity and pressure contours. Figure 7 shows on the left column a contour plot of the velocity field at certain instants of the oscillation cycle, where the last row (Figure 7g-n) is the flow field shown in the upper row (Figure 7a-h) after completing one cycle; the right column of Figure 7 shows the corresponding pressure contour plots with superimposed velocity vectors, whose relative length is proportional to the velocity magnitude. The behavior of the velocity field has been extensively described from an experimental point of view [14] and it will not be repeated here; thus, only a brief synopsis will be given. Looking at the right column of Figure 7, it is noticed that pressure grows in the mixing chamber when the jet impinges on the exit nozzle walls; it is maximum when the jet impacts the exit nozzle wall with the steepest angle (e.g., Figure 7b,i) and minimum when the flow directly reaches the exit nozzle. The high pressure values at the impinging points are clearly noticed in Figure h,i,k,l, and the large vortex formed in the outlet branches in Figure 7h,l presents the minimum pressure values. In this way, an oscillation pattern of the pressure field is developed, which is connected with the mass flow cycle (in fact, they have opposite phases), a fact that has been referred to as "axial pumping" [14] in the fluid actuator. Interestingly, previous studies usually report only the behavior of the velocity field, so that the present work is one of the first studies to show explicitly the pressure field cycle inside this kind of device. Figure 8 displays the turbulent kinetic energy contour field along one oscillation cycle at the corresponding times of Figure 7. It can be noticed that the highest values of are attained in the mixing layers generated by the jet in the mixing chamber and in the outlet branches, as well as in the recirculation zones that flank the jet between the upper and lower blocks. It is precisely in such areas where large velocity gradients are generated, leading to substantial turbulent kinetic energy production. However, in the surroundings where the jet impinges the exit nozzle walls, adopts moderate values because in such zones, turbulent dissipation is high. The large vortical structures developing in the mixing chamber and in the outlet branches are zones of high values of turbulent kinetic energy (see Figure 8a,c,d,f,g), as already observed in the study Guzmán et al. [3], which also exhibits low pressure values. In Figure 8b  Overall, it can be said that the main flow characteristics are reproduced by the present numerical simulations (compare Figure 7, left column, with Figure 9 of [14]). In summary, the jet flow inside the mixing chamber shows a meandering behavior in which the main jet is flanked by two vortices whose size augments and reduces as the oscillation proceeds. The jet impinging on the inner walls of the exit nozzle-on the lower wall (Figure 7a,b,h,i) in the first part of the cycle and on the upper wall (Figure 7c,d,j,k) in the second part-divides the flow towards the feedback loops and the exit nozzle. Flow in the recirculation channels is complex, displaying several separation zones that force the fluid to follow defined paths (see light green and blue areas along the loops in Figure 7a) resembling those of a car in a racing track. The pressure field in the mixing chamber in such situations is consequently positive, due to the existence of a stagnation point in one of the exit walls (which presents the maximum value of pressure); however, downstream of the exit nozzle, the pressure field shows mainly negative values. Referring to Figure 3, the incoming flow at outlet 1 reaches a maximum when the exit flow in the second outlet is also maximum. In that moment the pressure inside outlet branch 1 is markedly negative and its flow is characterized by a large and intense vortex (Figure 7d,k), which causes the ingestion of fluid by outlet 1. Such fluid is eventually conducted towards the other branch to be expelled from the domain. The roles of branches 1 and 2 are interchanged when the exiting flow through outlet 1 reaches a maximum, which is reflected in the minimum in the curve of Figure 3, because flow out of the domain is considered negative (Figure 7b,i). Therefore, the outlet branches exhibit net flow ingestion during roughly one third of the oscillation period.
On the other hand, the fluid entering the feedback loops from the rear of the mixing chamber emerges just downstream of the exit of the inlet (or power) nozzle pushing the main jet towards the opposite side (see Figure 7b,e), thus self-sustaining the oscillation cycle. In the flow through such feedback loops, some short-circuits can be noticed, due to the generation of some vortices, which are also visible in the experiments. There are two specific instants in which the main jet flow in the mixing chamber directly reaches the exit nozzle causing the flow to switch up the exit branch (from upper to lower or vice versa), producing mainly low negative pressures in the mixing chamber (Figure 7j,m).
Looking at the right column of Figure 7, it is noticed that pressure grows in the mixing chamber when the jet impinges on the exit nozzle walls; it is maximum when the jet impacts the exit nozzle wall with the steepest angle (e.g., Figure 7b,i) and minimum when the flow directly reaches the exit nozzle. The high pressure values at the impinging points are clearly noticed in Figure h,i,k,l, and the large vortex formed in the outlet branches in Figure 7h,l presents the minimum pressure values. In this way, an oscillation pattern of the pressure field is developed, which is connected with the mass flow cycle (in fact, they have opposite phases), a fact that has been referred to as "axial pumping" [14] in the fluid actuator. Interestingly, previous studies usually report only the behavior of the velocity field, so that the present work is one of the first studies to show explicitly the pressure field cycle inside this kind of device. Figure 8 displays the turbulent kinetic energy k contour field along one oscillation cycle at the corresponding times of Figure 7. It can be noticed that the highest values of k are attained in the mixing layers generated by the jet in the mixing chamber and in the outlet branches, as well as in the recirculation zones that flank the jet between the upper and lower blocks. It is precisely in such areas where large velocity gradients are generated, leading to substantial turbulent kinetic energy production. However, in the surroundings where the jet impinges the exit nozzle walls, k adopts moderate values because in such zones, turbulent dissipation is high. The large vortical structures developing in the mixing chamber and in the outlet branches are zones of high values of turbulent kinetic energy (see Figure 8a,c,d,f,g), as already observed in the study Guzmán et al. [3], which also exhibits low pressure values. In Figure 8b Next, the sensibility of the results regarding the turbulence model is presented for the base case. Apart from the mentioned SST k-ω turbulence model, the standard k-ε, Transition SST and Reynolds Stress Model (RSM) formulations have been tested at ≈ 16,000. The standard k-ε model [19] is by far the most used two-equation turbulence model due to its numerical robustness and stability; however, it has well-known limitations regarding the description of rotating and adverse pressure gradient flows with boundary layer separation. The RSM [20] is a seven-equation turbulence model which solves differential equations for the Reynolds Stress tensor components, avoiding the second moment closure by means of the Boussinesq approximation; it is claimed to be superior to classical two-equation models for describing anisotropic flows and flows with strong curvature and/or swirl effects and adverse pressure gradients. Finally, the transition version of the SST k-ω model [21] has been evaluated. This model includes a for- Next, the sensibility of the results regarding the turbulence model is presented for the base case. Apart from the mentioned SST k-ω turbulence model, the standard k-ε, Transition SST and Reynolds Stress Model (RSM) formulations have been tested at Re ≈ 16, 000. The standard k-ε model [19] is by far the most used two-equation turbulence model due to its numerical robustness and stability; however, it has well-known limitations regarding the description of rotating and adverse pressure gradient flows with boundary layer separation. The RSM [20] is a seven-equation turbulence model which solves differential equations for the Reynolds Stress tensor components, avoiding the second moment closure by means of the Boussinesq approximation; it is claimed to be superior to classical two-equation models for describing anisotropic flows and flows with strong curvature and/or swirl effects and adverse pressure gradients. Finally, the transition version of the SST k-ω model [21] has been evaluated. This model includes a formulation for taking into account the effects of the transition from laminar to the turbulent boundary layer. It solves two extra equations, one for the transition momentum thickness Reynolds number and another for the intermittency; therefore, it is really a four-equation turbulence model that has shown its potential in a wide range of engineering applications providing good agreement with experiments (e.g., [22][23][24][25]). Figure 9 shows the normalized mass flow profiles at outlet 1 in the base case at Re ≈ 16, 000. As can be seen in Figure 9a, there are some differences between the additional tested turbulence models, but all of them provide a shorter period than the SST k-ω model. This situation is similar to what happened in [15], where the other employed turbulence models (k-ω, k-ε, RNG k-ε and k-ε EARSM) also provide larger oscillation frequencies than the SST k-ω. In the present study, all the alternative models provide mass flow curves which are very similar in shape but showing some noticeable differences in their peak values. Moreover, the three investigated turbulence models present curves without inflection points in the ascending and descending branches, while the SST k-ω model shows two of such points in these branches (i.e., around times 0.02 s and 0.04 s). Although not shown, such particular points appear when no flow crosses one of the two outlets. As a comment, the value of the time step necessary to obtain a converged solution using the k-ε and RSM turbulence models is notably smaller than for the SST models (k-ω and transition).  Figure 9b shows the effect of modifying the turbulence intensity level at the inlet of the fluidic oscillator (SST k-ω turbulence model). As it can be noticed, only very marginal changes happen when it is increased from a very low level (0.05%) to a more than moderate level of 10%. The reason is connected with the fact that the flow spends some time being accelerated in the power (inlet) nozzle. During such transit the turbulence kinetic energy decays until the flow reaches the mixing chamber, where strong shear layers are generated in the sudden expansion (see Figure 8).

Effect of Feedback Channels Geometry Modification on Output Frequency
The operational parameters of the studied fluidic oscillator are governed by its geometric configuration, which means that its output frequency is fixed for a given incoming Reynolds number. Therefore, an opportunity to modify such frequency arises from the alteration of the internal geometry, which would allow a customization of the fluidic actuator for specific applications.
A previous study evaluated the effect of the width of the inlet and outlet sections of the mixing chamber, as well as the angles of the inner walls of the exit nozzle on the frequency and amplitude of the mass flow signal through the upper outlet [16]. They found a certain dependence of frequency and amplitude over the investigated parame-  Figure 9b shows the effect of modifying the turbulence intensity level at the inlet of the fluidic oscillator (SST k-ω turbulence model). As it can be noticed, only very marginal changes happen when it is increased from a very low level (0.05%) to a more than moderate level of 10%. The reason is connected with the fact that the flow spends some time being accelerated in the power (inlet) nozzle. During such transit the turbulence kinetic energy decays until the flow reaches the mixing chamber, where strong shear layers are generated in the sudden expansion (see Figure 8).

Effect of Feedback Channels Geometry Modification on Output Frequency
The operational parameters of the studied fluidic oscillator are governed by its geometric configuration, which means that its output frequency is fixed for a given incoming Reynolds number. Therefore, an opportunity to modify such frequency arises from the alteration of the internal geometry, which would allow a customization of the fluidic actuator for specific applications.
A previous study evaluated the effect of the width of the inlet and outlet sections of the mixing chamber, as well as the angles of the inner walls of the exit nozzle on the frequency and amplitude of the mass flow signal through the upper outlet [16]. They found a certain dependence of frequency and amplitude over the investigated parameters, as mentioned in the introduction. Here, however, we focus on the influence of the geometry of the feedback loops on the output operational frequency of the fluidic actuator. For that, two extra configurations have been considered, which are shown in Figure 10. The grids of the modified configurations are similar in size and quality to the intermediate mesh of the original fluidic oscillator. In the C1 case, it consists of about 328,000 cells, whereas in the C2 case, it is around 314,000 cells. The numerical set up employed for these new configurations was the same as in the original case described in Section 2, including the turbulence model. Figure 11 shows the normalized mass flow through outlet 1 versus time for the three geometrical configurations: original, C1 and C2 at ≈ 16,000. It can be seen that as the feedback channels get rounder, the frequency of the oscillation decreases (i.e., the period is larger) and the mass flow curve is flattened in some parts of the increasing and decreasing branches (see areas inside the ovals in Figure  11) containing the inflection points; even in the C2 case, a plateau seems to be formed. On the other hand, the minimum and maximum values of the normalized mass flow in the three curves are very similar.  Figure 12 shows velocity contours in the three geometrical configurations at approximately the same moment in the oscillation cycle, when the jet is bended downwards and the fluid is mainly divided between the exit nozzle of the mixing chamber and the As can be appreciated from Figure 10, the feedback loops connecting the rear area of the mixing chamber and its front part are rounder than in the original configuration, i.e., the radius of curvature of the channel elbows is larger. Moreover, the corners of both blocks, upper and lower, in the new configurations are not sharp 90 • angles but have a certain curvature radius to preserve the width of the channels, and even in the case of C2, two of the corners have disappeared in each feedback loop. Such geometrical changes are aimed at easing the fluid circulation along the feedback loops.
The grids of the modified configurations are similar in size and quality to the intermediate mesh of the original fluidic oscillator. In the C1 case, it consists of about 328,000 cells, whereas in the C2 case, it is around 314,000 cells. The numerical set up employed for these new configurations was the same as in the original case described in Section 2, including the turbulence model. Figure 11 shows the normalized mass flow through outlet 1 versus time for the three geometrical configurations: original, C1 and C2 at Re ≈ 16, 000. It can be seen that as the feedback channels get rounder, the frequency of the oscillation decreases (i.e., the period is larger) and the mass flow curve is flattened in some parts of the increasing and decreasing branches (see areas inside the ovals in Figure 11) containing the inflection points; even in the C2 case, a plateau seems to be formed. On the other hand, the minimum and maximum values of the normalized mass flow in the three curves are very similar. Figure 12 shows velocity contours in the three geometrical configurations at approximately the same moment in the oscillation cycle, when the jet is bended downwards and the fluid is mainly divided between the exit nozzle of the mixing chamber and the upper feedback channel. Very similar flow patterns for the three geometries are observed at the outlet branches, with a large vortex in the lower branch and with very similar intensity. However, differences can be noticed in the upper recirculation loop, where the original (Figure 12a) and C1 (Figure 12b) configurations clearly show separated flow at the inlet corner of the feedback channel, but the fluid in the C2 configuration (Figure 12c) attaches to the outer wall of the loop without presenting any separated vortices. Moreover, the velocity magnitude within such vortices is lower in the C1 geometry, its shape being more elongated than in the original case; in other words, the vortex is stretched to conform to the rounded walls. Another vortex is visible in the original configuration, close to the lower right corner of the upper block and near the entrance of the upper feedback channel; such a vortex is not so evident in our configurations, possibly due to the rounded corners.
328,000 cells, whereas in the C2 case, it is around 314,000 cells. The numerical set up employed for these new configurations was the same as in the original case described in Section 2, including the turbulence model. Figure 11 shows the normalized mass flow through outlet 1 versus time for the three geometrical configurations: original, C1 and C2 at ≈ 16,000. It can be seen that as the feedback channels get rounder, the frequency of the oscillation decreases (i.e., the period is larger) and the mass flow curve is flattened in some parts of the increasing and decreasing branches (see areas inside the ovals in Figure  11) containing the inflection points; even in the C2 case, a plateau seems to be formed. On the other hand, the minimum and maximum values of the normalized mass flow in the three curves are very similar.   Therefore, as Figure 12 shows, the fluid motion inside the feedback channels is substantially altered by their geometry, which is reflected by a change in the operating frequency. As a result, such a frequency is reduced as the feedback loops get rounder.

Preliminary Three-Dimensional Results
In this section, some preliminary results of the three-dimensional simulation are shown at ≈ 16,000. Figure 13 shows the structured grid, which consists of 2.85 million cells, a slightly refined mesh than in previous studies [15,16]. The simulation methodology and the numerical schemes used were the same as in the 2D computation. The chosen turbulence model was SST k-ω, because previously it was found that it provided the closest results to the experimental measurements. Additionally, the time step was set to ∆ = 10 s, which is one order of magnitude smaller than in the two-dimensional case. As a result, the total computational time in 3D is nearly two orders of magnitude higher than in 2D. With the previous settings, the average value Therefore, as Figure 12 shows, the fluid motion inside the feedback channels is substantially altered by their geometry, which is reflected by a change in the operating frequency. As a result, such a frequency is reduced as the feedback loops get rounder.

Preliminary Three-Dimensional Results
In this section, some preliminary results of the three-dimensional simulation are shown at Re ≈ 16, 000. Figure 13 shows the structured grid, which consists of 2.85 million cells, a slightly refined mesh than in previous studies [15,16]. Therefore, as Figure 12 shows, the fluid motion inside the feedback channels is substantially altered by their geometry, which is reflected by a change in the operating frequency. As a result, such a frequency is reduced as the feedback loops get rounder.

Preliminary Three-Dimensional Results
In this section, some preliminary results of the three-dimensional simulation are shown at ≈ 16,000. Figure 13 shows the structured grid, which consists of 2.85 million cells, a slightly refined mesh than in previous studies [15,16]. The simulation methodology and the numerical schemes used were the same as in the 2D computation. The chosen turbulence model was SST k-ω, because previously it was found that it provided the closest results to the experimental measurements. Additionally, the time step was set to ∆ = 10 s, which is one order of magnitude smaller than in the two-dimensional case. As a result, the total computational time in 3D is nearly two orders of magnitude higher than in 2D. With the previous settings, the average value The simulation methodology and the numerical schemes used were the same as in the 2D computation. The chosen turbulence model was SST k-ω, because previously it was found that it provided the closest results to the experimental measurements. Additionally, the time step was set to ∆t = 10 −6 s, which is one order of magnitude smaller than in the two-dimensional case. As a result, the total computational time in 3D is nearly two orders of magnitude higher than in 2D. With the previous settings, the average value of y + was around 2, still an acceptable value that allows resolving the boundary layer, and the average Courant number was lower than 0.1 (maximum 0.9), which guarantees both the stability and accuracy of the numerical method.
The normalized mass flow through outlet 1 as a function of time is shown in Figure 14. As a result, the computed operating frequency in the 3D case was somewhat higher than in 2D, similar to the findings of Krüger et al. [15]. It can be seen that the peaks of mass flow are smaller in the three-dimensional simulation than in the two-dimensional one, a fact probably related to the strong secondary flows observed in the outlet branches of the former, which are out of reach for the latter. Moreover, the position of the minimum of the curve is noticeably displaced in time in both simulations; however, the location of the maximum is quite similar in both curves. Another fact to be noticed is that the 3D simulation curve does not show the inflection points previously seen in the 2D computation.
Appl. Sci. 2022, 12, x FOR PEER REVIEW 14 of 17 mass flow are smaller in the three-dimensional simulation than in the two-dimensional one, a fact probably related to the strong secondary flows observed in the outlet branches of the former, which are out of reach for the latter. Moreover, the position of the minimum of the curve is noticeably displaced in time in both simulations; however, the location of the maximum is quite similar in both curves. Another fact to be noticed is that the 3D simulation curve does not show the inflection points previously seen in the 2D computation. The three-dimensional simulation provides many more details about the flow structure inside the fluidic oscillator than 2D computations. In particular, the flow at the entrance of the feedback channels shows the typical counter rotating vortices, known as Dean vortices, that appear in pipes and channels when the flow moves along a 90° elbow [26,27]. In such a configuration, the progress of the primary flow is strongly influenced by the generation of a secondary flow in the cross section of the duct, which results from a balancing among the existing pressure gradient between the inner and outer walls and the centrifugal acceleration experienced by the fluid. This is a typical 3D effect that is naturally absent in 2D and that intensifies the pressure loss of the flow within the recirculation loops. Just as an illustration, Figure 15 shows the vorticity in the stream-wise flow direction in the mid plane of the fluidic oscillator (Figure 15a) as well as some of the vortices present in the flow, represented by the grey structures (Figure 15b). The three-dimensional simulation provides many more details about the flow structure inside the fluidic oscillator than 2D computations. In particular, the flow at the entrance of the feedback channels shows the typical counter rotating vortices, known as Dean vortices, that appear in pipes and channels when the flow moves along a 90 • elbow [26,27]. In such a configuration, the progress of the primary flow is strongly influenced by the generation of a secondary flow in the cross section of the duct, which results from a balancing among the existing pressure gradient between the inner and outer walls and the centrifugal acceleration experienced by the fluid. This is a typical 3D effect that is naturally absent in 2D and that intensifies the pressure loss of the flow within the recirculation loops. Just as an illustration, Figure 15 shows the vorticity in the stream-wise flow direction in the mid plane of the fluidic oscillator (Figure 15a) as well as some of the vortices present in the flow, represented by the grey structures (Figure 15b). Figure 15a shows the vorticity in the flow direction (x-axis) in the mid plane of the device. Moreover, some vertical planes cutting the feedback loops and the mixing chamber are included, which are also colored by the x-component of the vorticity, ω x . In them, some streamlines of the secondary flow are drawn, illustrating the complex three-dimensional flow structure in both the mixing chamber and the recirculation channels. The highest absolute values of ω x appear in the mixing layer formed when the jet enters the mixing chamber, the strong recirculation zones in the mixing chamber, the entrance in the feedback channels, and upstream and downstream of the exit nozzle, including the corresponding outlet branch. a balancing among the existing pressure gradient between the inner and outer walls and the centrifugal acceleration experienced by the fluid. This is a typical 3D effect that is naturally absent in 2D and that intensifies the pressure loss of the flow within the recirculation loops. Just as an illustration, Figure 15 shows the vorticity in the stream-wise flow direction in the mid plane of the fluidic oscillator (Figure 15a) as well as some of the vortices present in the flow, represented by the grey structures (Figure 15b).  Figure 15b presents the fluid velocity contour plot in the oscillator mid plane; in it, the curved jet impinges on the upper wall of the exit nozzle, dividing the flow between the lower outlet branch and the upper recirculation channel. Vortical structures are also shown in the same plot as grey surfaces; they have been identified with the help of the Q criterion [28] for the value Q = 2 × 10 5 s −2 . In that figure, a long vortex perpendicular to the impinging jet (A) is observed, which extends towards the entrance of the upper feedback channel and crosses the exit nozzle (whose fingerprint is also seen in the vertical plane at the exit nozzle, Figure 15a). Another similar, counter-rotating vortex is formed on the other side of the device, hindered by the midplane. Such vortices arise from the interaction of the flow emerging from the impingement point with the top and bottom walls (orthogonal to the z-axis) that determine the width of the device. Moreover, vortical structures are clearly visible at the inlets of the recirculation loops (B), owing to the aforementioned Dean vortices, in the large recirculation area inside the mixing chamber (C), above the jet, in the two outlet branches (D) and in the separated flow zone behind the second elbow of the upper recirculation loop (E). Altogether, such flow features evidence the complexity of the turbulence dynamics inside the fluidic oscillator.
It has to be said that the graphs in Figure 15 are a representation of the flow at a certain instant of time and that such flow structures will evolve along the oscillation cycle. Therefore, future work will be devoted to elucidating the different stages of the vorticity and coherent structures development in the fluidic actuator along the complete cycle.

Conclusions and Future Work
A numerical study of the flow inside of an existing microfluidic oscillator has been carried out in this contribution. Two-dimensional simulations have been implemented to determine the oscillation frequency of such a device, which has been validated against experimental measurements and previous 2D and 3D computations in a range of Reynolds numbers. The computed flow field compares favorably with the experimentally measured velocity fields obtained with particle image velocimetry (PIV). Additionally, in order to improve the understanding of the flow dynamics inside the actuator, the pressure and turbulent kinetic energy fields along the oscillation cycle have been analyzed; such fields are not usually shown in the related literature. The main conclusion is that the 2D computations suffice to predict the operating oscillation frequency of the fluidic actuator. Moreover, it was found that the normalized mass flow curve in the oscillator versus non-dimensional time is independent of Re for the considered range; therefore, the non-dimensional frequency is also independent of the Reynolds number. Further, the effect of the modification of the feedback channels' geometry on the output frequency of the oscillator has been studied; it was concluded that the rounder designs of such feedback loops implied a reduction in frequency, a fact connected with the reduction in the number of vortices appearing in the flow field. Finally, preliminary results of a full three-dimensional simulation are presented, which allows appreciating the development of coherent structures and the vorticity field, characterizing the complex turbulent flow inside the fluidic oscillator.
In the near future, it is planned to perform a thorough analysis of the behavior of the three-dimensional vorticity and coherent structures along the complete cycle of oscillation of the fluidic actuator. Such a study would shed light on the complex flow inside this kind of device that will improve their performance in specific applications.