HYSTERESIS OF PLANAR DOUBLE SLOT IMPINGING AIR JETS

In present study, the bistability and hysteresis of the non-isothermal flow behind two slot nozzle is numerically investigated. Theoretical review of bistability of isothermal stagnation flow is presented and both flow field patterns that exist in the region of bistability are commented on. Hysteresis is found in number of conducted simulations and effects of various parameter changes are discussed. Moreover, the mechanism of transition from one flow field pattern to another is discussed for both transition directions. In conclusion, validity of 2D simulation results against real 3D problems is of concern and general contribution of this research is discussed. Throughout this work, emphasis is on applications regarding heat transfer.


Introduction
In many engineering applications, heat and mass transfer are of great importance. Due to the perpendicular orientation of the flow towards the wall, impinging jets achieve the highest heat transfer rates ̇ [ 2 ] among all the convective processes. [1] The heat transfer rate is given by Newton law of cooling which separates the contribution of temperature difference of the heated surface and the fluid (at sufficiently large distance from the surface) and all other factors represented by heat transfer coefficient [ 2 ]. The equation therefore is ̇= = ( − ∞ ).
Compared to flow parallel to the cooled surface, jet impingement produces heat transfer coefficients that are up to three times higher at a given maximum flow speed. [2] In the case of cooling of electronic devices, it is common requirement to keep the temperature of cooled component below given value. Since the temperature of cooling substance is, in case of an air jet, usually given by external conditions and can be therefore quite high, it is desirable to employ cooling process that produces as high heat transfer coefficient as possible.
An impinging jet is configured by high momentum flow discharged through a nozzle (or series of nozzles) and its subsequent impingement on the target surface.
The greatest drawback in the endeavor to achieve highest possible heat transfer coefficient is the existence of so-called viscous sublayer, which forms on the impingement wall. Fluid in this sublayer is effectively stationary and through its thickness, the heat is therefore transferred only by conduction governed by Fourier law ̇= − , (2) where [ ] is the coefficient of thermal conductivity. Due to favourable orientation of the flow relative to the target surface, perpendicular impinging jets achieve lowest possible viscous sublayer thickness.
Impinging jets have therefore large potential in heat and mass transfer applications and are subject of extensive research. Over past years, many reviews were published, [1,3,4,5,6] to name a few. This paper deals with double slot impinging jets and aims to shed some light on hysteresis phenomenon which this configuration exhibits. Hysteresis describes the response to a quasi-steady change of one of the parameters that influence the flow field. If the response of a system depends on the direction of parameter change, hysteresis is present and the relationship between specific parameters forms the hysteresis loop. Bistability, which is closely related to hysteresis, occurs when boundary conditions alone can not determine the flow field uniquely. This leads to the existence of two stable steady states.
Hysteresis in double slot isothermal impinging jets was already observed by Travnicek and Krizek [7] and some theoretical background was given by Marsik et al. [8]. Shariatmadar et al. presented a study concerning heat transfer characteristics of laminar slot jet arrays impinging on a constant target surface temperature. [9] In their paper, they also investigated double slot impinging jet configuration and they obtained results corresponding to both flow patterns of concern in this paper. However, they did not comment on hysteresis nor bistability.
Experimental research of hysteresis in impinging jets is recently focused on annular jets, which exhibit similar flow patterns as the double slot jets (main difference being that the vortices in double slot impinging jets are in fact vortex rings in annular impinging jets).
Bistability and hysteresis on annular impinging jets was theoretically predicted in 1991 by Kokoshima et al. [10].
Travnicek et al. [11] used fluidic control (the term fluidics means a technology of flow handling based on the fluid flow interaction without action of moving components) and achieved two stable flow field patterns of annular impinging jet both experimentally and numericaly. Moreover, hysteresis was observed in their experiments.
Furthermore, Travnicek and Tesar found two different regimes of the time-mean flow field, differing in the size of the central recirculation zone, with either the single stagnation point or the stagnation circle. They used acoustic excitation to switch between the two regimes. [12] In another publication, Travnicek and Tesar experimentaly investigated annular impinging air jet [13]. In their work, the bistability and hysteresis were found, but only at small Reynolds numbers, and at small annular nozzle slot widths. Also, they concluded that an extension of the nozzle center-body promotes the tendency to the hysteretic behaviour.
On the other, hysteresis in annular impinging jets is not universal phenomenon as is suggestested by work of Tesar et al. [19]. They reported on five flow regimes of annular jet, however, no hysteresis was observed.
While there is noticeble effort to study hysteresis in impinging jets experimentaly, there is significant lack of data from numerical simulations. Similar approach as in this paper was used by one of the authors to model hysteresis in annular impinging isothermal flow [15]. Other than that, no contribution employing computational fluid dynamics to investigate hysteresis of impinging jets was found in available literature.
In present study, hysteresis is evaluated via nozzle moving with relatively low velocity relative to the target surface. Therefore the hysteresis was evaluated only to some extent. That is because hysteresis is defined as response to the directional change of some parameter alone. Here, the movement of the nozzle has also some influence.
We attempted to control the flow structure be means of motion of the nozzle. This approach is, off course, not suitable for use in applications in heat transfer where some hypothetical device would need to function efficiently continuously for relatively long period of time. However, it gives some insight to the flow field pattern change.
Effects connected with hysteresis might become highly important in microelectronics to cool thermally highly loaded components. Length scales of impinging jets suitable for applications in microelectronics are generally too low for turbulence to take place. Jets are therefore effectively laminar, heat transfer rates are than not enhanced by turbulence and need rises for other mechanisms to compensate for it.

