Study of Airflow Behavior for Duplex Circular Cylinders

The modeling of atmospheric ice accretion on duplex cylinders received a limited attention, with modeling carried by Wagner and Qing et al. The publicly available experimental data about the ice accretion on the duplex cylinders is limited to experiments of Qing et al. and Veerakumar et al. When comparing the data of Wagner and Qing et al. with the results of Veerakumar et al., the major difference is the airflow behavior in the wake of the windward cylinder, the extent of the wake and recirculation bubble, and the velocity distribution in the wake. Thus, its needed to study the effect of the turbulence model on the airflow behavior of duplex cylinders, with focus being the behavior of the wake of the windward cylinder. This study reports the simulation results of the complex airflow behavior of duplex circular cylinder bundle obtained using several turbulence models employed by commercial CFD code.


INTRODUCTION
The study of atmospheric icing on structures is a well-established field consisting of several decades of accumulated knowledge.For example, the analytical models in use for the power line icing specifically originate as early as 1980's [1].These models, in turn, are based on the theoretical work on atmospheric icing of structures, dating back to the works of [2] and other research, conducted at the Mt.Washington Observatory in the same timeframe.The resultant aggregated theoretical knowledge has been incorporated in the ISO 12494 standard "Atmospheric Icing of Structures" [3].ISO 12494 modeling framework received widespread attention when it comes to the analytical modeling of ice accretion on simple geometries, which can be approximated by circular cylinder, such as simplex power line conductors, tubular telecommunication masts, etc. [3].
Consequently, the modeling of atmospheric icing on structures with other geometrical configurations has received limited attention.Partially, this due to the analytical model of Finstad [4] and Finstad et al. [5], which is a basis of the ISO 12494 modeling framework [3] only being applicable to the simplex circular geometries.However, the modeling of other geometries for the purposes of studying the atmospheric icing process has major practical applications, with one practical example of such a case is modeling of atmospheric icing on bundled conductors -duplex, triplex, hexa, etc. bundled configurations.Such geometric configurations are of significant importance in the modeling of atmospheric icing on the overhead transmission lines, as the majority of high-voltage transmission networks consist of bundled conductors.However, to the best of authors' knowledge only limited amount of work has been done in modeling of atmospheric ice accretion on bundled cylinders and/or conductors, those primarily being the works of Wagner [6], Qing et al. [7] and Veerakumar et al. [8].
Thus, there is a need to better understand the process of ice accretion, in both the airflow behavior and droplet impingement, although, it is not certain how or if ISO 12494 can handle the modeling of the atmospheric ice accretion on bundled conductors, nor how one could easily validate the obtained results.
Therefore, the use of Computational Fluid Dynamics (CFD) solvers is employed, which have been steadily increasing in popularity for the purposes of modeling atmospheric ice accretion [9], [10].The results of CFD modeling of atmospheric icing on structures have been extensively validated, primarily in the in-flight icing studies, for example, works by Papadakis et al. [11], Ratvasky et al. [12], Wiberg et al., [13], and more recently -Sokolov and Virk [14], etc.However, when it comes to the geometric configurations different from the either standard circular collector or simple airfoil, there is an open question and a knowledge gap as to how well the existing CFD models can recover the airflow solution, needed to study the droplet impingement afterwards in the icing model.The focus of this work is to ascertain the applicability and "correctness" of some widely used numerical turbulence models, employed in CFD code, for the flow around in-line duplex circular cylinders.
The setup of the numerical simulations, performed in this study, is aiming to replicate the experimental conditions in Veerakumar et al. [8] who have carried out as series of icing wind tunnel experiments on bundled conductors, corresponding to the typically used power line conductors, under typical icing conditions.For simplicity, in this study, the usage of circular cylinders is employed, and they are deemed an acceptable approximation to the actual stranded conductor with ribbed bare surface, especially after initial accretion smoothens the conductor surface.The details and setup of numerical simulations are given in subsequent section.

