CHEMICAL ENGINEERING TRANSACTIONS VOL. 52, 2016 A publication of The Italian Association of Chemical Engineering Online at www.aidic.it/cet Guest Editors: Petar Sabev Varbanov, Peng-Yen Liew, Jun-Yow Yong, Jiří Jaromír Klemeš, Hon Loong Lam Copyright © 2016, AIDIC Servizi S.r.l., ISBN 978-88-95608-42-6; ISSN 2283-9216 Computational Analysis of Swirling Pipe Flow Miloslav Dohnal*, Jiří Hájek Brno University of Technology, Faculty of Mechanical Engineering, Institute of Process Engineering, Technická 2, 616 69 Brno, Czech Republic dohnal@upei.fme.vutbr.cz Swirling pipe flow was investigated using advanced scale-adaptive turbulence model (𝑘 − 𝜔 SST SAS). The analysis of turbulent swirling fluid motion in pipe focuses mainly on the model prediction capabilities. For validation purposes were used experimental data from (Pashtrapanska et al., 2006). Global instability was produces by a contraction ring placed downstream of swirl generator. Predicted data were sampled at several distances from the swirl generator to compare them with experiment. Repeated simulations were performed on several computational grids involving unstructured and structured meshes. Grid tests revealed strong influence of computational grid on fluid flow behaviour predicted by the model. Most distinguishing feature of the predictions is asymmetry of mean velocity profiles, which is contrary to results obtained from experiment. 1. Introduction This paper is devoted to a numerical study of swirling pipe flows. Swirling pipe flows find application in many practical processes as reviewed in (Mitrofanova, 2003). Example of swirling pipe flow that does not employ twisted tape or any other tube insert is a novel process for calcination of diatomite (Hajek et al., 2015) or (Papoulias and Lob, 2015). The focus in this work is specifically on configurations with swirling flow prepared in controlled laboratory environment, where the inlet conditions closely resemble solid body rotation. Such flow has been studied experimentally by laser Doppler velocimetry e.g. in (Rocklage-Marliani et al., 2003), where the flow has been confined in a straight pipe without any obstructions. Very similar configuration was investigated by the same technique in (Pashtrapanska et al., 2006), where the authors however added a so-called tripping ring close after the swirl generator. The two swirling flows possess many similarities, but there are also notable differences, mainly the presence of large-scale instability generated by the contraction in the second configuration. Above mentioned experimental results provide a solid database for the validation of CFD simulations of swirling pipe flows. Experiments uncovered significant turbulence anisotropy in swirling pipe flows, and the character of these flows has been analysed in detail. Several computational studies have attempted to predict the flows in the configuration of (Rocklage-Marliani et al., 2003), while no computational studies have been so far performed to predict the configuration of (Pashtrapanska et al., 2006). This work provides a validation study for several often used as well as emerging turbulence models based on the latter configuration. One of the focal points of this analysis is the experimentally observed asymmetry and unsteady behaviour of the flow. To this end, advanced scale-resolving model with moderate CPU requirements has been used to ascertain its capability to predict the complex features of swirling pipe flows. The validation of such models that promise good accuracy at acceptable computational costs is of primary importance for engineering CFD practitioners. 2. The experimental platform Due to strict space limitation, only short review of studied geometry is presented here. For additional information about experiment, readers should consult (Pashtrapanska et al., 2006). DOI: 10.3303/CET1652127 Please cite this article as: Dohnal M., Hájek J., 2016, Computational analysis of swirling pipe flow, Chemical Engineering Transactions, 52, 757-762 DOI:10.3303/CET1652127 757 Experimental device consists of a straight pipe with 50 mm internal diameter and length of 5,000 mm (𝑙 𝑑⁄ = 100). Downstream the swirl generator, there is a contraction ring which generates global flow instability and reduces development length of the fluid flow. Swirling fluid flow is generated by rotating honeycomb with the same diameter as the pipe has. The honeycomb is made of numerous small tubes 4 mm in diameter. The strength of the swirl is denoted as rotation rate, which is defined as a ratio of circumferential velocity of the pipe wall to the bulk flow velocity (Pashtrapanska et al., 2006): 𝑁 = (𝑈𝜃 )𝑤𝑎𝑙𝑙 𝑈𝑚 (1) where 𝑁 [-] is the rotation rate, (𝑈𝜃 )𝑤𝑎𝑙𝑙 [m/s] is wall tangential velocity and 𝑈𝑚 [m/s] stands for bulk velocity. Flow rate and fluid properties are summarized in Table 1 below. Table 1: Experimental setup (Pashtrapanska et al., 2006) Experiment Working fluid – Diesel oil at 20 °C Re [-] N [-] 𝜌 [kg/m3] 𝜈 [m2∙s-1] 30,000 1 830 3.9∙10-6 Note that the data published in (Pashtrapanska et al., 2006) were before use normalized by volume flow rate. Similarly all predicted cross-sectional line profiles presented in Figure 4 were normalised by the same procedure. 3. Modelling methodology As was mentioned in previous paragraph, the experimental device consist of two main parts – rotating swirl generator and steady straight pipe. Simulation of the whole experimental device as in (Vaidya et al., 2011) would enormously increase computational time. To reduce the computational cost, flow domain was split into two separate sections, swirl generator and the pipe. 3.1 Boundary conditions The swirl generator was simulated first to obtain average turbulence intensity for the pipe inlet boundary condition. Average value of turbulence intensity was taken at outlet from the swirl generator. For preliminary simulation purposes, steady 𝑘 − 𝜔 model was used. Considering there are several interpretations of turbulence intensity, the formula for its determination is shown in Eq(2) (Wilcox, 2006). 𝐼 = 𝑢′ �̅� (2) where 𝐼 [-] is turbulence intensity, 𝑢′ [m/s] stands for root-mean-square velocity fluctuations and �̅� [m/s] denotes average velocity value. During simulation of the main pipe flow, the swirl generator was replaced by solid body rotation in the inlet boundary condition. Initial turbulence intensity value was obtained from preliminary simulation of the swirl generator and its value was set to 10.48 %. The outlet conditions should involve zero pressure gradient due to unequal pressure distribution throughout the generated vortices while SAS model is used for simulation (Menter, 2012). In this case, the outlet boundary condition is far away enough from the orifice (flow is expected to be fully developed and no significant fluid movements occur), pressure outlet with constant pressure distribution was chosen. This choice provided improved convergence of the computation at the cost of “artificial” stability in the near-outlet part of flow. Near wall treatment was realized using wall functions and non-slip boundary condition was used at all the walls. 3.2 Computational grid and discretization To analyse the impact of grid topology, which has been previously reported to have a noticeable impact on SAS model predictions in (Shim et al., 2009), several computational grids (structured and unstructured) were generated. According to (Menter, 2012), the 𝑘 − 𝜔 SST based SAS model is written in 𝑦+ - insensitive way; and sufficient grid resolution should be characterized by: Δ𝑚𝑎𝑥 ≤ 0.05𝐷 (3) where Δ𝑚𝑎𝑥 [mm] is maximal cell edge length and [mm] is representative length (here the tube diameter). The computational grids were generated in agreement with Eq(3). There was no special treatment near wall boundary since the SAS model is 𝑦+ - insensitive. 758 Two unstructured grids (medium density and fine grid) with different cell density were generated using polyhedral cells. The grids consist of 1 and 2.2 million cells, respectively. Additionally, two hexahedral computational grids with different topologies were employed; body-fitted structured hexahedral grid (1 million cells) and unstructured “hexa-core” grid (950,000 cells). Some of the computational grids are illustrated in Figure 1 below. (a) Fine polyhedral (b) Structured hexahedral (c) Hexacore (pseudo-cubic) Figure 1: Examples of used computational grids Numerical schemes used for spatial discretization of governing equations were set in line with the recommendations in (Menter, 2012). Thus turbulent quantities, pressure, and time derivative were discretized using second order differencing, while momentum was discretized using bounded central differencing. 3.3 Unsteady simulation and data sampling A preliminary steady simulation of swirling pipe flow was performed using 𝑘 − 𝜔 SST model to calculate initial flow field for unsteady simulation, performed using SAS model. However, steady simulation was able to converge only with first order discretization. From that point, the simulation were ran as unsteady without any preliminary calculated flow field. To reduce error and influence of initial conditions, the start-up period of the simulations was three residence times (6.4 seconds). After that time, the simulation was continued for next seven residence times (14.95 s) to obtain convergence of time-averaged flow field variables. Data were sampled at 8 distances from the swirl generator (l/d = 3, 10, 17.3, 37.3, 44.8, 52.3, 81.7 and 98.4) consistently with experimental measurements of (Pashtrapanska et al., 2006). However, only first three positions are presented in this work due to strict space limitation. 4. Results In this section, simulation results are presented and divided into several sections to cover all aspects of simulation. 4.1 The energy spectrum One of the basic aspects that every developed turbulent flow should exhibit is the Kolmogorov’s -5/3 law. The energy spectrum was calculated by sampling velocity magnitude from point positioned on the centerline; mean velocity value was subtracted from the signal to obtain fluctuating velocity components. (a) Medium polyhedral mesh (b) Hexahedral mesh (c) Hexacore mesh Figure 2: Resolved energy spectrum at l/d = 10 It should be noted that better understanding of wave’s dissipation could probably be provided by two-point correlation to calculate dissipation in space (not in one point) (Davidson, 2016). Comparison of resolved energy spectrum on the three grids is shown in Figure 2. The energy spectrum was normalized by Kolmogorov length scale 𝜂 and velocity scale 𝑢𝜂. 759 4.2 Unsteady flow characteristics One of the main advantages of SAS model is to resolve turbulent eddies in LES-like manner, when global instabilities are strong enough. On the other hand, when the instability is not strong enough, the model falls back to (U)RANS-like solution. The source of global instability is here represented by contraction ring. To recognize where SAS model resolving capability is turned on, contours of blending function are shown in Figure 3. The blending function can reach values in interval 〈0,1〉. When value of 1 is reached, the artificial turbulence term in governing equation disappears; and SAS model acts like 𝑘 − 𝜔 SST model. (a) Inlet and orifice area (b) Outlet part Figure 3: Instantaneous values of blending function F1 on hexacore mesh Contours of the blending function shows that contraction ring generates flow instability strong enough to switch on the SAS model (value of blending function is near zero). Red color in front of the contraction ring (orifice) is caused by homogeneous turbulence of the fluid prescribed on inlet boundary. On the other hand, situation near outlet from the domain is completely different. Fluid motion in this area is rather stabilized and SAS model partially falls back to (U)RANS-like solution. This situation starts to occur between l/d = 52.3 and 81.7. 4.3 Flow Asymmetry In Figure 4, normalized average values of axial and tangential velocity components are shown and compared with experimental data. From experimental data of axial velocity (black curve), it can be easily seen that the swirling pipe flow is perfectly symmetrical behavior near the orifice. This is quite expected behavior, as the geometry is axially symmetrical and this naturally includes also the wake behind the contraction ring. In the case of this flow, symmetry of mean values therefore is a primary validation criterion for CFD predictions. In contrast to that, simulation results predict asymmetrical flow behavior. Even more striking is the fact that the time-averaged solutions on grids of the same cell size but differing in grid topology display marked qualitative features that may be directly linked to the grid topology, which is clearly seen in Figure 5. The impact on line samples is that they display unexpected arbitrary-looking behavior. This asymmetry is noticeable near the orifice where flow instability is very strong; and boundary layer wake plays significant role in developing flow field. The best agreement with experimental data was achieved on fine unstructured polyhedral mesh. Results obtained from structured grid are rotationally periodic rather than axially symmetrical, which is however apparent only from the two-dimensional cross-sectional profiles in Figure 5 and not from the graphs in Figure 4. These results are very similar to results presented in (Shim et al., 2009) where SAS model prediction capabilities are better on unstructured grid. Only in that publication, authors did not consider this behavior erroneous. Tangential velocity component exhibits similar qualitatively unreasonable behavior. It seems to be less affected by computational grid; however more detailed investigation is necessary. It is important to note that RANS simulations (not documented in the graphs in this manuscript) performed for this flow display very little impact of particular grid topology. It thus seems that the SAS model may suffer from a previously undocumented deficiency, which causes the results to be influenced by grid topology to an unprecedented degree. 760 (a) Mean axial velocity at l/d = 3 (c) Mean tangential velocity at l/d = 3 (b) Mean axial velocity at l/d = 10 (e) Mean tangential velocity at l/d = 10 (c) Mean axial velocity at l/d = 17.3 (f) Mean tangential velocity at l/d = 17.3 Figure 4: Comparison of normalized average axial and tangential velocity components 4.4 Contours of average velocity magnitude Discrepancy of fluid behavior (related to computational mesh topology) with 𝑘 − 𝜔 SST SAS model is documented in Figure 5. Structured hexahedral mesh yields square-shape of the flow pattern in the central part of the pipe. Hexa-core mesh provides slightly better results, however the fluid is still moving with unrealistic features. Completely different (also unrealistic) fluid motion exhibit the results on medium-density polyhedral mesh. It should be noted that velocity magnitude is superposition of axial, radial and tangential velocity components; velocity profiles shown in Figure 4 were sampled on x-axis. Comparing Figure 5 and Figure 4 reveals that different sampling position would results in completely different shape of velocity profiles. 761 (a) Medium polyhedral mesh (b) Hexahedral mesh (c) Hexacore mesh Figure 5: Contours of averaged velocity magnitude at l/d = 3 5. Conclusions In this work, swirling pipe flow was investigated using advanced scale-adaptive model (𝑘 − 𝜔 SST SAS). The behavior of SAS model was investigated on several meshes including structured and unstructured hexahedral and polyhedral grid. Strong influence of computational grid on the fluid motion was found. The model predicts strong asymmetry, whereas the experiment shows symmetrical flow field. No attempt to resolve the boundary layer grid was performed, as the focal point of this analysis did not include precise prediction of near-wall behavior. This might have led to slight improvement of the results, but could hardly change the character of the turbulent flow core. Acknowledgements The authors gratefully acknowledge financial support within the project NETME CENTRE PLUS (LO1202) which was co-funded by the Ministry of Education, Youth and Sports within the support programme „National Sustainability Programme I”. References Davidson L., 2016, Fluid mechanics, turbulent flow and turbulence modeling, accessed 19.04.2016. Hajek J., Jurena T., Baier B., Bagherzadeh A., Schnick T., Smejkal Q., 2015, CFD Process Modeling of the Industrial Calcination of Diatomite, Chemie Ingenieur Technik, 87(7), 991–997. Menter F.R., 2012, Best practice: Scale-resolving simulations in ANSYS CFD, accessed 02.09.2015. Mitrofanova O.V., 2003, Hydrodynamics and Heat Transfer in Swirling Flows in Channels with Swirlers (Analytical Review), High Temperature, 41(4), 518–559. Papoulias D., Lob S., 2015, Advances in CFD Modelling of Multiphase Flows in Cyclone Separators, accessed 05.09.2016. Pashtrapanska M., Jovanović J., Lienhart H., Durst F., 2006, Turbulence measurements in a swirling pipe flow, Experiments in Fluids, 41(5), 813–827. Rocklage-Marliani G., Schmidts M., Ram V.I.V., 2003, Three-Dimensional Laser-Doppler Velocimeter Measurements in Swirling Turbulent Pipe Flow, Flow, Turbulence and Combustion, 70(1–4), 43–67. Shim Y. M., Sharma, R. N., Richards P. J., 2009, Numerical study of the flow over a circular cylinder in the near wake at Reynolds number 3900, in 39th AIAA Fluid Dynamics Conference, 4160-4173. accessed 14.04.2016. Vaidya H.A., Ertunç Ö., Genç B., Beyer F., Köksoy Ç., Delgado A., 2011, Numerical simulations of swirling pipe flows- decay of swirl and occurrence of vortex structures, Journal of Physics: Conference Series, 318(6), 62022. Wilcox D.C., 2006, Turbulence Modeling for CFD, 3rd ed., DCW Industries, Inc., San Diego, CA, USA. 762