Abstract
This work focuses on the evolution of a free plane laminar jet in the near-nozzle region. The jet is buoyant because it is driven by a continuous addition of both buoyancy and momentum at the source. Buoyancy is given by a temperature difference between the jet and the environment. To study the jet evolution, numerical simulations were performed for two Richardson numbers: the one corresponding to a temperature difference slightly near the validity of the Boussinesq approximation and the other one corresponding to a higher temperature difference. For this purpose, a time dependent numerical model is used to solve the fully dimensional Navier-Stokes equations. Density variations are given by the ideal gas law and flow properties as dynamic viscosity and thermal conductivity are considered nonconstant. Particular attention was paid to the implementation of the boundary conditions to ensure jet stability and flow rates control. The numerical simulations were also reproduced by using the Boussinesq approximation to find out more about its pertinence for this kind of flows. Finally, a stability diagram is also obtained to identify the onset of the unsteady state in the near-nozzle region by varying control parameters of momentum and buoyancy. It is found that, at the onset of the unsteady state, momentum effects decrease almost linearly when buoyancy effects increase.
1. Introduction
Jet flow occurs in a variety of industrial applications such as pulverization, thermal isolation, film cooling, solid straightening, welding in aerodynamic, and hydrodynamic fields. For these applications, there exist two kinds of jet flow configurations: the axisymmetric jet, characterised by a circular nozzle, and the plane jet, characterised by a rectangular nozzle. Such types of jet flows suddenly are three-dimensional and show symmetric conditions that are object of different studies. The study of a two-dimensional jet flow configuration generally is one of the stages to better understand some phenomena that govern three-dimensional jet flows. In terms of numerical modelling, three-dimensional jet flows clearly are extremely more complex and demand long computing time.
Some works on laminar jets are limited to analytical or numerical solutions of the governing equations in the developed region. The solutions use a change of variables with the aim of ignoring the discharge conditions at the nozzle orifice [1]. For nonisothermal free laminar jets, works have been devoted to study the plume away from its source [1–5] and in the initial near-nozzle region [1, 6].
In reality, jet flows are turbulent and have been analysed both experimentally [7–9] and numerically [10–12]. Other experimental studies have analysed turbulent jet flow with variable density by considering two mixed fluids [13]. An analytical solution for the transition region, with buoyancy forces and inertial forces of same order of magnitude, has been addressed by Yu et al. [14] in laminar regime. These authors have introduced new parameters and used a change of variables. The equations they obtained were solved by taking into account the boundary conditions and the two integration constants derived from the momentum conservation and the heat flux. These two constants replace the discharge conditions at the nozzle orifice and were verified by Mhiri et al. [15].
Despite the number of studies on plane jet flow, its complexity still makes the necessity of studying the influence of various physical parameters related to the jet evolution. When taking into account density variations, the problem becomes even more complex due to the interaction of buoyancy-momentum. Therefore, the jet flow with variable density is still a large open subject and this work focuses on the study of the evolution of such a flow in the near-nozzle region. The jet flow configuration studied here is a free vertical plane jet, discharged from a rectangular nozzle, where the aspect ratio width/thickness is considered significant. For a nozzle of width and thickness , the lateral expansion of the jet can be neglected when the aspect ratio is . Thus, the flow can then be considered two-dimensional on average [16]: the average flow values in a given plane are identical in all other parallel planes. The present jet flow configuration is a kind of embedded or immersed jet, as it emerges in an environment constituted of the same fluid but colder. The jet flow is also free, since the boundary conditions theoretically are considered to be at the infinity.
In this study, the fully dimensional Navier-Stokes equations are solved and no simplifying assumption is considered on the role of pressure. The numerical simulations were also reproduced by using the Boussinesq approximation to find out more about its pertinence for this kind of flows. Different authors have focused on the criteria to validate the Boussinesq approximation. One of the most cited works about this subject is the work of Gray and Giorgini [17], who have shown that the Boussinesq approximation is well adapted to natural convection flows under small temperature gradients and negligible compressibility effects. For the case of air, regarding density variations with temperature, Gray and Giorgini have shown that the Boussinesq approximation is true for temperature differences approximately less than 28.6 K. More recently, other authors have carried out some studies with the aim of defining the limitations of the Boussinesq approximation [18, 19].
A time dependent two-dimensional numerical model is used in the present work, where density variations are given by the ideal gas law and flow properties as dynamic viscosity and thermal conductivity are considered nonconstant. The numerical approach is based on a second-order finite difference formulation and the coupled set of equations is solved by using an iterative predictor-corrector procedure at each time step. Some numerical tests were also carried out to define the appropriate boundary conditions since stability and flow rates control must be ensured.
2. Numerical Model
2.1. Governing Equations
The governing equations for conservation of mass, momentum, and energy, in Cartesian coordinates, are used in the following form: where is the velocity vector in the direction (in m/s), is the pressure (in Pa), is the density (in kg/m3), is the temperature (in K), is the specific gas constant (in J/(kg K)), is the specific heat at constant volume (in J/(kg K)), is the specific heat at constant pressure (in J/(kg K)), is the ratio of heat capacities, and the viscous stress tensor (in N/m2) and the heat flux (in W/m2) along the direction are defined, respectively, as where is the dynamic viscosity (in kg/(m s)), is the thermal conductivity (in W/(m K)), and .
Viscous dissipation in energy equation is ignored according to low values of velocity. Additionally, the ideal gas law is used to compute density variations as where J/(kg K) for air. Because of the temperature differences, the dynamic viscosity () and thermal conductivity () are computed by the Sutherland law: where and are the dynamic viscosity and the thermal conductivity at . For air, the following values are considered: kg/(m s), W/(m K), K, and the Sutherland constants are K and K. This choice is valid for a range of temperature values from K to 1900 K approximately.
2.2. Numerical Approach
The solution method is an adaptation of a fully implicit formulation based on a time-accurate algorithm proposed by Sewall and Tafti [20]. The local density variation is coupled with both temperature and pressure variations, whereas other properties are only coupled with temperature variations. No low Mach number assumption is used in this formulation. Sewall and Tafti [20] successfully applied this method to natural convection in a differentially heated square cavity flow and in a Poiseuille-Bénard flow with large temperature differences.
The numerical approach herein is based on a second-order finite difference formulation with a fully staggered arrangement. The primitive variables as pressure, density, and temperature are located at the middle of the cell and velocity components are shifted at the middle of the cell sides. The coupled set of equations is solved using an iterative predictor-corrector procedure at each time step. During the prediction step, the momentum conservation equation (2) and the energy equation (3) are advanced implicitly using the Crank-Nicolson method for the advection and diffusion terms. The divergence term in (3) is treated implicitly. The numerical solution of these equations is considered satisfactory when the residual, summed over the whole calculation domain, is smaller than .
The coupling between the mass and the momentum conservation equations is completed through a pressure correction (), determined from the following elliptic Helmholtz equation: where is the discrete time step (in s) and the superscripts and denote, respectively, time (or outer) and inner iteration. The symbol () means an intermediate or predicted value. The solution of the elliptic Helmholtz equation (8) is obtained using a V-cycle multigrid method with the strong implicit procedure ILU as smoother [21]. The numerical solution of this equation is assumed to be converged when the normalized residual, summed over the whole calculation domain, is smaller than . The different steps of the algorithm are summarized in Table 1. The numerical approach is further documented in Barrios-Piña [22] and a study in a backward-facing step configuration is shown in Barrios-Piña et al. [23]. The numerical simulations were carried out on a cluster IBM-xeon, 14xbi-quadX3550 3 GHZ, 16 and 32 GbRAM.
3. Geometry of the Flow and Initial Conditions
Figure 1 illustrates the configuration of the jet flow under study. The flow is a laminar coflowing plane jet from a nozzle of thickness . The length and the height limit the study domain. The fluid of the jet and the ambient is considered air. The jet temperature is greater than the ambient temperature (positive buoyancy). The Reynolds and Richardson numbers which characterise the flow are defined as where is the gravity, is the jet kinematic viscosity (in m2/s), , is the jet velocity (in m/s), and is the coflow velocity (in m/s). The relative jet density is defined by the relationship jet density to reference density (in kg/m3). The coflow velocity must be smaller than the jet velocity to not affect its natural expanding into the ambient. Therefore, we considered a relationship [10]. A coflow velocity is considered here to ensure numerical stability and to control jet propagation out of the domain.
For the initial thermal state, the aerostatic equilibrium of the atmosphere is considered, where pressure, temperature, and density change with the elevation. From the equation of fluid statics, and from the following expression of density, the vertical variations of pressure, density, and temperature as function of elevation are written as where the index implies properties at ground level defined at (properties of reference), the index denotes properties of the adiabatic atmosphere in steady state, and is defined as .
This state of reference, which is adiabatic and aerostatic, considers an environment where pressure, density, and temperature are not constants and vary with elevation, as happens in the real lower atmosphere.
4. Boundary Conditions
From Figure 1, three different kinds of boundary conditions are identified: the inflow boundary condition (AD), the outflow boundary condition (BC), and the lateral boundary conditions (AB and CD). The treatment of boundary conditions is not trivial in this kind of flow configuration. The control of flow rates through the open boundaries and the pressure condition at these boundaries are the cause of the difficulty. Numerous works in the literature have focused on this topic; however, there is not a general methodology to ensure good prior boundary conditions.
The characteristic boundary conditions applied by Thompson [24, 25] and later discussed and modified by Poinsot and Lele [26] are widely used for the treatment of open boundaries, especially when acoustic waves appear in the flow like in supersonic flow regimes [27]. The basic idea of characteristic boundary conditions is to split the convective terms in the boundary-normal direction into several waves with different characteristic velocities and then express unknown incoming waves as a function of known outgoing waves. Thus, to apply this kind of boundary conditions, the strength of the waves entering the computational domain must be determined, which is not trivial [28]. Moreover, the characteristic boundary conditions particularly work well in numerical algorithms based on density, which generally use explicit schemes. For algorithms using implicit formulations based on pressure, like the present numerical model, the implementation of the characteristic boundary conditions is more complex.
Another kind of boundary conditions widely used for the treatment of open boundaries is the absorption boundary condition known as PML (Perfectly Matched Layer) [29]. The PML technique considers an absorbing zone which is theoretically reflectionless for multidimensional linear waves at any angle of incidence and any frequency. Although the PML technique itself is relatively simple when it is viewed as a complex change of variables in the frequency domain, it is important to note that, to yield stable absorbing boundary conditions, the phase and group velocities of the physical waves supported by the governing equations must be consistent and be in the same direction. Thus, special treatment needs to be performed to improve absorption of waves or the different flows that impinge the base of the PML at grazing incidence and travel inside the PML over long time periods. Moreover, to implement the PML boundary condition, additional mesh points are necessary to consider the absorbing region. The fact of considering more mesh points will result in a CPU-time increase, which is not desirable, especially for fine mesh numerical simulations. Some other techniques have been developed to improve the implementation of the PML boundary conditions, for example, convolutional formulation of the Perfectly Matched Layer [30].
Given the complexities mentioned above, the subsonic regime of the present jet flow, and expected negligible acoustic wave effects, some numerical tests using typical boundary conditions like Neumann, Dirichlet, and convective were performed to define the one to be imposed on each side of the domain (see Figure 1). Thus, the boundary conditions retained here are the following.
Inflow Conditions (AD). For the vertical velocity component and for temperature, hyperbolic tangent profiles were imposed: where is the jet temperature (in K), is the reference temperature at ground level (in K), , is the momentum thickness (in m), and is imposed to be . In this way, the initial profiles adopt a flattened shape and are symmetric with respect to the jet axis. The horizontal velocity component is considered to be zero.
For pressure, the equation of fluid statics is imposed: and density is obtained from the ideal gas law.
Outflow Boundary (BC). At the upper boundary, convective conditions, also known as Sommerfeld radiation conditions, are considered. Convective conditions relate the temporal derivative and the normal derivative of the unknown in case of boundaries far away from sources and normal to the propagating waves. Its application in cases of nonperpendicular wave incidence and boundaries close to sources is doubtful, and therefore, this boundary condition has to be understood as an approximate boundary condition to avoid wave reflection in such cases. However, the numerical implementation of this boundary condition is not difficult compared with the boundaries mentioned before. The upper boundary condition is considered as follows: where is the variable computed at the boundary and is the phase velocity. To estimate the phase velocity, the Orlanski method [31] is used. First the phase velocity is computed locally with where denotes the last mesh point along the direction, normal to the boundary, and is the time iteration. Once the phase velocity is computed, the following expression is used to calculate at : where must be satisfied. At the limit , expression (17) yields and when , or In this way, no information comes from the internal solution; therefore, this condition must be considered as information imposed from outside (by imposing an appropriate value or a value from the previous time step). It is noticeable that for each variable corresponds to an expression of form (15); thus, there exist different phase velocities that must be computed separately for each variable.
To keep consistency with the inflow boundary condition for pressure, the equation of fluid statics is still imposed: and density is obtained from the ideal gas law.
Lateral Conditions (AB and CD). The lateral conditions implemented here allow the contribution of mass across boundaries and also the pressure control. These boundary conditions are known in the literature as traction-free boundary conditions. This kind of boundary conditions has been used for incompressible jet flows [32, 33]. In the present work, the traction-free boundary conditions have been reformulated to consider compressibility. The principle is where the stress tensor is and is the unit normal to the boundary. Thus, the boundary conditions for the velocity components are The boundary condition for pressure is obtained from the following mass conservation equation combined with (24) as follows: The traction-free boundary conditions are stable when the cell Reynolds number is less than 2 at the boundary [32]. This restriction can be satisfied by mesh refinement at the near boundary region. However, it may not be necessary for the present configuration because of the low velocity values obtained away from the jet influence.
The lateral boundary condition imposed for temperature is the homogeneous Neumann boundary condition of the form Finally, density is obtained from the ideal gas law.
5. Initial Conditions and Simulation Parameters
For the numerical simulations, and are and , respectively. The mesh resolution is 195 × 259 cell points and is refined in shear zones and near boundaries (see Figure 2).
The initial conditions given by (12) are imposed for initial pressure, density, and temperature, respectively. The velocity components are zero all over the computational domain: . Convergence is supposed to be reached when is less than and the time step is s, in all cases.
The state of reference at ground level defined at is Pa, K, , and m/s2.
6. Results and Discussions
The numerical simulations (cases S1 and S2) focus on the analysis of the jet evolution by taking into account two different Richardson numbers. The first one corresponds to a temperature difference where the pertinence of the Boussinesq approximation may be questionable. The second case corresponds to a temperature difference significantly higher, where the use of the Boussinesq approximation is no more applicable. The interest of this analysis is to test the behaviour of the boundary conditions and to evaluate by comparison of results the relevance of using the present numerical model against the Boussinesq approximation.
Results for the analysis of instabilities at the near-nozzle region are finally shown, where a stability diagram was obtained to specify the onset of the unsteady state. For this reason, some numerical simulations were carried out by considering the Reynolds and the temperature difference , as control parameters.
6.1. Analysis of the Jet Evolution
The physical parameters considered in the numerical simulations are shown in Table 2.
We observe that for the physical parameters considered in the numerical simulations the flow is unsteady and presents two-dimensional sinusoidal oscillations, known as sinuous mode. Thus, the following discussion is based on averaged values calculated with an integration time of 10 s.
Figure 3 shows the time-averaged vertical velocity at the jet centreline for cases S1 and S2, Figures 3(a) and 3(b), respectively. In both cases, the centreline vertical velocity remains constant in a region close to the nozzle. In this region, called potential cone, the inertial forces are dominant over the buoyancy forces. The length of the potential cone in case S2 is smaller than in case S1. To a higher height, the centreline velocity increases to a maximum magnitude, which implies that buoyancy forces become dominant over the inertial forces: this zone is known as plume zone. The maximum velocity obtained in case S2 is greater than in case S1. Beyond the plume zone, the centreline velocity decreases with height, which implies that inertial forces then become preponderant. By comparing these results against those obtained with the Boussinesq approximation, we observe that the velocity is overestimated when using this approximation. The difference between the maximum velocity magnitudes in case S1 (Figure 3(a)) is less important than in case S2 (Figure 3(b)).
(a) Case S1,
(b) Case S2,
Figure 4 illustrates the decay curve at the centreline of the time-averaged temperature. We observe in both cases a region where the temperature remains constant. This is the region of the temperature cone. According to Figures 4(a) and 4(b), the height of the temperature cone is almost the same for both cases. The Boussinesq approximation slightly overestimates this height in case S1 and it is even more significant for case S2.
(a) Case S1,
(b) Case S2,
Concerning the effects of boundary conditions, in case S1 both centreline velocity and temperature normally develop to the outflow (the upper boundary). However, in case S2 where the temperature difference is large, there is an evident problem with the boundary conditions due to the unexplained increase of velocity and temperature. This effect, inherent to such boundary conditions, could be dissipated either by mesh refinement near boundary or by implementing a buffer zone.
Figure 5 shows the time-averaged vertical velocity profiles at three different heights: at the nozzle level, , and . A nearly flat shape characterises the lower velocity profile. Away from the nozzle, the velocity profiles show a Gauss-bell form, which gradually expand laterally as height increases. For a height of , the velocity increases at the jet centre, and then it decreases, as can be observed in the velocity profile.
(a) Case S1,
(b) Case S2,
For both cases, in the near-nozzle region up to , the velocity profiles obtained with the Boussinesq approximation are similar to those computed without approximation. However, as height increases beyond , the maximum velocity magnitude at the jet centre is overestimated by the Boussinesq approximation for case S1. Moreover, for case S2 the profile obtained by the Boussinesq approximation is more flattened than that computed without approximation.
Longitudinal temperature profiles are shown in Figure 6. Temperature profiles imposed at the inflow adopt similar flat shape compared to the velocity profile. In both cases S1 and S2, the temperature profiles gradually spread as height increases; see Figures 6(a) and 6(b), respectively. The temperature at the jet centre also decreases with height. The temperature profiles obtained with the Boussinesq approximation are similar to those obtained without approximation in the near-nozzle region (); however, an overestimation of temperature at the jet centre is evident, especially when using Boussinesq approximation for case S2. Far away from the nozzle, the temperature profiles obtained with and without the Boussinesq approximation are closed for case S1. On the other hand, a significant difference between both profiles obtained for case S2 can be observed in Figure 6(b). The profile obtained with the Boussinesq approximation shows two peaks or maximum temperature values, whereas the profile without approximation has a Gaussian shape with a single maximum temperature value.
(a) Case S1,
(b) Case S2,
Finally, Figure 7 shows the instantaneous temperature fields obtained without approximation for both cases. The presence of the temperature cone which rises in the near-nozzle region is observed. In the plume zone, the instabilities develop in a sinuous mode, where the plume is more extended for case S2.
(a) Case S1,
(b) Case S2,
6.2. Analysis of the Onset of Instabilities at the Near-Nozzle Region
This section focuses on the unsteady-state onset at the near-nozzle zone. For this purpose, a number of simulations were carried out by considering two control parameters: the Reynolds number, which characterises the intensity of inertial forces, and the temperature difference , which characterises the intensity of buoyancy forces. To observe the occurrence of instabilities on the flow, a time signal analysis of the velocity components was carried out at different points distributed along the symmetry jet axis.
The Reynolds number for which the flow becomes unsteady, to a given temperature difference, is assumed when the maximum amplitudes of the velocity time series were greater than % of the jet velocity at the nozzle, . Within the range of Reynolds numbers considered, only the sinuous mode was observed on the flow.
Figure 8 shows the stability diagram obtained from the numerical simulations. The diagram graphs the Reynolds number versus the temperature differences. The critical Reynolds number is about for . Then, as the temperature difference increases, the buoyancy forces become progressively dominant against the inertial forces. This fact causes the onset of the unsteady state for Reynolds numbers considerably lower than the isothermal critical Reynolds number. For a temperature difference about 25 K, the flow is unsteady for a very low value of the Reynolds number (). This allows us to conclude that beyond K, which corresponds to a high Richardson number (), the flow is governed by buoyancy forces, whereas the inertial forces are nearly negligible. In this case, the Reynolds number does not seem to be a relevant parameter to characterise the flow. Consequently, since the driving force is mainly due to buoyancy (pure natural convection), the use of the Grashof number seems to be more appropriate to study such a flow.
The onset of instabilities shown in Figure 8 is supposed to be reached by considering some tolerance. This tolerance is illustrated by black dots and unfilled triangles. The black dots denote the Reynolds number for which the flow is steady, whereas the unfilled triangles show the Reynolds number for which the flow is unsteady. The region under the black dots is characteristic of the steady flow regime and the region above the unfilled triangles corresponds to an unsteady flow regime. It can be observed that the Reynolds number based on the jet velocity () and the nozzle thickness () loses its relevance in a system of pure natural convection.
7. Conclusions
A free vertical plane jet under large temperature gradients is presented. The first tests allowed us to define the best adapted boundary conditions for the flow under study. At the outflow boundary, the convective conditions were implemented and the Orlanski method was used to compute the phase velocity. The lateral boundaries were treated by using traction-free boundary conditions and were reformulated to consider density variations.
Two numerical simulations were conducted for studying the jet evolution and boundary conditions by using different Richardson numbers. The convective conditions imposed at the upper boundary (outflow) show an undesirable effect in a zone close to the boundary, especially when the Richardson number is important. However, this undesirable effect can be attenuated either by a mesh refinement near boundary or by implementing a buffer zone. In any case, this boundary should be pushed as far as possible from the nozzle.
A comparison with results from the Boussinesq approximation was also conducted to discuss the relevance of the present numerical model which solves the complete governing equations. It appears that the Boussinesq approximation overestimates the velocity and temperature. Disagreements between the results obtained with and without Boussinesq approximation were observed, especially when the temperature difference becomes significant (on the order of 100 K).
Finally, an analysis of the onset of unsteadiness in the near-nozzle region shows that the two control parameters need to be considered: the Reynolds number and the temperature difference . For the unsteady state, the critical Reynolds number decreases almost linearly when the temperature difference increases. Beyond K, the flow is essentially governed by the buoyancy forces and the inertial forces are nearly negligible.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.