DESIGN OF EXPERIMENT
Table 1 shows the operating conditions in this study with Figure 1 giving a schematic overview of the duplex bundled cylinders setup.The choice of cylinder diameter, bundle separation and the operating follows those of the Veerakumar et al. experiments [8].This is done in order to permit direct comparison of the experimental and the numerical results.Furthermore, in order to streamline the subsequent analysis some assumptions and simplifications have been made in this study.First, the choice of circular cylinder over stranded conductor is deemed an acceptable simplification, based on the previous results and experiences in the field of atmospheric icing.This extends to the assumption of airflow behavior also.Second, the comparison of drag coefficient and drag forces ratio is done for the reported values in the Veerakumar et al. [8], before the spray nozzles are turned on in the experiments (corresponding to the experimental duration of less than 100 s).This is to exclude the influence of the accreted ice shape on the values of drag coefficient.Finally, the velocity contours are compared with the time-averaged experimental velocity contours.This is due to the fact that Veerakumar et al. do not present additional Particle Imaging Velocimetry (PIV) contours in their work.

Numerical Setup
The numerical CFD simulations were carried out using ANSYS Fluent, which is based on the Eulerian-Lagrangian approach in which the Eulerian approach is used for the analysis of the carrier phase dynamics while the Lagrangian approach is used for the analysis of the dispersed phase.The Eulerian method treats the carrier phase as a continuum and develops its conservation equations on a control volume basis and in a similar form as that for the fluid phase.The Lagrangian method considers particles as a discrete phase and tracks the pathway of each individual particle.By studying the statistics of particle trajectories, the Lagrangian method is also able to calculate the particle concentration and other phase data.On the other hand, by studying particle velocity vectors and its magnitudes in Eulerian method, it is possible to reconstruct the pathways and trajectories of particles in a phase.
The equation for conservation of mass, or continuity equation, can be written as follows [10]: Equation ( 1) is the general form of the mass conservation equation and is valid for incompressible as well as compressible flows.The source is the mass added to the continuous phase from the dispersed second phase (for example, due to vaporization of liquid droplets) and any user-defined sources.
Conservation of momentum in an inertial (non-accelerating) reference frame is described by [15]: where  is the static pressure,  is the stress tensor (described below), and  and  are the gravitational body force and external body forces (for example, that arise from interaction with the dispersed phase), respectively. also contains other model-dependent source terms such as porous-media and user-defined sources.
The stress tensor  is given by [15]: where μ is the molecular viscosity,  is the unit tensor, and the second term on the right-hand side is the effect of volume dilation.Ansys Fluent is a RANS-based (Reynolds-Averaged Navier-Stokes) solver.In Reynolds averaging, the solution variables in the instantaneous (exact) Navier-Stokes equations are decomposed into the mean (ensemble-averaged or time-averaged) and fluctuating components.For the velocity components [15]: where u � i and u i ' are mean and fluctuating velocity components ( = 1, 2, 3).Likewise, for pressure and other scalar quantities: where  denotes a scalar such as pressure, energy, or species concentration.Substituting expressions of this form for the flow variables into the instantaneous continuity and momentum equations and taking a time (or ensemble) average (and dropping the overbar on the mean velocity, ū) yields the ensemble-averaged momentum equations.They can be written in Cartesian tensor form as [15]: ∂ ∂t Equations ( 6) and ( 7) are called Reynolds-averaged Navier-Stokes (RANS) equations.They have the same general form as the instantaneous Navier-Stokes equations, with the velocities and other solution variables now representing ensemble-averaged (or timeaveraged) values.
Additional terms now appear that represent the effects of turbulence.These Reynolds stresses, -ρu i ' u j ' ����� , must be modeled in order to close Equation (7).
The Reynolds-averaged approach to turbulence modeling requires that the Reynolds stresses in Equation ( 7) are appropriately modeled.A common method employs the Boussinesq hypothesis to relate the Reynolds stresses to the mean velocity gradients [15], [16]: The Boussinesq hypothesis is used in the Spalart-Allmaras model [17], the k-ε models [18], and the k-ω models [19], [20], [21].The advantage of this approach is the relatively low computational cost associated with the computation of the turbulent viscosity, μ t .In the case of the Spalart-Allmaras model, only one additional transport equation (representing turbulent viscosity) is solved.In the case of the k-ε models, and the k-ω models, two additional transport equations (for the turbulence kinetic energy, k, and either the turbulence dissipation rate, ε, or the specific dissipation rate, ω) are solved, and μ t is computed as a function of k and ε or k and ω.The disadvantage of the Boussinesq hypothesis as presented is that it assumes is an isotropic scalar quantity, which is not strictly true.However, the assumption of an isotropic turbulent viscosity typically works well for shear flows dominated by only one of the turbulent shear stresses.In many cases, models based on the Boussinesq hypothesis perform very well, and the additional computational expense of the Reynolds stress model (RSM) is not justified [15].
The computational mesh used for all CFD simulations is a hybrid mesh, consisting of structured quad elements near the cylinders and the unstructured tri elements elsewhere.The first cell height at both cylinders is 1×10 -6 m with exponential growth factor of 1.1 and a total of 100 inflation layers.In addition, the length of the cylinder "wall" itself is divided into 100 nodes.This results in 30,000 structured quad elements per cylinder (as a three-cell extrusion in z-direction is used), for a total of 60,000 structured quad elements.The rest of the mesh is filled with unstructured tri mesh, for a total cell count of 105,960 cells in the computational domain.Furthermore, fine quad and tri elements were used in cylinder wake region, whereas coarse tri elements were used for the rest of the outer domain.The computational mesh is also shown in Figure 2. Detailed mesh sensitivity analysis was carried out to accurately determine the boundary layer characteristics (shear stress and heat fluxes), a y+ values of less than 1 is used near the cylinder wall surface.Number of mesh elements and y+ value was selected based upon the heat flux calculations, where a numerical check was imposed that the heat flux computed with the classical formulae / should be comparable with the heat flux computed with the Gresho's method.
The calculation of  + value is performed in the following way [22]: where ρ f and μ f are the density and dynamic viscosity of the continuous phase (air), U ∞ is the freestream velocity, L ∞ is the characteristic length, i.e., cylinder diameter, C f is the skin friction coefficient, τ wall is the shear stress at the wall, U fric is the friction velocity and Δs is the wall spacing (first cell height).These computations are based on the flat-plate boundary layer theory from [22].Based on the highest operating wind speed in the Table 2 being equal to 20 m/s, this gives the  = 12.2 × 10 -6 m for + = 1.Thus, the first cell height used in the meshes in this study, and being equal to  = 1 × 10 -6 m yields a  + value of + = 0.08 for the highest freestream wind speed value of 20 m/s.
The turbulence models chosen for this study are the k-ω Shear Stress Transport (SST) model [19], the Transition SST (also known as The Langtry-Menter 4-equation Transitional SST Model [20], [21]), the Spalart-Allmaras [17], and the Realizable k-ε models, [18] [15].This choice of these models was governed by a few factors.First, in the authors opinion, these models are commonly used for this particular type of modeling, when considering other available turbulence models in Ansys Fluent.Second, is to test the hypothesis of Wagner [6] of the potential overestimate of wake extension of k-ω models and the potential implication of this on the resultant airflow behavior.Third, is to test the performance of the k-ω SST model itself for this type of modeling, as the k-ω SST is a widely used turbulence model, which combines the robustness of k-ω model in near-wall and boundary layer region with the reliability of k-ε model in the far field region.Last, is to ascertain, if the modeling of the bundled conductors can be carried within the constraints of one commercial CFD package, without coupling the solution procedure to other CFD packages, e.g., Fluent or the in-house code.
All simulations were carried out in Ansys Fluent in transient (unsteady) mode.The time step chosen was one millisecond, based on the analytical calculations using the Strouhal number for circular cylinder (St = 0.20).Finally, the total simulation duration is 15,000 timesteps, which corresponds to the total flow time of 15 s.This value is deemed sufficient to achieve the solution convergence and flow reaching stead-state, following similar procedure to one in Sokolov and Virk [14].

