Comparison of turbulent models in the case of a constricted tube

The validation of a proper solution is an indispensable phase of every numerical simulation. Nowadays, many turbulent models are available, whose application leads to slightly different solution of flow behaviour depending on the boundary conditions of a specific problem. It is essential to select the proper turbulence model appropriate for the given situation. The aim of this study is to select the most suitable two-equation eddy-viscosity model, which can be further used during calculations of airflow in human airways. For this purpose, geometry of a constricted tube with well-documented experimental measurements was chosen. The flow in the constricted tube was calculated using Spallart-Almaras, komega, k-epsilon and SST model approach using commercial software. The outcome of the comparison is a choice of the suitable model which is capable of simulating the transition of the boundary layer from laminar to turbulent flow. This transition typically arises in the upper part of the respiratory system, where the airways are constricted, specifically in the area, where the oral cavity continues through the glottis to trachea. The simulations were performed in a commercial solver Star-CCM+.


Introduction
The human airways comprise extremely complicated system of branches with many variable parameters.The route for the oral inhalation of pharmaceutical aerosols consists of the oral cavity, larynx, pharynx, glottis, trachea and bronchi, where the flow divides in the system of bifurcations with branches of decreasing dimensions.
A part of airways that serves for delivery of aerosol into the lungs can be modelled as a tube with a constriction in the area of glottis.This constriction significantly influences the downstream flow [1].The most common flow rates that are being used for the breathing simulation are 15 and 30 l/min.The flow in the upper airways is laminar under these conditions, however, in glottis, where the vocal cords are located, is the flow accelerated above the critical Reynolds number due to the constriction, and therefore the transition to turbulent flow occurs.This condition must be reflected in the process of selection of the turbulent model during the numerical simulation.
It is therefore necessary to select such a model, which accurately simulates flow in a circular tube with low Reynolds numbers, but at the same time is capable of precise modelling of the boundary layer transition from laminar to turbulent flow.Previous studies aiming at similar topics stated, that it is necessary to use the models of turbulence, which do not utilize wall functions, because they are not suitable for our type of flow [2].
The aim of this study is to select the most suitable two-equation eddy-viscosity model integrated in the commercial solver Star-CCM+ (CD Adapco), which can be further used during calculations of airflow in human airways.Accordingly, the models supporting the low wall y+ approach were selected based on the recommendations given by [3] for human airway flow simulations.
The results of this study will bring valuable information for the numerical simulations of human airways.According to Longest [4] is the accurate solution of the flow field necessary for calculations of local deposition characteristics and therefore the correct turbulent model is a necessary condition.

Methods
The geometry of a circular tube with 50 % and 75 % constriction was created in a software Rhinoceros (McNeel, Seattle, USA) for the numerical simulation.Several models of turbulence were used and compared on the geometry.

Model geometry
The test geometry comes from the paper of Ahmed and Giddens [5].A smooth pipe with diameter of 50.8 mm is constricted to 50% of its diameter.The length of the constriction is 101.6 mm.The upstream length of the pipe is 4877 mm to allow the velocity profile to fully develop.The downstream length of the pipe is 1829 mm (see Fig. 1).Owing to the use of turbulent models based on low wall y+ approach, close attention was payed to the boundary layer treatment.A near wall layer consisting of 10 layers of prism cells was created.The value of wall y+ was approximately 1 on whole length of used geometry.The computational mesh was made of polyhedral elements and the prism layer contained 13.2 million cells.(see Fig. 2). .

Numerical method
The Star-CCM+ solver contains many turbulence models, whose application depends on the specific requirements of a given problem.Only the models supporting the low wall y+ approach were eligible from the range of available models in our case.Furthermore, just the models suitable for the internal pipe flow with a constriction could have been selected on the basis of their description given in the Star-CCM+ manual [6].
The following turbulence models were chosen: Standard Spalart-Allmaras model is a good choice for applications in which the boundary layers are largely attached and separation is mild if it occurs.The Spalart-Allmaras turbulence models solve a single transport equation that determines the turbulent viscosity.In its standard form, the model is applied without wall functions.According to the formulation of the model, the entire turbulent boundary layer, including the viscous sublayer, can be accurately resolved and the model can be applied on fine meshes.
Standard k-İ model.Generally, k-İ models provide a good compromise between robustness, computational cost and accuracy.The k-İ turbulence model is a twoequation model in which transport equations are solved for the turbulent kinetic energy and its dissipation rate.
In its original form, the k-İ turbulence model was applied with wall functions, but was later modified to use the Low Re approach for resolving the viscous sublayer.The above mentioned models were applied in default settings without changes in their specific constants.Axial velocities (i.e. the velocity component in the direction parallel with the axis of the tube with the stenosis) were observed.The flowing media in the tube was 63 % solution of glycerol in water at 33 °C, while the Reynolds number was Re = 2000.

