Numerical and experimental study of a jet impinging with axial symmetry with a set of heat exchanger tubes

The main purpose of this research is to predict the flow field of impinging turbulent jet with cylinder in a 3-d state. Nowadays, shear stress has multiple uses in the industry. The study has been conducted in both compressible and incompressible states with output velocities between 100 and 150 m/s, in different eccentricities of jet with respect to the first cylinder, and in various nozzle outlet orifice distances from front edge of the first cylinder. Pressure distribution and shear stress on cylinder surfaces have been determined and efficiency of jet in cleaning the heat exchangers pipes has been analyzed. Also, an experimental investigation has been conducted in order to verify the accuracy of the numerical results. Results show that if the distance between nozzle outlet orifice and the front edge of the first cylinder (L ) equals 1.52 D (D is the diameter of the cylinder), the jet has the highest cleaning effect.


Introduction
In 1992, the impinging of laminar jet with a cylinder was studied [1]. In 1996, the impinging of turbulent jet with a circular cylinder was also studied. In this study, the optimal distance from the cylinder to obtain the maximum level of cleaning, and the position of the separation points in states with eccentricity and without eccentricity, obtained theoretically and experimentally [2].
Also, in 1998, a impinging of two dimensional turbulent jet with a cylinder was studied at different eccentricities using various turbulence models [3]. In the same year, impinging of two dimensional turbulent jet with two cylinders, which was relative to the cylinder axis with a 45degree angle, was studied [4]. Another researcher in 2000 examined numerically and experimentally the collision of two-dimensional turbulent jet with two parallel cylinders [5]. In 2001, one of the researchers studied the collision of a jet with axial symmetry with a cylinder in the direction of the perpendicular to the cylinder axis and studied numerically and experimentally at the various eccentricities. He concluded that in the state of L = 2.51 D, jet has the most cleaning effect [6].
In the present work, the impact of a jet with axial symmetry with four cylinders in the direction perpendicular to the axis of the first cylinder and in the different eccentricities to the first cylinder is numerically and experimentally studied. The main advantage of this study is that the flow field is three-dimensional and the effect of the cylinders on the cleaning of the other cylinder surfaces has been studied.
In this paper, high speed jets are considered, which they can be used as separating jets. The main purpose of this work is to analyze the efficiency of jet in cleaning the surfaces, particularly the surfaces of heat exchangers pipes. Heat exchangers are devices that provide heat transfer between two fluids with different temperatures, which are separated by a solid wall. Heat exchangers have applications in building heating systems, air conditioners, plants, petroleum industries, waste heat recovery, and chemical processes.
In using heat exchangers, after a working period, heat transfer surfaces are covered by sediments caused by fluids of the system. These sediments cause an additional heat resistance, which may decrease the performance of the heat exchanger to half of its initial value. Therefore, the quality of cleaning the heat transfer surfaces is of importance for performance recovery of the exchanger. The aim of cleaning effect is the amount of shear stress on the surface of cylinders.
Usually, the fluid that has a tendency to precipitate is considered inside the exchanger pipes; because high velocity and smooth surfaces reduce the formation rate of the sediment layer, but if the conditions and factors in the design of the exchanger cause the fluid to flow on the outer side of the outer shell, a mechanical cleaning method is used to remove sediments.