RESULTS AND DISCUSSION
The comparison of results will be carried out using a few different parameters -the velocity magnitudes, and their distributions; Root Mean Square Errors (RMSEs) of velocity magnitudes, and their distributions; the pressure coefficient, C  , values, and their distributions; the drag coefficient, C  ; and the drag coefficient ratios.From the results in Figure 3 a few observations can be made.First, is that the Realizable k-ε and Transition SST models predict markedly different airflow behavior in the wake of the windward cylinder when compared to the k-ω SST and Spalart-Allmaras models.The Realizable k-ε model predicts the shortest extent of the windward cylinder's wake and recirculation zone, while the Transition SST model predicts them to be the largest out of all models.The k-ω SST and Spalart-Allmaras models have a decent agreement between each other as to the extent of these effects, and the length of the windward wake is roughly "inbetween" the results of the Realizable k-ε and Transition SST models.Furthermore, it is of interest to note that the hypothesis of Wagner [6] of the potential overestimate of wake extension of k-ω models and the potential implication of this on the resultant airflow behavior holds with these results.
Second, is the values of the velocity magnitudes and the velocity distributions.The Realizable k-ε model predicts the highest velocity magnitudes and the sharpest velocity gradients, while the results from the Transition SST model have opposite trends.Again, the obtained results from the k-ω SST and Spalart-Allmaras models are "in the middle" though the agreement between them is poorer than in the case of wake behavior, with the Spalart-Allmaras model having higher velocity magnitude values and sharper velocity gradients than the k-ω SST model.
In order to better understand the differences in the wake behavior, velocity magnitudes and their distributions, Figure 4 shows the values of RMSE for all tested models in this study.4 shows a few noteworthy observations.First, is that all tested models show a high value of RMSE in the wake of the leeward cylinder.This indicates that there is a large uncertainty in the results and the flow behavior in the wake of the leeward cylinder.Second, is that all models, with the exception of the Realizable k-ε model feature the almost the same RMSE values and their distributions, indicating that this behavior is mostly consistent across the tested models.The primary hypothesis as to this is the extent and the influence of wake of the windward cylinder on the results of the leeward cylinder.Last, is that the behavior Realizable k-ε model is notably different from the rest of the models -the maximum RMSE value in the wake of leeward cylinder is approximately half, when compared to the other models, however, the Realizable k-ε model also features high RMSE values in the wake of the windward cylinder, indicating an uncertainty in the results for that wake also (for the remaining models the RMSE in the wake of the windward cylinder is approximately zero).
The combined results from Figures 3 and 4 suggest that such differences in the airflow behavior may be caused due to inherent differences in the treatment of the production and dissipation of the Turbulent Kinetic Energy (TKE) and the momentum conservation across the different models, however, such a detailed analysis is deemed to be outside the scope of the current paper.When it comes to comparison with the results of PIV measurements of Veeramuar et al. [8], Figure 5 shows the velocity distributions, as measured in the experiments.
Figure 5. Experimental values of the PIV measured velocity distributions in Veerakumar et al [8].
When comparing the experimental results from Figure 5 with the results from the numerical simulations from Figure 3 a few observations can be made.First, all tested models do overestimate the length of the recirculation bubble behind the leeward cylinder, and all, bar the Realizable k-ε model overestimate the wake length and the recirculation zone of the windward cylinder also.However, the Realizable k-ε model erroneously recovers the velocity fully in the wake of windward cylinder.The Realizable k-ε model shows the worst agreement with the experimental measurements, when it comes to the velocity magnitudes, generally featuring considerably higher velocity magnitudes than those measured in the experiments, along with being the only model with large RMSE in the wake of windward cylinder.
The main issue of the Transition SST model is the larger extent of the wake behind the windward cylinder, which significantly affects the flow around the leeward cylinder, although, this model has the closest agreement with the experimental values when it comes to the values of the velocity magnitudes (however, all models overestimate them, to an extent).Overall, from the tested turbulence models, the standard k-ω SST model shows the closest agreement with the experimental measurements, followed by the Spalart-Almaras model.
Before proceeding with the results and discussion of the obtained drag coefficients results, the pressure coefficient and its distribution have to be looked at, in order to establish the baseline for comparison.Figure 6 shows these values, as obtained from the results of CFD simulations.These results suggest that the drag coefficient values, and their ratios will be quite different for all tested models.These values are shown in the Figures 7 and 8, respectively, while Table 2 shows some aggregated values of interest from the preceding discussion, for the readers convenience.From the values of Figure 7 and Table 3 it can be seen that the Spalart-Allmaras model predicts the closest value of the drag coefficient compared to the reference value (this being C  = 1.17 for the windward cylinder) while k-ω SST overestimates it slightly.On the other hand, both the Realizable k-ε and the Transition SST models significantly under-and overestimate it, respectively.Compared to the drag coefficient ratio between the leeward and the cylinder, as reported by Veerakumar et al. [8], and it being approximately equal to 0.65, all tested models, with exception of the Transition SST model do overestimate this ratio.The primary reason behind this is believed to be high fluctuation in the drag coefficient, as evidenced by the Figure 7, with these fluctuations believed to be the function of the high RMSE in the wake of leeward cylinder, as shown in Figure 4.The obtained results show that all tested models do overestimate the length of the recirculation bubble behind the leeward cylinder.Moreover, all tested models do have a high RMSE value for the velocity in the wake of leeward cylinder, with the Realizable k-ε model featuring a large RMSE in the wake of windward cylinder.Moreover, the Realizable k-ε model erroneously recovers the velocity fully in the wake of windward conductor.In addition, the Realizable k-ε model has the highest velocity magnitudes recorded, in excess of other models and experimental measurements.
Furthermore, all models tend to predict largely the same positive values of C  and their distributions, around the windward cylinder, with significant spread in the negative C  values with the with the Realizable k-ε model having the highest values and the Transition SST model having the lowest, around both the windward and the leeward cylinder.Consequently, this also extends to the values of the drag coefficients and their ratios, with the Realizable k-ε model having significant underestimation of the windward drag coefficient with overestimation of the ratio, while the Transition SST model features an opposite situation.
Summarizing, from the tested turbulence models, the standard k-ω SST model showed the closest agreement with the experimental measurements, followed closely by Spalart-Allmaras model.The Realizable k-ε model shows the worst agreement with the experimental measurements, along with being the only model with large RMSE in the wake of windward cylinder.Finally, the main issue of the transition SST model is the larger extent of the wake behind the windward cylinder, which significantly affects the flow around the leeward cylinder.