Single slot impinging jet
Single slot impinging jet is, in this case, composed of single high aspect ratio rectangular nozzle and a target plate below it. Here, the flow is modelled as 2D with the gradients along longer side of the nozzle profile being omitted.

Transition to turbulence
It was found by Lemanov et al. [16] that submerged air jet exiting from plane channel of width 600 exhibits transition to turbulence roughly at the distance of 30 from the exit. Since the slot width in our simulations was 0.5 and the target surface was not further than 14 from nozzle exit, it is justified to assume that the flow is laminar.

Double slot impinging jets
Double slot impinging jets are formed by two parallel single impinging jets in proximity to each other.

Two stable flow field patterns
The existence of two flow patterns in double slot impinging jet was already investigated by Marsik et al. [8]. In their paper, theoretical guidelines for the bistability of double slot nozzle jet were given. They even provided some experimental results, however they did not succeed in reproducing the results by CFD simulation.
Double slot impinging jets exhibit two stable flow field patterns. Existence of specific pattern is influenced mostly by the nozzle to wall distance. Other parameters, such as distance between the two slots, width of the slot, flow rate through the nozzles and others play significant role in influencing the hysteresis itself. The two main flow field patterns are illustrated in figure 1.

Fig. 1. Flow field patterns of type A and B [8]
Here, the two main flow field patterns of the annular impinging jets will be referred to as A and B. The main characteristic of the pattern A is the existence of the reattachment point R. Below the point R, the flow field has a similar character to the flow field behind a single slot nozzle. If flow pattern A occurs, the region upstream of point R will be called the intermediate merging zone and the region downstream of point R will be called the fully merged zone. Flow pattern A is generally found for larger nozzle-to-wall distances.
Flow field pattern B is, on the other hand, characterized by the tendency of the flow to bend away from the nozzle axis and is generally found for lower nozzle-to-wall distances.

Hysteresis Evaluation
In present paper, the hysteresis of double slot impinging jet is evaluated using the computational fluid dynamics (CFD) simulation. This approach was chosen mainly because it provides information about the entire flow field. Moreover, numerical simulations have the advantage over experiment in investigating wide range of geometry modifications and other input parameters easily.

Mathematical model
Hysteretic behavior of annular impinging jets was investigated numerically in Fluent using the finite volume approach for solving appropriate conservation equations of fluid mechanics. Hysteretic behavior was confirmed through the dynamic mesh simulation where the nozzle moved relatively to the target surface with the velocity magnitude significantly smaller than the magnitude of the flow nozzle exit velocity.