Governing equations of fluid flow
The equations governing fluid flow are the three equations of conservation or transfer, which can be briefly stated as: r r Equation (1), is a continuum equation that derives from the conservation principle of mass, equation (2), are the Navier-Stokes equations derived from the conservation principle of momentum and equation (3), is an energy equation that results from the conservation principle of energy. Turbulent flows are determined by fluctuating velocity field. Continuity equation and Navier-Stokes equations with unsteady term are momentarily true for turbulent flows. In the turbulent flow, every quantity, like f, contains average and fluctuating components, where the average component is independent of time, but the fluctuating component depends on time. Therefore, it can be stated as: Here, the quantity of f can be velocity, pressure, temperature or every other desirable quantity. Such as: v Now, if the continuity equation in the incompressible fluid is considered, taking average of the total equation results in: Then, if (5) be subtracted from the main equation (momentary), it can be shown as: Therefore, the continuity equation is separately true for both the average and fluctuating components. But, if the fluid is compressible, the aforementioned content cannot be true. Because the expressions containing r 0 u 0 i , relate the two equations to each other. We mainly consider (5).
In order to obtain momentum equation in the state of turbulent flow, Navier-Stokes equations for incompressible fluid are considered; instead of u, v, w, p and f, average and fluctuating components are used; and time average is taken from the total equation. Regarding that (6) is true, equations below can be written.
The left sides of equations in (7) are like Navier-Stokes equations for laminar flow, where instead of u, v, w, the quantities u, v, w are used. The pressure and friction expressions in the right side of equations are also similarly changed. The additional expressions appeared in the equations are due to turbulent fluctuating motions. This component can be defined as a stress tensor as shown below: The above stress tensor is called Reynolds stress tensor and the equations in (7) are called Reynolds equations in turbulent flow.

Turbulence models
A turbulence model is a computational procedure for making a system of average flow equations so that about a wide range of flow problems can be solved. In the twoequation models, two-equation of transfer, the partial differential equation is solved for turbulence kinetic energy k in one of the equations, and for the rate of dissipation of turbulence kinetic energy e in the other [7]. From the four turbulence models to be mentioned in the following sections, the first three models are two-equation, and the last one is seven-equation.
In the transfer equations of k-e models, G k is the generation of turbulence kinetic energy caused by average speed gradients; G b is the generation of turbulence kinetic energy caused by buoyancy effects; and Y m shows the effect of fluctuating expansion on compressible turbulence with total rate of dissipation. The calculation of each of the above-mentioned parameters is shown as follow: where S is the absolute value of average strain rate tensor and is obtained as . S ij is obtained as below: where Pr t is turbulent Prandtl number for energy. For the standard and realized k-e models, the value of Pr t equals 0.85; and for the k-e RNG models, it is obtained as where M t is the turbulent Mach number and is defined as is the speed of sound [8].

Standard k-e model
This model is a two-equation and semi-experimental model based on equations of transfer model for the turbulence kinetic energy k, and rate of dissipation e, which has been proposed by Jones and Launder [9]. In its extraction two important conditions are considered: 1) the flow must be completely turbulent and 2) molecular viscosity effects must be negligible. This model has an acceptable accuracy for industrial flows; its pros and cons are completely obvious; and has wide application in the industry. Transfer equations for k and e are as follow: where G k is the generation of turbulence kinetic energy caused by gradients of average velocity, G b is the generation of turbulence kinetic energy caused by buoyancy effects and Y m is fluctuation expansion effect in compressible turbulence with total rate of dissipation. Calculation of the above parameters has been mentioned in the previous section. C 1e , C 2e , C 3e are constant values and s k and s e are turbulent Prandtl numbers for k and e, respectively. Turbulent or eddy viscosity, m t , is calculated combining k and e as follows: where C m is a constant value.
The constant values of this model are given in Table 1.

RNG k-e modeling
This model is obtained applying a difficult statistical method (called group theory of re-nondimensionalizing). This model is similar to standard k-e model in terms of structure, but also contains the following modifications: -RNG model has an additional expression in the equation of e, which improves the accuracy especially for the flows with rapid strain; -In this model, the effect of rotation on turbulence is considered, which enhances the accuracy for rotational flows; -RNG theory gives an analytical formula for turbulent Prandtl numbers, whereas standard k-e model uses specified constant values; -While standard k-e model is a model with high Reynolds number, RNG theory gives a differential equation for effective viscosity, which is obtained analytically and considers effects of low Reynolds numbers. Though, effective use of this property depends on an appropriate treatment of the near-wall area.
These properties have caused this model, compared to standard k-e model, to be accurate and logical for wide variants of flows. A comprehensive explanation of RNG theory and its application in turbulence have been given in [10]. The transfer equations for k and e, in this model, are given as follow: where G k , G b and Y m are as the quantities presented for standard k-e turbulence model. a k and a e are inverse effective Prandtl numbers for k and e, respectively. These quantities, using RNG theory, are obtained as follows: where a 0 = 1.0. In limit state with high Reynolds number m mod m eff << 1 , we have a k = a e ≈ 1.393. The process of scale elimination in RNG theory, leads to a differential equation for turbulent viscosity as follows: whereỹ ¼ m eff m and C y ≈ 100. In order to obtain the quality of effective turbulent transfer change with effective Reynolds number (with eddy scale), (18) is integrated, therefore, this model can be used satisfactorily for the flows with low Reynolds number and the flows in the near-wall. The limit state, i.e. high Reynolds number, leads to m t ¼ rC m k 2 e , where using RNG theory, the value C m = 0.0845 is obtained. It is observable that this value is very close to the experimental value obtained for C m in standard k-e model as 0.09.
The main difference of this model with standard k-e model is in the R expression, which is obtained as follow: where h ¼ Sk e , h 0 = 4.38 and b = 0.012. The constants of this model are given in Table 2.

