Different numerical methods in the study of passive scalar transport in a pipeline x-junction

Different numerical methods in the study of passive scalar transport in a pipeline x-junction W. Vicente * , M. Salinas-Vazquez, C. Chavez, E. Carrizosa Instituto de Ingenierı ´a, Universidad Nacional Auto ´noma de Me ´xico, Ciudad Universitaria, 04510 Me ´xico DF, Mexico Received 1 April 2006; received in revised form 1 January 2008; accepted 23 January 2008 Available online 8 February 2008 Abstract A computational fluid dynamic model is used to analyze the transport processes of a passive scalar generated in the mixing of two fluids or flows in a pipeline x-junction. Turbulent flow field is computed for the merging of streams using two and three-dimensional simulations, which are achieved employing Cartesian coordinates, BFC and a cut cell method. These different numerical solving methods are compared. The numerical model is validated through an experimental setup. Different parameters are measured for various operating conditions. The influence of the angle between pipe inlets is studied to establish the optimal condition in which the passive scalar concentration in both outlets is similar. Ó2008 Elsevier Inc. All rights reserved. Keywords: Concentration; Pipes; Potable water; x-Junction; CFD 1. Introduction Turbulent mixing of many fluid streams in pipelines is of much interest due to its wide application in industrial and power generation processes. Of particular interest for this work is the study of potable water distribution systems, where water from different sources is mixed. In addition, the dispersion of some disinfectant or polluting substances in the water is also analyzed. Different T-junction type pipelines have been studied numerically and experimentally [1–4]. However, a reduced number of studies are available for a confined cross-flow, which requires a more complex geometry and turbulent mass transfer simulation. The geometry studied here is a two-pipe configuration with an x-like crossing point. Both pipes transport clean water, but a specified concentration of a passive scalar is added in one of the inlets. An experimental set-up was designed for this work, consisting of an x-junction with an angle between the two pipe inlets a=90°. The volumetric flow rate in the clean water inlet was varied to study its effect on concentrations in both pipe outlets. The passive scalar used in the experimental set-up represents the hipochlorine 0307-904X/$ - see front matter Ó2008 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2008.01.022 * Corresponding author. Tel.: +52 55 56 23 35x1113; fax: +52 55 56 16 21 64. E-mail address: [email protected] (W. Vicente). Available online at www.sciencedirect.com Applied Mathematical Modelling 33 (2009) 1248–1258 www.elsevier.com/locate/apm in the water. Different parameters were obtained experimentally and results were used for validating the numerical model. In the numerical model, the hydrodynamics was represented with the Reynolds-averaged continuity, momentum and conservative species equations. The standard kemodel [5,6] was used to consider turbulence effects within the flow. The hydrodynamic equations were solved using a finite-volume method. When flow simulations involve curved geometries or boundaries that are not aligned with the grid orientation, the use of Body-Fitted Coordinates (BFC) or staircase steps approximation in Cartesian coordinates is common. However, many problems are observed, such as: inaccurate results, high computational costs, low accuracy or iteration divergence. As an alternative approach, cut cell methods have become popular, because they have the simplicity and accuracy of the Cartesian coordinates for internal cells (cells with fluid flow), and a more elaborate method only for the boundary surfaces. Many different cut cell methods have been developed, such as: immersed boundary or surface method, embedded boundary method, or the ghost fluid method. In this work, the ASAP method (Arbitrary Source Allocation Procedure) [7,8], which determines the intersections of the complex geometry with the Cartesian grid lines and calculates the free areas and volumes of the partially blocked cells for modifying the transport by source terms, is used. Considering that numerical simulations give rise to the study of more parameters than those studied experimentally, the present work is only focused on analyzing the influence of the aangle between the pipe inlets in the outlet passive scalar concentrations. The paper is organized as follows: in the next two sections, the characteristics of the experimental set-up and the numerical model are discussed, and results and performances of different solving methods are compared. Section 4compares the numerical and experimental results for validation of the mathematical model. In Section 5, the effect of the angle between both pipe inlets over the passive scalar concentration in the outlets is analyzed. Finally, in Section 7, an outline of the most relevant results is presented as conclusion. 2. Experimental set-up The experimental set-up of the x-junction was designed and constructed in the Hydraulic Laboratory of the Engineering Institute of the National University of Mexico (UNAM). A schematic view of the experimental set-up is shown in Fig. 1. The junction was built in acrylic tube with 90°cross-shaped pipe. The pipe length was 1.23 m with an inner diameter of 0.01905 m (3/400). The junction had two inlets and two outlets, labeled inlet1,inlet2 and outlet1, outlet2 respectively. Inlet1 was supplied with a mixture of clean water and a dye substance, while inlet2 was delivered with only clean water. Fluid for both inlets was supplied by two storage tanks, provided with one 372.85 W (1/2 HP) pump each. The dye substance was blue ink. Measurements of volumetric flow rates at the inlets were taken with flow-meters, while at the outlets, volumetric tanks were employed. The concentrations were measured through a Hach spectrophotometer, which employs a color technique in terms of platinum–cobalt units. At the branches, pressure was measured with a piezometer. Tests were performed for different volume flow rates at inlet2. The volume flow rate at inlet1 was approximately equal to 0.4 dm 3 /s for all cases. The Reynolds number (based on mean inlet velocity and pipe diameter) at inlet1 was Re = 28,000. The volumetric flow and Reynolds number for inlet2 varies from 0.30 to 0.80 dm 3 /s, and from 23,000 to 50,000, respectively. The experimental uncertainties associated with the measurement of volume flow rates and the concentrations of tracer species were 3% and 4%, respectively. 3. Numerical model 3.1. Numerical procedures Different procedures were performed to represent the cross-flow. The simplest form is based on a Cartesian grid through the application of porosity, Fig. 2a, where values go from 0, blocked cell, in the solid zone without fluid flow, and 1 in the fluid zone where the entire cell is open to the flow. The second form consists of a BFC grid, where only the flow zone is considered. For this procedure, three different zones were represented and interconnected, Fig. 2b. In Table 1, we present the CPU time and memory used for an example of cross- W. Vicente et al. / Applied Mathematical Modelling 33 (2009) 1248–1258 1249 Fig. 1. Experimental set-up. Fig. 2. Different numerical procedures. 1250 W. Vicente et al. / Applied Mathematical Modelling 33 (2009) 1248–1258 flow with an angle a=90°. The resolution inside the pipes is the same for the three cases. The third case will be explained later. We observed similar results in all cases in velocity and concentration fields and we could not conclude on the best method. However, for the ASAP method, CPU time is more than 2.5 and 1.5 times the CPU time for BFC and Cartesian methods, respectively. These calculations were performed in a PC-Pentium 4 with Windows-XP operative system. The problems appear when the angle a–90°, the connection between pipes is non-orthogonal and/or the walls are not parallel to the Cartesian axis. For the blocking Cartesian grid, the pipe wall is represented in a staircase fashion, Fig. 2c. This wall representation is inadequate and results in an over predicted pressure loss and spurious predictions close the wall. In the case of the BFC grid, the problem is located in the wall connection, which increases convergence time, or simply makes it impossible to establish converged results (mainly in 3D cases). Due to these problems, and indeed a higher CPU time and memory (not restrictive for a PC-Pentium 4 computer), an alternative procedure was considered, which consists of the use of the ASAP method [7,8]. This is a cut cell method that uses a Cartesian grid for all cells except for those that are intersected by a complex–boundary pipe surface. All the advantages of using the conventional Cartesian grid method are restricted to the internal cells and a more elaborate treatment is only required for the boundary cells. In the cases where a finite volume method was used in this work, the intersections of the complex geometry with the grid lines are determined and the free areas and volumes of the partially blocked cells are calculated, then, the transport equations are modified (for these partially blocked cells) by source terms that account for the local non-orthogonal intersections between the grid and the pipe geometry. This ensures a higher degree of precision and accuracy in representing the walls’ true geometry. This method is used in this work for 2D and 3D cases. 3.2. Numerical solution The cross-flow hydrodynamics was modeled via Reynolds-averaged continuity, momentum and conservative species equations (with the over-bars which usually denote the average dropped for simplicity). The turbulent–convection term, which results from the averaging of the convection term in the conservation equations, was modeled with the standard kemodel [5,6]. For the turbulent flow of homogeneous, incompressible, viscous fluid with constant properties, mass and momentum conservation equations in a Cartesian coordinates system may be written as follows [9]: oui oxi ¼0;ð1Þ oui otþuj oui oxj ¼1 q op xi þo oxj ðmþmtÞoui oxj  ;ð2Þ where qis the fluid density, u i is the time-average velocity vector, Pis the time-average pressure, mis the kinematic viscosity, and m t is the eddy viscosity. The concentration was modeled through the conservative specie equation: oc otþuj oc oxj o oxj m Sc þmt rc  oc oxj  ¼0;ð3Þ where cis the conservative specie, Sc and Sc t are the laminar and turbulent Schmidt numbers, respectively. For the present work, the following values were assigned: Sc = 1.0 and Sc t = 0.6 [10]. Table 1 Numerical results comparison Numerical method Cells number CPU memory (kB) CPU time (s) Cartesian 12,100 9.7 153 BFC 2100 12.3 229 ASAP (Cartesian) 12,100 10.3 389 W. Vicente et al. / Applied Mathematical Modelling 33 (2009) 1248–1258 1251 3.3. Numerical details The x-junction was simulated in two and three dimensions. The hydrodynamic equations were solved using a finite-volume method. The mesh in BFC coordinates used for the 2D calculations is shown in Fig. 3. In this geometry, the General Collocated Velocity (GCV) method was used for solving the Navier–Stokes equations [11]. This method uses a segregated pressure-based solver strategy with an additional correction of cell-centre momentum velocity components. The pressure–velocity coupling was based on a linearization which is similar to the well known SIMPLE procedure. The GCV allows the creation of an unstructured multiblock grid from relatively simple structured blocks, which can be generated individually. In the mesh of Fig. 3, three blocks were made and these were linked in the center. When the angle between inlets differed from 90°, two blocks were added to form the joint, avoiding highly non-orthogonal cells. 3D simulations using BFC grids were full of complexities and the grid deformation for a–90°cases did not allow easy convergence of simulations. Cartesian coordinates provide an alternative for the 3D simulations. The Cartesian approach consists of using the ASAP method [7,8]. The ASAP allows to work with complex geometries without creating meshes that are difficult to attain convergence, by the attachment of source terms to the appropriate cells. This simplification was employed for the tri-dimensional simulation of the pipeline x-junction. Fig. 3. Numerical grid. Table 2 Flow characteristics 2D 3D Grid points 4200 360,000 Pipes length 1.23 (m) 1.23 (m) Vol. flow (inlet1) 0.43 (dm 3 /s) 0.43 (dm 3 /s) Vol. flow (inlet2) 0.03–0.82 (dm 3 /s) 0.03–0.82 (dm 3 /s) Concentration (inlet1) 1.0 1.0 Concentration (inlet2) 0.0 0.0 1252 W. Vicente et al. / Applied Mathematical Modelling 33 (2009) 1248–1258 In both cases, BFC and Cartesian grid refinement tests were carried out; the results proved to be grid independent for the above mentioned meshes. The characteristics of simulations are summarized in Table 2. Flow enters the pipe with a uniform longitudinal velocity. To study the effect of high and low velocities in the passive scalar transport, simulations with inlet volumetric flow values lower than 0.1 and higher than 10.0 dm 3 /s are performed. The incoming values of kand eare calculated from the empirical relations [12]: kinlet ¼ckð0:3uinletÞ2;ð4Þ einlet ¼ceð0:090:75k1:5 inletÞ=ð0:1DÞ;ð5Þ where Dis the pipe diameter, u inlet is the inlet velocity, c k and c e are constants. The value for both constants in empirical relationships is 1.0. To test whether different inlet kand evalues influence the concentration transport, the constants c k and c e are changed from 0.1 to 10.0. For simplicity in the presentation of results, several parameters were defined. The predicted mean volumetric flow rate Qand the mean bulk concentration Cat the outlets are computed as: C¼RucdA RAudA;ð6Þ Q¼ZA udA;ð7Þ where Ais the pipe transversal area, uand cthe mean velocity and mean concentration respectively. The inlet and outlet volumetric flows for pipe 1 and 2 are, respectively: Q i1 ,Q i2 ,Q o1 and Q o2 . They are normalized by the total inlet mean volumetric flow (Q i1 +Q i2 ). The parameter r Q is defined as the inlet volumetric flow ratio r Q =Q i2 /Q i1 . Finally, the mean bulk concentrations are normalized by the total inlet mean bulk concentration (C i1 ). 4. Validation case Several experimental measurements have been performed for a pipeline x-junction with an angle a=90°so as to validate the numerical model. Experimental results consisted of mean pressure gradients, concentrations, and volumetric flows. To test the validity of a 2D approximation, 2D and 3D simulations were performed with the characteristics shown in Table 2. A large number of simulations for much higher and lower volumetric flows than those shown in Table 2 were performed, maintaining similar values of r Q . The results differ less than 3% from the ones obtained using the original values specified in Table 2. This shows independence of the scalar concentration with the range of the here studied volumetric flows. The velocity field for the same inlet volumetric flows at both pipes, r Q = 1.0, is illustrated in Fig. 4. The momentum exchange at the junction deforms the velocity profiles and forces the passive scalar toward outlet2. In spite of this strong momentum exchange, the mean bulk velocity remains the same before and after the junction in every pipe. As expected, turbulence increases at the junction and reaches a maximum at the recirculation zone, with the consequent passive scalar mixing. Fig. 5 shows the longitudinal evolution of the mean turbulent kinetic energy. Regardless of the high level of turbulence at the junction zone, here the scalar transport from one pipe to another seems to be independent of it. This turbulence independence has been tested in simulations with values of k and eten times higher and lower than those obtained with the empirical relations, Eqs. (4) and (5). As a result of high turbulence, the concentration is more homogeneous in the fluid just after the junction, but its value remains basically the same at the outlets. For this validation case, more than 70% of the concentration incoming from inlet1 is evacuated though outlet2. After the junction, the scalar in pipe 1 travels mainly along the high velocity zone. This is illustrated in Fig. 6, where numerical and experimental results show very good agreement. When the value of r Q is increased, the flow at outlet2 is also enhanced, with the corresponding higher concentration values, C o2 .Qand Cprofiles at both outlets as a function of r Q are shown in Fig. 7. A critical point W. Vicente et al. / Applied Mathematical Modelling 33 (2009) 1248–1258 1253 is observed in Fig. 7-top, where the concentration of both outlets is the same. This point is found for a value r Q 0.23, for an angle a=90°. Working around this point, it is expected to obtain the concentration homogeneity between two flow branches. For many applications, this critical point is of relevant utility. Fig. 4. Velocity vector field for the 3D simulation. a=90°and r Q = 1.0. Fig. 5. Validation case, a=90°. Turbulent kinetic energy (k) along pipe2 longitudinal direction, r Q = 1.0. Fig. 6. Validation case, a=90°. Numerical (left) and experimental (right) concentration field. r Q = 1.0. 1254 W. Vicente et al. / Applied Mathematical Modelling 33 (2009) 1248–1258 The Qprofiles show a similar behavior to the Cprofiles. However, the critical point found here, where volumetric flows at both outlets are the same, is displaced to r Q = 1.0, since both pipes transport the same fluid with equal volumetric flow. Analyzing the passive scalar profiles, it may be concluded that most of the fluid evacuated by outlet2 came through inlet1 (more than 70%) and vice versa. The pressure losses, defined as DP=(P inlet P outlet )/P inlet , as a function of r Q , are illustrated in Fig. 8. The difference between both numerical configurations is comparable with measurement errors in experimental data, lower than 4%. Higher differences are mainly found for the highest r Q values studied. The two dimension approximation proved to be a good way to save computational resources for this kind of study, which requires a high number of numerical data. Just for the validation case, more than 50 simulations were performed. The major discrepancies between numerical and experimental data were found to be lower than 4% in terms of scalar concentration, which establishes the numerical model accuracy. 5. Influence of the junction angle a As a consequence of the good performance of the numerical model shown in the previous section, the study of configurations with different aangles was achieved considering four other angles: a=45°, 67.5°, 112.5°and 135°. Flow conditions and number of grid points were similar for every case. Only an increase of 5% in the Fig. 7. Validation case, a=90°. Normalized mean bulk concentration (top) and normalized mean volumetric flow (bottom) as a function of r Q at both outlets. outlet2 and outlet1 experimental data. —– 2D numerical simulation,  3D numerical simulation. W. Vicente et al. / Applied Mathematical Modelling 33 (2009) 1248–1258 1255 grid resolution at the junction was necessary for cases a=45°and 135°. The whole range of r Q values studied for a=90°was not possible for every angle. For a= 135°, convergence was not attained for low and high r Q values. Many other authors have also reported convergence problems in flows with more than one exit [3], due to mass conservation convergence. The normalized bulk volumetric flow and passive scalar concentration profiles are shown in Figs. 9 and 10 (a=90°included) as a function of r Q . The critical point in the bulk volumetric flow profiles takes place at r Q = 1.0 for every case, although they present profiles with different slopes. The curves with maximum and minimum slope are cases a=45°and 135°, respectively. This strong flow interaction is also observed in the passive scalar concentration profiles. For low avalues, the concentration at outlet2 increases substantially. Hence, the concentration at outlet2 tends to a value close to 1.0 faster. For high avalues, this process of passive scalar evacuation by outlet2 is slower. For the passive scalar, the value of the critical point, r Qcrit , is not the same for every junction angle. Before this point, outlet1 transports a higher concentration than that of outlet2. After this point, the behavior of the scalar transport is inversed. In Fig. 10, it is possible to observe that as aincreases, so does r Qcrit value for crit- Fig. 8. Validation case, a=90°. Mean pressure loss as a function of r Q at both pipes. outlet2 and outlet1 experimental data. —– 2D numerical simulation,  3D numerical simulation. Fig. 9. Different junction angle. Normalized outlet mean volumetric flow as a function of r Q .—–a=45°; a= 67.5°;---a=90°;– a= 112.5°and – – a= 135°. 1256 W. Vicente et al. / Applied Mathematical Modelling 33 (2009) 1248–1258