Governing equations and boundary conditions
Mass, momentum and energy conservation equations together with the equation of state of ideal gas form a closed set of 6 coupled scalar equations for 6 unknown variables. The equations are [17] ∂ρ ∂t + ∇. (ρv ⃗ ) = 0. ( Although the dimensions and air density are small in value, gravity term in momentum balance must be included in order to account for buoyancy forces acting on air at different temperatures. Equations (3) to (8) together with appropriate boundary conditions give a solvable problem.
The velocity inlet profile is given by the analytical solution to a steady 2D flow between two infinite planes driven by pressure difference in one direction [18], i.e.
where x is distance from the jet axis. The ratio of pressure difference to the slot length is the source term in equation (9). Although it is custom to give fluid inlet velocity constant profile [1], a fully developed channel profile was used here. This is because long channel of the nozzle can help the flow achieve steady state without any disturbances which would influence the flow downstream from the nozzle and which would be therefore an unknown effect in future experimental study.

Computing
As was already mentioned, the flow was modelled as two dimensional. Also, in respect of low velocities and low dimensions of the investigated flow, simple model of laminar flow was used. Furthermore, the fluid thermophysical properties (i.e. viscosity and thermal conductivity) were assumed constant. This assumption is justified if there the temperature differences in the domain are kept small. Pressure-based solver was used to calculate the fields of variables that were solved for. Coupling between continuity and momentum equations were solved for iteratively by means of the PISO algorithm [19]. Second order discretization schemes are used for pressure and second order upwind schemes are used for solution of momentum and energy equations. Whether convergence has been reached was evaluated based on scaled residual evolutions. Limits for scaled residuals were set to 10 −3 for continuity and momentum equations and to 10 −6 in case of energy equation.

Hysteresis evaluation method
Several configurations were investigated. For future reference, those will be denoted Jet1 to Jet8. Table 1 gives details of individual configurations in terms of dimension and table 2 gives main characteristics of the flow.
For the purpose of characterization of the flow, Reynolds number is defined with b as the characteristic dimension and with the average nozzle exit velocity.  (10), the Reynolds number for all the cases is 50.

Re
The dynamic mesh case with layering was simulated, there have been sporadic step-wise deviations from the smooth trend of investigated variables. To produce a nice plots, data was smoothed by fitting successive sub-sets of adjacent data points with low-degree polynomial. This method is known as the Savitzky-Golay filter. [20]

Results
In all following graphs, vertical lines of corresponding colour denote the point of greatest heat transfer coefficient for given configuration. Furthermore, all the graphs are non-dimensionalized.
It should be noted that in the non-dimensionalized plots is different for Jet6 and Jet7 as compared to the rest of the cases.

Streamlines of flow field patterns A and B
In this section, streamlines for both patterns as well as for chosen state of transition between individual patterns is shown.  Figure 2c) then show some state of transition between pattern B and pattern A (meaning that the nozzle now moves away from the surface). The transition occurs by splitting the two main vorticle structures of pattern B into four. Transition to pattern A is complete once the vertical structures next to target surface vanish completely. Similarly, the transition from pattern A to pattern B is again accompanied by splitting the two existing vortices into four.
The difference between the two transitions is that in case of transition from pattern B to pattern A, the two vortices are stretched in the direction perpendicular to the wall until they become instable and are then split into four roughly separate structures, two of which then decrease in size. On the other hand, in the case of transition from pattern A to pattern B, the two main vortical structures reach close enough to the target surface, the main flow does not reach to the intersection of double slot nozzle axis with the target surface and the pattern B is formed. No new vortex pair next to the target surface is formed in this case.
Above reasoning suggests that the bistability region is in fact occupied by patterns similar in nature to those shown in figure 2c and 2d rather than the two basic pattern A and pattern B.   shows configurations that only differ in relative velocity between the two slot nozzle and the target surface, therefore providing means to determine the effects of change in relative velocity. As the nozzle velocity gets lower, the flow field has more time to adapt to new circumstances. As a result, the hysteresis loop gets narrower. Figure 2b shows the effect on nozzle center-body extention and insertion on the hysteresis. It is clear that nozzle center-body extention shifts the hysteresis loop towards higher nozzle-to-wall distances. This is because the two main vortex structures are left with less space than in the case of Jet1. Interestingly, the hysteresis loop is shifted of its original position by a distance that is approximately equal to four times the extention value. The insertion makes the transition from pattern B to pattern A occur earlier. This is because the two vortical structures of pattern B are elongated thanks to the space made free by the insertion.