Realizable k-e model
The realizable k-e model is a rather recently developed model and has two main differences with standard k-e model: the realizable k-e model has a new formulation for turbulent viscosity; in this model, from a precise and complete equation for square average transfer of eddy fluctuation, a new transfer equation for the rate of dissipation of e has been obtained.
This model was proposed by Shih et al. [11]. The transfer equations of this model for k and e are as follow: where and Y m are as the quantities presented for standard k-e turbulence model. C 1e , C 2e , C 3e are constant values and s k and s e are Prandtl numbers for k and e, respectively. This model is valid for a wide range of flows such as homogeneous shear and rotational flows, free flows like jets and mixing layers, canal flows and boundary layer, and separation flows [11,12]. In realizable k-e model, the turbulent viscosity is modeled as in the standard k-e model and RNG model. The only difference is in the factor C m , which is not constant here and is obtained as follows [13]: where V ij is the average tensor of rotation intensity in a rotational framework with angular velocity as v k . Constants of the model, A 0 and A s , are obtained as:  Table 2. The constants of RNG k-e turbulence model.
It can be observed that C m is a function of average strain rate, average rotation rate, angular velocity of system's rotation and turbulence fields (k and e). The constants of this model are given in Table 3 [11].

Reynolds stress model (RSM)
This model strictly takes into account the effects of line curvature in flow, circulation, rotation and immediate changes in the intensity of strain. Compared to oneequation and two-equation models, this model has great capability in the accurate prediction of complex flows. This model has good efficiency in tornado flows, high circulation flows in combustion chambers, rotational flows in crossings and secondary flows excited by stress in tubes.
Transfer equations of Reynolds stresses, ru 0 i u 0 j , can be stated as follow [8][9][10]: where C ij , D L ij , P ij and F ij do not require to be modeled, but D T ij , G ij , F ij and e ij must be modeled.
Using scalar turbulent diffusion, turbulent diffusion transfer, D T ij , is simplified as follow [18]: Turbulent viscosity m t is obtained from m t ¼ rC m k 2 e . Lien and Leschziner have obtained the value s k = 0.82 in [18].
Considering the proposals by Gibson and Launder [14], Fu et al. [19], and Launder [15,20] the pressure strain term f ij can be separated and modeled as: where f ij,1 is slow pressure strain term (return to homogeneous state term); f ,ij,2 is rapid pressure strain term; and f v is wall reflection term and are modeled as follow: where C 1 = 1.8.
where C 2 = 0.6. P ij , F ij , G ij and C ij are defined as shown in (23) and some quantities are defined as P ¼ 1 where C 0 1 ¼ 0:5, C 0 2 ¼ 0:3, n k is k-th component of unity normal vector, perpendicular to wall, d is vertical distance to wall and C l ¼ C 3=4 m k in which C m = 0.09 and k = 0.41. This term is valid for re-distribution of vertical stresses near the wall and tends to decrease the stress perpendicular to wall and increase the stresses parallel to wall. Generation expressions caused by buoyancy effects are modeled as follows: where Pr t is turbulent Prandtl number for energy and it equals to 0.85. The turbulence kinetic energy, considering trace of Reynolds stress tensor, is calculated as follows: In order to obtain boundary conditions for Reynolds stresses, the transfer equation for turbulence kinetic energy is extracted as below: where s k = 0.82. Although the above equation is solved generally for the whole domain of flow, the values obtained for k are used just in boundary conditions. In other conditions, k is obtained from (30). Tensor of dissipation e ij is modeled as: where Y M ¼ re2M 2 t is an additional term as expansion dissipation [21] and M t is turbulent Mach number. The scalar quantity of rate of dissipation e is obtained as follows: where s e = 1.0, C e1 = 1.44, C e2 = 1.92 and C e3 is a function of local flow direction with respect to gravitational vector.