Results and discussion
The observed values were compared in the nondimensional coordinate Z = x/D, where x is the axial distance from the throat of the stenosis and D is the diameter of the unconstricted tube (see Fig. 3).The calculated values of the axial velocities for each turbulence model were normalized according to [5] and their values were compared.The figures 5 and 6 show the comparison of velocity profiles from the tube axis (0) to the wall (1) in the distances 0, 1 and 2 diameters (Fig. 5) and in the distances 3, 4 and 5 diameters from the stenosis throat (Fig. 6).Both the graphs illustrate gradual transition of the velocity profile from almost flat profile (Z = 0) to less turbulent profile.Fairly good agreement of all the models on the velocity profile was found in the location Z = 0, i.e. exactly in the stenosis throat.However, despite the fact, that all the models predicted the same velocity profile, the actual measured velocity profile was slightly different (Fig. 5, Z = 0).This small discrepancy was reported also by Tan [3].The velocity profile changes close to the wall due to the recirculation zone in the location Z = 1.The normalized velocity turns into negative values there.This phenomenon was captured by Spalart-Allmaras turbulence model, whose prediction corresponds most with the experiments.All other models also predict recirculation zone, however they significantly underpredict axial velocity values.
The Spalart-Allmaras model agrees well with the experiment also in the following location in the distance of 2 diameters from the throat, nevertheless, all the other models approximate the experimentally measured velocity profile as well.The main difference is that both the Spalart-Allmaras model and experiment show positive values of the normalized velocity.It means that Velocity profiles in distances 3 to 5 diameters downstream of the stenosis throat are plotted in Fig. 6.The graph confirms that k-Ȧ (KO), V2F, SKELR and SST models indeed predict larger recirculation region length compared to the actually measured values.They all predict the end of the recirculation zone in the Z = 5 distance.After that the transition to the ideal turbulent profile can be seen.Spalart-Allmaras model in the locations Z = 4 and 5 shows gradual transition of the velocity profile to the shape qualitatively similar to the experimental values, however, the axial velocity values are significantly underpredicted.

Conclusions
The results show that the Spalart-Allmaras model can well predicts the behaviour in near-wall region, but it is unable to predict flow in axis of pipe.This model was able to capture the separation of boundary layer with subsequent recirculation downstream of the constriction.Best models for predicting of flow in axis of pipe were models based on k-Ȧ approach.The least suitable models for this case can regarded a models based on k-İ approach, which are mostly applied to cases focused on free.Since the bronchial tree, which is part of the respiratory tract consists of a series of consecutive similar cases it is advisable to choose the model that best predicts this behaviour, thus preventing a significant increase of computational errors.

V2F
Low-Reynolds Number k-İ model is known to capture the near-wall turbulence effects more accurately, which is crucial for the accurate prediction of flow separation.This model solves two more turbulence quantities, namely the normal stress function and the elliptic function, in addition to turbulent kinetic energy and dissipation rate.This model is designed to handle wall effects in turbulent boundary layers and to accommodate non-local effects.Standard k-Ȧ model is similar to k-İ models in that two transport equations are solved, but differ in the choice of the second transported turbulence variable.The transport equations that are solved are for the turbulent kinetic energy and a specific dissipation rate.One reported advantage of the k-Ȧ model over the k-İ model is its improved performance for boundary layers under adverse pressure gradients.Perhaps the most significant advantage, however, is that it may be applied throughout the boundary layer, including the viscous-dominated region, without further modification.SST k-Ȧ model.The problem of sensitivity to freestream/inlet conditions was addressed by Menter[7], who recognized that the İ transport equation from the standard k-İ model could be transformed into an Ȧ transport equation by variable substitution.Menter suggested using a blending function that would include the cross-diffusion term far from walls, but not near the wall.This approach effectively blends a k-İ model in the far-field with a k-Ȧ model near the wall.In addition, Menter also introduced a modification to the linear constitutive equation and dubbed the model containing this modification the SST (shear-stress transport) k-Ȧ model.The SST model has fairly wide application in situations, where viscous flows are typically resolved and turbulence models are applied throughout the boundary layer.

Fig. 3 .
Fig. 3. Cross-sections used for comparison The results of the comparison of velocity profiles in the axis of the stenosis can be found in Fig. 4. It is obvious from the graph, that different turbulent models estimate the flow behavior with various accuracy.The most reliable models are k-Ȧ SST model (SST) and V2F Low-Reynolds Number k-İ model (V2F), which give the best agreement with experiment in terms of the axial velocity, as documented in the Fig. 4. The worst agreement was found for Spallart-Allmaras (SA) and Standard k-İ model, which predict faster decrease of velocity downstream of the stenosis in the distance 2 -3 diameters from the stenosis throat.The wrong prediction was probably caused by an inaccurate estimate of the shear layer separation.The figures 5 and 6 show the comparison of velocity profiles from the tube axis (0) to the wall (1) in the distances 0, 1 and 2 diameters (Fig.5) and in the distances 3, 4 and 5 diameters from the stenosis throat (Fig.6).Both the graphs illustrate gradual transition of the velocity profile from almost flat profile (Z = 0) to less turbulent profile.Fairly good agreement of all the models on the velocity profile was found in the location Z = 0, i.e. exactly in the stenosis throat.However, despite the fact, that all the models predicted the same