Stagnation pressure at the center of target surface
Finally, figure 2c demonstrates the effect of the double slot nozzle spacing. Naturaly, hysteresis region shifts towards lower nozzle-to-wall distances for lower double slot nozzle spacing. The central stagnation pressure for case of pattern A is greater as the slots are closer to each other. This is because the vortical structure below the center-body are smaller and they therefore require less energy from the main flow.
Whether any of the concidered configurations of the planar double slit jet exhibits bistability remains open question. However, some conclusions can be made for configuration Jet1. In this case, as the velocity of the nozzle approaches zero (the relative nozzle velocity is reduced by half and by quarter in configurations Jet2 and Jet3 respectively), the region of hysteresis loop that is shown in figure 3a should decrease to the region of bistability. It is believed by the authors that the region of bistability would be achieved when the slope of the two transition lines would become 2 . Based on the character of the slope increase, it can be concluded that the region of bistability does in fact exist for this particular configuration.
Courtesy of the famous Courant-Friedrichs-Lewy condition [21], simulation takes longer and longer for lower nozzle velocities and investigating flow bistability through this approach is therefore not possible.

Maximum stagnation pressure
From the point of view of practical applications, much more interesting parameter is the maximum stagnation pressure. Figure 3 shows maximum stagnation pressure as a function of nozzle-to-wall distance. Since these maxima are not at constant distance from the jet axis, hysteresis loops such as the ones in figure 4 do not not emerge.
Here, discussion will be focused on local maxima obtained in transition from pattern A to pattern B.   Figure 4 shows that regions of maximum stagnation pressure are also regions of maximum heat transfer coefficient (as was previously mentioned, vertical lines indicate corresponding maxima in heat transfer coefficient). Figure 3a shows that, for larger nozzle velocities, higher maximum stagnation pressure is achieved. To illustrate the significance of this increase, table was set to calculate relative changes of important parameters between successive configurations. Since configuration Jet3 has no predecessor, nothing will be computed. About 1.7% increase of maximum stagnation pressure was found from case Jet1 to case Jet2 and about 2.5% increase of maximum stagnation pressure was found from case Jet2 to case Jet3.
Effect of nozzle center-body extention or insertion is not beneficial with respect to maximum stagnation pressure. But still, nozzle center-body extentions might be useful in cases where there is need to shift the maximum stagnation pressure position towards higher nozzle-to-wall distances. Figure 3c shows that the maximum of stagnation pressure increases with decreasing double slot spacing and, at the same time, also shifts towards lower nozzle-towall distances.

Heat transfer
Finally, maximum heat transfer coefficient, nondimensionalized to the form of Nusselt number, is plotted as a function of the nozzle-to-wall distance.
In order to determine the heat transfer coefficient, it had to be calculated from the total surface heat flux.
Goldstein et al. [22] carried out experimental tests and found that the convective heat transfer coefficient is independent of the temperature difference between jet and ambient as long as it is defined with the difference between the heated wall temperature and the adiabatic wall temperature. This gives the heat transfer coefficient as where q̇ is the total surface heat flux from the target surface.   Figure 5 shows the maximum local Nusselt number as a function of non-dimensionalized distance of the target surface from the nozzle exit. In respect of investigated configurations, overall best heat transfer rate is obtained for larger relative velocity of the nozzle, for the case where no insertion or extention of nozzle center-body is employed and for the case where the two slots are as close as possible. For all the cases, the increase is quite significant. Both with respect of comparison between individual cases and with respect to the heat transfer coefficient for other surface to nozzle distances.

Positions of maximum stagnation pressure and maximum heat transfer coefficient as a function of nozzle-to-wall distance for transition region from pattern A to pattern B
From the viewpoint of cooling applications, the most promising seems to be the region transition from pattern A to pattern B. This region has, as is shown in figure 4, local maximum of the stagnation pressure and consequently also local maximum of heat transfer coefficient.  Vertical lines show that the overall maxima in heat transfer coefficient coincide with the overall maxima for stagnation pressure. Interesting is that generally, in this transition regime, maxima of heat transfer coefficient for other nozzle-to-wall distances are not located at the same place as the stagnation pressure maxima. Maxima in heat transfer coefficient, with the exception of the overall maximum, tend to be further away from the jet axis than the stagnation pressure maxima. This suggests that the effect of overall maxima in heat transfer coefficient is achieved by some combination of impinging jet and wall jet.
As discussed above, the overall maxima in heat transfer coefficient can be shifted further away from the jet axis either by addition of nozzle center-body extention or by increasing the nozzle slot spacing. It is interesting to note that in case of lower nozzle velocities, the overall maximum is also shifted away from the axis. However, this effect is not as significant as the other two mentioned above.