The treatment of near-wall for turbulent flows surrounded by wall
In order to analyze near-wall treatment the non-equilibrium wall functions and logarithmic law of average velocity sensitive to pressure gradient are used, which are as follow: where y y is physical viscous sub-layer thickness and y Ã y ¼ 11:225. In order to calculate the generation and the dissipation of the turbulence kinetic energy in near-wall elements, the non-equilibrium wall function uses the concept of two-layer which is required in solving equation k in near-wall elements. It is assumed that the near-wall elements include viscous sub-layer and fully turbulent layer. For turbulence quantities, the following profile assumption is applied: where C l ¼ kC À3=4 m and y y is viscous sub-layer non-dimension thickness. Applying these assumptions, average generation (k) of element (G k ) and average dissipation rate of element (e) can be obtained from volumetric average of G k and e of near-wall elements. For quadrilateral and hexahedral elements, volumetric average is obtained as: where y n is the height of element defined as (y n = 2y p ). For the elements with other shapes (such as triangular and four-dimensional elements) appropriate volumetric averages are utilized.
In (37) and (38), generation and dissipation of turbulence kinetic energy for near-wall elements are sensitized effectively to the proportion of viscous sub-layer and fully turbulent layer which is widely changed from an element to other in the high non-equilibrium flows. It effectively undermines the assumption of local equilibrium (dissipation = generation), which is considered by standard wall functions in calculation of generation and dissipation of turbulence kinetic energy in near-wall elements. In fact, non-equilibrium wall functions, to some extent, consider neglected non-equilibrium effects of standard wall functions in calculations. Therefore, in complex flows, when the average flow and turbulence are subjected to severe pressure gradients and rapid changes and the flow contains separation, reattachment, and impingement, it is better to use non-equilibrium wall functions.

Numerical model
In this work, diameter of nozzle outlet orifice is d = 0.0125 m (0.5 inch), diameter of cylinders D = 0.0254 m (1 inch) and length of cylinders is 40 cm. Exit air velocity from nozzle in different stages is considered between 100 and 150 m/s and the flow of jet is in stable state. For the above nozzle velocities and assuming the sound velocity in air as 343.2 m/s, in ambient temperature as 20 C, the Mach number of output flow from nozzle is obtained according to Table 4, where in the analyses it is assumed that the flow is incompressible in velocity range of 100-120 m/s and it is compressible in the velocity range of 130-150 m/s.
In the study of impinging jets, the critical Reynolds number is calculated based on diameter of nozzle and its output velocity. This number for axis symmetric jets is in the range of 3000-14000, where in this study the calculated Reynolds number, based on output velocity as 100 m/s, is about 71000; therefore the jet is fully turbulent in all the states.

The arrangement of cylinders
In arranging the pipes, different factors play a role, where the heat transfer is one of the most important factors. In this work, the cylinders are arranged in way that to have the highest heat transfer. There are two arrangements for cylinders: 1) Linear 2) Staggered.
The set of pipes is specified by relative sloping step between axes of pipes. Based on the previous studies by researchers, from the viewpoint of heat transfer, the staggered arrangement is very effective. In this work, the values a and b are determined considering Reynolds number and the number of cylinders [16,17].

