Mathematical Problems in Engineering

Volume 2015, Article ID 891894, 11 pages

http://dx.doi.org/10.1155/2015/891894

## Analysis of Flow Evolution and Thermal Instabilities in the Near-Nozzle Region of a Free Plane Laminar Jet

^{1}Tecnológico de Monterrey, Avenida General Ramón Corona 2514, 45201 Zapopan, JAL, Mexico^{2}Instituto Mexicano del Petróleo, Eje Central Lázaro Cárdenas Norte 152, Colonia San Bartolo, 07730 Ciudad de México, DF, Mexico^{3}Centro de Desarrollo Aeroespacial del Instituto Politécnico Nacional, Belisario Domínguez 22, 06010 Ciudad de México, DF, Mexico

Received 6 November 2014; Accepted 13 December 2014

Academic Editor: Oluwole Daniel Makinde

Copyright © 2015 Hector Barrios-Piña et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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/m^{3}), 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/m^{2}) and the heat flux (in W/m^{2}) 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.