Figure 2 .
Figure 2. Computational mesh.Bottom image shows the closeup of the windward cylinder wake region.

Figure 3 .
Figure 3. Velocity magnitude distributions in the CFD simulations for different models.From top to bottom: k-ω SST, Spalart-Allmaras, Realizable k-ε and Transition SST.

Figure 4 .
Figure 4. RMSE of velocity magnitude in the CFD simulations for different models.From top to bottom: k-ω SST, Spalart-Allmaras, Realizable k-ε and Transition SST.

Figure
Figure4shows a few noteworthy observations.First, is that all tested models show a high value of RMSE in the wake of the leeward cylinder.This indicates that there is a large uncertainty in the results and the flow behavior in the wake of the leeward cylinder.Second, is that all models, with the exception of the Realizable k-ε model feature the almost the same RMSE values and their distributions, indicating that this behavior is mostly consistent across the tested models.The primary hypothesis as to this is the extent and the influence of wake of the windward cylinder on the results of the leeward cylinder.Last, is that the behavior Realizable k-ε model is notably different from the rest of the models -the maximum RMSE value in the wake of leeward cylinder is approximately half, when compared to the other models, however, the Realizable k-ε model also features high RMSE values in the wake of the windward cylinder, indicating an uncertainty in the results for that wake also (for the remaining models the RMSE in the wake of the windward cylinder is approximately zero).The combined results from Figures3 and 4suggest that such differences in the airflow behavior may be caused due to inherent differences in the treatment of the production and dissipation of the Turbulent Kinetic Energy (TKE) and the momentum conservation across the different models, however, such a detailed analysis is deemed to be outside the scope of the current paper.When it comes to comparison with the results of PIV measurements of Veeramuar et al.[8], Figure5shows the velocity distributions, as measured in the experiments.