Numerical solution field and computational network
The numerical solution field used in this work is a 3-d field, which is shown in Figure 2. Considering Figure 2 and division values of field boundaries given below, the computational network is generated Table 5. Figure 3 shows the computational network, the 486 elements, and the nodes.
It should be mentioned that in selecting the distance of the nodes, firstly the problem is solved for a large network, and then gradually the network is made smaller until the field dimensions of 0.635 Â 0.4 Â 0.4 is determined as an appropriate numerical solution field; and this field was used in numerical solutions. Also, different networking was applied on the model, and finally the determined network was used. This network satisfies the important condition that solutions and results must be independent of the network.

Boundary conditions
Applying appropriate boundary conditions plays an important role in the process of numerical solution. Conditions in field input i.e. nozzle outlet are very important in the prediction of velocity in the center of jet and shear stresses. Velocity profile in the nozzle outlet is affected by the designing properties of nozzle such as proportion of length to diameter, shape of nozzle, possibility or impossibility of the presence of transfer region in upper-hand flow of nozzle. Here, the velocity profile in the nozzle outlet is assumed to be uniform [22][23][24] and the conditions of flow input in computation domain are stated as [25]: where R e DH is Reynolds number based on hydraulic diameter, A is cross section and P is perimeter length. Table 6 gives the input conditions for the velocities between 100 m/s and 150 m/s. Therefore, boundary condition of input velocity is applied on the G surface.
Whereas a pressure equal to atmospheric pressure is applied to the surfaces A, B, C and D, therefore, the boundary condition applied to the surfaces A, C and D is defined as input pressure boundary condition and the boundary condition applied to the surface B is defined as output pressure boundary condition. These conditions are applied considering the real conditions of the test system. Also, taking into account the no-slip condition, the velocity on rigid surfaces such as the surface of cylinders and their holder walls, i.e. E and F, are considered as zero [26]. Therefore, wall boundary condition is defined for them.

The stages of performing numerical analyses
The problem under study is analyzed in the following states (Fig. 4): -The nozzle is put in front of cylinder, and considering L = 1.52 D, the problem is analyzed in 3-d state for the values of U n equal to 100, 110, 120, 130, 140, 150 m/s.   -Like the third stage, the jet, in five steps and in each step to the amount of 0.2 of cylinder radius has been moved to the bottom, in a way that in the last step, the axis of symmetry of the jet is tangential to the lower surface of the first cylinder. In this stage L = 1.52 D and U n = 100 m/s. -Considering L = 1.52 D, and without eccentricity, besides RSM turbulence model (which was considered as turbulence model in all the above states), the problem is analyzed with the other three two-equation models as standard k-e, RNG k-e and realizable k-e. In this state, U n is considered as 100 m/s.

Experimental model
The experimental model presented in Figure 5 is established to analyze the accuracy of the results obtained from the numerical solution. Before starting the experiments, the test platform must be located on a flat place, and be completely leveled. Also, connecting hoses related to nozzle and manometers must be controlled to completely establish the connections. After ensuring that the nozzle is placed in its place, the valve of air input from the compressor is opened, and considering the mounted outlets a b c  Magnitude and quality of networking 1, One-sided 1, One-sided 0.6, Two-sided 1, One-sided 1.1, One-sided 1, One-sided 1, One-sided

The stages of performing experimental analyses
Experimental tests have been performed in three stages. In all the stages, the distance between the nozzle outlet orifice and the front edge of the first cylinder, specified by L in Figure 4, has been chosen as 1.52 D. Also, in all the stages, the velocity of the exit air from the nozzle orifice has been chosen as 100 m/s. The stages of performing the test are as follows.
-The nozzle is located exactly in front of the first cylinder.
At this stage, there is no eccentricity and the pressure distribution on the surface of the first cylinder is measured and recorded, exactly in front of the jet. In order to measure the pressure, the first cylinder, in each part to the amount of 10 degrees, is circulated around its axis of symmetry. The process is repeated for the other cylinders. Therefore, in each testing, 4 Â 36 numbers of information are required; -In the second stage, the jet (i.e. nozzle), in five steps and in each step to the amount of 0.2 of cylinder radius is moved to the top, in a way that in the last step the axis of symmetry of the jet is tangential to the upper surface of the first cylinder. In this stage, in each step, 4 Â 36 numbers of information are recorded; -In the last stage of experimental tests, the nozzle, in five steps and in each step to the amount of 0.2 of cylinder radius is moved to the bottom, in a way that in the last step, the axis of symmetry of the jet is tangential to the lower surface of the first cylinder. In this stage, too, in each step, 4 Â 36 numbers of information are recorded.