Stability and Validity of 2D model
For lower nozzle exit velocities, investigated flows were all stable and results from corresponding simulations were provided in previous text. However, for larger nozzle exit velocities, flow became unstable.
Generally, flow might be instable when there are destabilizing factors present. Those factors are namely the presence of shear region within viscous fluid and adverse pressure gradient. [18] Instable flow was obtained for the case when the domain height was diminishing, while, for the same nozzle exit velocity, the flow was stable when the domain height was stretching. This is because the relative movement of the nozzle toward the target surface is a source of the adverse pressure gradient, which is, as mentioned, one of the destabilizing factors.
The other destabilizing factor is the shear region within viscous fluid and it is present in all flows of real fluids, double slot jets not excluded. Therefore there is also a Reynolds number threshold at which, flow with nozzle moving away from the target surface becomes instable.

Assumption of 2D flow
Real 3D double slot jet with slot lengths much larger than slot widths might be mathematically modeled as 2D flow, taking into account only single section, supposing the same flow picture for any value of the invariant parameter (i.e. the position along the length of the slot). Such assumption is justified everywhere except at the ends of the slots.
Uruba shows [15] that even for low Reynolds numbers, omitting 3D structure of the flow could lead to misinterpretation of the flow field and wrong indentification of physical mechanisms involved in flow process. This was demonstrated for the case of a wake behind a cylinder which is geometrically similar to the double slot jet case presented here (it can be considered a nozzle centerbody in a cross-flow). However, the results might be correct on average.

Conclusions
Present paper deals with the hysteresis of 2D double slot jets. Hysteresis was confirmed for range of configurations and some quantitative understanding, applicable to real word applications, is obtained. Moreover, the process of transition of one flow field to another was describes and some interesting observations were made. To the best knowledge of the authors, no article on numerical evaluation of hysteresis in double slot impinging jets was yet published in available literature.
Provided results do not aim to give a complex view on the problem of hysteresis of planar double slot jets. However, some insight was provided to the effect of nozzle extention/insertion and slot spacing and possible means to identify the region of bistability were proposed.
Motivation for research concerned with bistability and hysteresis of jets lies in the development of more efficient systems for heat and mass transfer technologies. Presented results suggest a way to significantly increase local heat transfer coefficient. Such result is extremely important for the field of microfluidics where turbulence is generally not accessible to improve heat transfer.
The aim to achieve highest possible heat transfer rate is generally not the only consideration. From construction point of view, location of the heat transfer maximum as well as the effect of nozzle to wall distances are of concern. Insight for both these concerns and more are given by results provided in this paper.
Greatest drawback of present work is the fact that presented results are not numerical grid resolution independent. In the hysteresis region, one configuration is generally stable while the other is instable. The transition from one configuration to the other can therefore be caused by both the physical effect as well as numerical effect relevant to the actual solution of governing equations. The main focus regarding future work will therefore be in arriving at the numerical grid independent solution. In light of this paragraph, it is further emphasized that the value of this article lies mainly at the application of method of dynamic mesh for investigating hysteresis of double slot impinging jets and in comparison of various simulation results which were obtained under similar numerical conditions.
As for 2D flow, the instability thresholds with respect to at least some parameters are still to be investigated for the case of stretching domain as well as for the case of diminishing domain. Moreover, effect of different nozzle exit velocities is still to be examined. This is made difficult by the strong dependence of the resulting hysteresis loop on the nozzle velocity. In this paper, nozzle velocity was referred to as some small fraction of the average nozzle exit velocity. However, there is no guarantee that the effect, even approximately and for low nozzle velocities, depends linearly on the average nozzle exit velocity. Therefore, in order to really understand the effect of nozzle exit velocity on hysteresis, it is appropriate to investigate the bistability, which corresponds to the limit of velocity of the nozzle approaching zero. One of the topics of future research therefore should be to propose a method to