Figure 6 .
Figure 6.Pressure coefficient values in the CFD simulations for different models.From top to bottom: k-ω SST, Spalart-Allmaras, Realizable k-ε and Transition SST.

Figure 6
Figure6shows an interesting behavior of both SST models tested in this work.For both of them, the influence of the front wake on the leeward cylinder is strong enough to prevent the formation of clearly defined stagnation point.Furthermore, they do feature larger extent of the distribution of negative pressure coefficient values, although the negative   values are smaller than those for the Realizable k-ε and Transition SST models.Furthermore, the extent of the suction side of the windward cylinder for the Transition SST model is large enough to call the physicality of these results into question.Moreover, unlike with results for the velocity distributions and velocity RMSE, the pressure coefficient values for different models display the large spread in the negative   values from -0.91 (in the results of the Transition SST model) to the -2.01 for the Realizable k-ε model, although the maximum positive C  values do not differ by much (1.01 for the Realizable k-ε model, and 1.02 for the rest of the models).

Table 2 .
Aggregated simulation data.Figure8shows the drag coefficient ratios of bundled cylinders.This ratio is obtained by dividing the C  of leeward cylinder by C  value of windward cylinder.Moreover, the aggregated values from Figures7 and 8are shown in the table form in Table3.

Table 3 .
Values of drag coefficients in the CFD simulations