Presentation of results
Contours of velocity are shown in Figure 7 for the two states a) without eccentricity and L = 1.52 D, U n = 100 m/s and b) without eccentricity and L = 1.52 D, U n = 150 m/s, in yz plane. As it can be observed, the velocity of the flow from the output of nozzle to impinging point is decreased gradually in the central line of jet, which is close to zero in the resting point. At this point, the pressure reaches to its maximum value and also in the state without eccentricity the velocity profile has rather uniform distribution, which is begun from the value close to zero and is reached to a maximum value, and then is approached to the zero value again. Also, in the eccentricity states, due to the Coanda effect, in spite of departing, the jet is pulled towards the cylinders, which is shown by Figures 8a and 8b.
In Figures 9 and 10, the procedures of shear stress distribution and pressure distribution, in 3-d, on the surface of first cylinder are shown, respectively. Also, in Figures 11 and 12, shear stress and pressure distributions, on the surface of first cylinder, along the axis of symmetry of jet, are presented, respectively. It should be mentioned that in all the figures, the shear stress and pressure have been non-dimensioned by velocity square of jet in nozzle outlet orifice and density. Pa is the pressure of atmosphere. It can be deduced from the figures that when the shear stress has the lowest value, the pressure has the highest value, because at these points velocity gradients have their     Figures 14a and 14b show the comparison between shear stress on the first cylinder and exactly in front of the jet for the states with different eccentricities between zero and 1R, where L = 1.52 D and Un = 100 m/s, and for the movement of the nozzle to the bottom or top, respectively. It can be observed that the shear stress, compared to other states, is the highest in the without eccentricity state. Also, considering the curves of the velocity level, it can be seen that applying eccentricity, the flow passing region moves to the back of the cylinder.  models give predictions close to each other, and the difference appears after the separation point and in the flow passing region, where the comparison between numerical and experimental results shows that the RSM model has higher accuracy among all the models, however, its convergence speed is lower than other models, because in RSM model, more equations are solved. Also, results show that the standard k-e model is weak in predicting the reverse and circulation flows, and the flow passing region.
In performing the tests and measuring the pressure some errors, such as error in adjusting the exit air velocity in the nozzle orifice, error in accurate adjusting of eccentricity, error in reading values and etc. occur that must not be neglected. Though, it is tried to lower the above errors as much as possible.

Conclusion
Considering the results obtained from the state of without eccentricity, it can be observed that for the optimal distance L = 1.52 D, in 3-d, the shear stress has its maximum value. In other words, in the stated distance, the jet has the highest cleaning effect, and by making nozzle distant from or close to the first cylinder, the shear stress becomes lower than the optimalstate. This is true for both the states of incompressible fluid (U n =100 m/s) and compressible fluid (U n = 150 m/s).
The results show that when applying eccentricity to the first cylinder, in the state of without eccentricity the shear stress is lower than other states. Also, applying different eccentricities to the first cylinder, the separation point moves to the back of the cylinder and so the jet affects a large part of the first cylinder. In other words, the areas that were in the flow passing region and had not been cleaned, applying the eccentricity, went out of the flow passing region and were cleaned. The presence of other cylinders shrinks the domain of the flow passing region and causes a large part of the back of the cylinder to be cleaned.
As it is obvious from the results, for the first cylinder in a specified distance of nozzle to the first cylinder, increasing the output velocity of the nozzle causes the cleaning effect of jet on the surface of first cylinder to be increased. Because, increasing the nozzle output velocity, the shear stress and consequently the cleaning effect on the surface of the first cylinder increased. Therefore, it can be stated that in order for cleaning the heat exchangers pipes, a fluid of jet (such as air jet) with high velocity can be used, so that it moves parallel with the diameter of the cylinder and sweeps the surface of the cylinder.