Numerical investigations of the vortex feature-based vorticity confinement models for the assessment in three-dimensional vortex-dominated flows

In order to improve the vortex resolution in aerodynamic wakes, a locally normalized vortex feature-based vorticity confinement method is implemented into the multi-block, structured computational fluid dynamics solver (ROSITA). In this method, the second vorticity confinement (VC2) scheme with two well-known vortex feature detection methods (non-dimensional Q criterion, non-dimensional λ2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _2$$\end{document} criterion) is employed to counterbalance the truncation error introduced by the numerical discretization of the convective term. The flow field of two benchmark three-dimensional steady vortex-dominated cases, the NACA0015 wing and the Caradonna–Tung hovering rotor, is simulated with the implemented method. The improvements in aerodynamics prediction, vorticity preservation, computational stability, and efficiency are demonstrated. From the numerical results, the vortex feature-based confinement models significantly improve the computational stability, the aerodynamic loads prediction and vorticity preservation capability, especially for the λ2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda _2$$\end{document}–based VC2 model. In addition, it allows the use of higher confinement parameters on a coarse grid with a relatively higher computational efficiency to obtain better results than those of a finer grid.


Introduction
In vortex-dominated flows, like the airplane wake and helicopter rotor wake, the numerical solution of vortices tends to dissipate at a higher rate than what is observed in real life due to the diffusion and dissipation from the Navier-Stokes equations discretization. This significantly affects the evaluation of the aerodynamic performance. Over the past decades, numerous numerical methods have been proposed to help reduce vortices dissipation and increase the capability of vorticity preservation, including the high fidelity methods (LES [1], IDDES [2], ⋯ ), high-order flux reconstruction schemes (5th-order WENO [3], 4th-order MUSCL [4], ⋯ ), and mesh adaptation algorithms (AMR [5]). Despite the improved accuracy of these approaches, such techniques are computationally expensive and increase the complexity of numerical codes. In contrast, the vorticity confinement (VC) methods, which allow to conserve vortices over a Abstract In order to improve the vortex resolution in aerodynamic wakes, a locally normalized vortex feature-based vorticity confinement method is implemented into the multi-block, structured computational fluid dynamics solver (ROSITA). In this method, the second vorticity confinement (VC2) scheme with two well-known vortex feature detection methods (nondimensional Q criterion, non-dimensional 2 criterion) is employed to counterbalance the truncation error introduced by the numerical discretization of the convective term. The flow field of two benchmark three-dimensional steady vortex-dominated cases, the NACA0015 wing and the Caradonna-Tung hovering rotor, is simulated with the implemented method. The improvements in aerodynamics prediction, vorticity preservation, computational stability, and efficiency are demonstrated. From the numerical results, the vortex feature-based confinement models significantly improve the computational stability, the aerodynamic long period of time by adding an anti-diffusion source term into the Navier-Stokes equations, present a remedy with a moderate computational cost. For this reason, the use of VC methods in support of computational studies for the vortex-dominated flows has led to an interesting topic in the last few decades.
Since Steinhoff and co-workers proposed the original VC (VC1) method in their works [6][7][8], many studies have been carried out on this topic. In the early research, several attempts have been carried out to extend this method to compressible flows [9][10][11]. After G.Hu et al. [11] presented a more stable VC scheme for the compressible flows by interpreting Steinhoff's formulation as a body force, the primary research on the VC method changed to the development of high-order schemes and the determination of the confinement parameter .
In terms of the high-order VC schemes, Costes et al. proposed a series of high-order VC formulations based on the second vorticity confinement (VC2) scheme and applied them to the rotor wake simulations [12][13][14][15]. M. Costes and G. Kowani also derived a dynamic confinement parameter which is related to the local vorticity in Ref. [16]. Robinson [17] developed an expression of the confinement parameter based on helicity. Hahn and Iaccarino [18] introduced a new adaptive VC parameter related to the difference between central and upwind discretisation of the convection terms. In more recent research, the VC method is coupled with the TVD (Total Variation Diminishing) limiters to reduce the sensitivity of the solutions to the value [19].
Although many studies on the VC approach were carried out in the past decades, there still remains an intrinsic issue. In most of these studies, the non-zero vorticity magnitude is used as the factor for vortex identification. It means that the vorticity confinement term is computed at each point in solution domain where the vorticity magnitude is not equal to zero. However, using vorticity magnitude to define a vortex structure is not adequate. It will not provide the correct results in the areas where the vorticity magnitude is non-zero but there is no vortex, like in the boundary layer, for instance. To eliminate this problem, R. Boisard and co-workers [20] used the Q criteria to avoid the application of confinement inside the boundary layer. Feder et al. [21] employed the 2 criterion as a limiter to restrict the vortices region in tracking the tip vortex of NACA0012 wing. Mohseni [22] introduced a class of hybrid methods that combine four different vortex feature detection approaches with the VC1 method and investigated the effects of the hybrid techniques on the computational precision and speed with a 2D stationary single vortex. However, the investigation of the above hybrid VC methods on three-dimensional flows are still ignored. In addition, the performance of different vortex feature detection VC approaches need to be evaluated. This paper presents two vortex feature-based vorticity confinement models based on the VC2 scheme and the well-known vortex detection formulations, Q-criterion and 2 -criterion, for the simulation of three-dimensional vortex-dominated flows. Both of these two vortex detection methods are normalized by the local shear-strain rate. At first, the benchmark NACA0015 wing case is calculated; these two formulations are compared with the original VC2 model regarding the aerodynamics prediction, vortex profile and computational stability. Then, the influence of each VC model on typical rotational flows is investigated by calculating the case of the Caradonna-Tung rotor in hovering flight mode. The original VC2 method and the two vortex feature-based vorticity confinement methods are briefly referred as OVC2, FVC2-Q, and FVC2-L2.

Flow solver (ROSITA)
This work is performed using a multi-block, structured Reynolds-Averaged Navier-Stokes (RANS) solver developed at Politecnico di Milano, ROSITA (ROtorcraft Software ITAly) [23]. ROSITA is a threedimensional, MPI-parallel, finite-volume solver with the one-equation Spalart-Allmaras turbulence model and CHIMERA technique. The Navier-Stokes equations are formulated in terms of the absolute velocity, expressed in a relative frame of reference RF linked to each component grid to simplify the flow field solution in the overset grid system. The finite-volume forms could be written as: xx n x + yx n y + zx n z xy n x + yy n y + zy n z xz n x + yz n y + zz n z Φ x n x + Φ y n y + Φ z n z + K T∕ n donates the vector of conservative variables inside the flow domain. The expressions, f c and f d , represent the convective flux tensor and diffusive flux tensor, respectively. n is the outward normal unit vector and v is the entrainment velocity vector. V ijk is the cell volume and S ijk is the cell surface, f s stands for the source term due to the movement of the relative reference frame. Ω x , Ω y , Ω z are the angular velocity components of the RF in the absolute frame.
In ROSITA, the space discretization leads to a system of ordinary differential equations for the rate of change of the conservative flow variables associated with the centers of the cell volumes, the equation (1) then reads: where R ijk stands for the flux balance across surface S ijk of the hexahedral cell (i, j, k). The flux balance can be written as: where (Q c ) ijk is the convective flux balance (convective and pressure effects), ( d ) ijk is the diffusive flux balance (viscous effects), (F s ) ijk is the source terms. The convective flux is approximated by the use of Roe-MUSCL 2nd-order discretization with a modified version of Van Albada limiter [24], and a standard second-order central discretization calculates the diffusive flux; the components of stress tensor are discretized by the application of the Gauss theorem.
An implicit dual-time method is used for the time derivative; the equation (2) can be replaced by an implicit 2nd-order backward differential formula: where the state vector W n+1 ijk is solved by sub-iterations in pseudo-time at each physical time step Δt . In the sub-iteration steps, a generalized conjugate gradient (GCG) [25], in conjunction with a block incomplete lower-upper (BILU) pre-conditioner [26], is implemented.
The CHIMERA technique is based on the modified Chesshire and Henshaw algorithm [27]. The Oct-tree and alternating digital tree data structure are adopted to speed up the tagging process.

Vorticity confinement method
In the past, Steinhoff proposed two famous vorticity confinement schemes, called VC1 [6] and VC2 [28], respectively. In this work, the VC2 scheme, which has the advantage of making momentum conservative and not singular at the vortex center, is adopted for the preservation of vorticity in vortex-dominated flows.
The implementation of the VC2 scheme is based on the experience accumulated with the VC1 scheme, that is, the confinement term is added to the momentum equation alone as a body force term f b . This approach is preferred because far better results are obtained when the vorticity confinement term is removed from the energy conservation equation [29].
The expression of f b may be written as: where the quantity is a small positive value to prevent division by 0 in w . The vector w can be interpreted as the locally normalized vorticity vector multiplied by the harmonic mean of vorticity magnitude. The harmonic mean is calculated over a localized stencil of N cells ; for a uniform hexahedral mesh, N = 7, which involves the center cell and six is the confinement parameter which is a positive coefficient. The values used in this work come from a trial and error procedure. The vorticity magnitude is defined as follows: With the expression of f b , the source term f s for the compressible confinement formulation is modified as:

Vortex feature detection method
In this section, two well-known vortex feature detection methods are introduced: Q-criterion, and 2 -criterion methods. Both methods are presented in non-dimensional forms by imposing a normalization with the local shear-strain rate. The vortex feature detection approaches establish a threshold function, f threshold , for identifying the vortex. This function is evaluated at each grid cell, and a vortex is recognized if the f threshold value overtakes a predetermined value.

Non-dimensional Q
Ref. [30] identified the vortex structure in flow regions with a positive second Galilean invariant (Q) of velocity gradient tensor ∇u . Q is expressed as the difference of Frobenius norm between the asymmetric tensor Ω and symmetric tensor S with the following expression: where S and Ω can be respectively written as As the Q value is still dependent on the local characteristic length and velocity, a suitable non-dimensional form is obtained by dividing the equation (9) by ∥ S ∥ 2 F . The resulting threshold function becomes According to the criterion, vortices are identified where the value of the threshold function f threshold is greater than zero.

Non-dimensional 2
The 2 criterion for the identification of vortex was first proposed by Ref. [31]. By neglecting the viscosity and transient terms, the incompressible Navier-Stokes equation turned into an eigenvectoreigenvalue problem as The 2 criterion evaluates the second largest eigenvalue of the symmetric tensor S 2 + Ω 2 , namely 2 . A vortex is then detected with the condition of 2 < 0 . Similar to the non-dimensional Q criterion, ∥ S ∥ 2 F is used to normalize the non-dimensional 2 criterion. The threshold function can be written as where a vortex is obtained with a positive threshold value ( f threshold > 0 ). In this work, since the S 2 + Ω 2 is real and symmetric, the eigenvalues can be readily calculated by the non-iterative algorithm, presented in Ref. [32].

Results and discussion
This section considers two practical test cases to analyze the second vorticity confinement method with different vortex feature detection formulations in detail. First, in the case of the NACA0015 wing, the performance of the vortex feature-based vorticity confinement methods is evaluated in terms of aerodynamics prediction, vorticity preservation, and computational stability. The second part contains the solutions of the helicopter rotor flow with the application of implemented VC models.
The flow around the NACA0015 wing with square tips is now considered to demonstrate the performance of the implemented VC models with the Navier-Stokes equations. Numerical simulations without any VC models are also undertaken for comparison. The determination of the value of the confinement parameter is the critical issue in applying the VC models for a given discretization, a too large constant value may lead to some robustness problems in the simulation. On the other side, a small value may have an insignificant effect on the vorticity confinement. In this work, the values were carefully determined by a trial and error procedure, seeking, for the given discretization, the largest value which allows to obtain a reduction of the residuals of at least three orders of magnitude. The values so obtained were considered as optimal values.
In Ref. [33], an experimental campaign was conducted in the 7 × 10-Foot subsonic wind tunnel at NASA Ames to measure the NACA0015 wing pressure and trailing vortex. The experiment refers to a NACA0015 wing with aspect ratio of AR = 6.6 and chord length of c = 0.52m at different operation conditions. In this work, the selected test case is operated under the condition of = 12 • , M = 0.1235 , and Re = 1.5 × 10 6 . An overset grid system generated using the Chimera technique was adopted for the study of the NACA0015 wing case. It consists of a background grid, a vortex grid, and a near-body wing grid (see Fig. 1 ). The background grid was discretised with an H topology. In stream-wise direction, the inflow section was located at 12c from the leading edge of the wing, the outflow section was placed at a distance of 19c from the trailing edge of the wing. The far-field boundary was extended up to 9c from the wingtip in the span-wise direction and to 7c from the wing surface in the normal direction. An intermediate vortex grid was extended up to 6c from the trailing edge to maintain the integrity of the stream-wise vortex. The near-body wing grid was built with a C-H topology. The wall distance of the first layer of body surfaces was set to 1 × 10 −5 c so that the y + value was less than 1. A non-slip boundary condition was applied on the NACA0015 wing surface. The zipped-grid technique [34] was employed on overlapping surface grids to deal with the wing root configuration with no gap between the wing root section and the wind tunnel wall.
To account for the grid spatial resolution effects on aerodynamic loads prediction, different wing meshes Fig. 1 Computational domain and detailed view of NACA0015 wing grid system: a overset grid system, b cross-section view of the overset grid system on the same geometry and topology were employed with cell densities increasing from 1.33 to 5.96 million. Details of grid discretization information are reported in Table 1 for background and vortex meshes and in Table 2 for wing mesh.

Grid sensitivity study
A grid sensitivity study was conducted first for the wing mesh with three different grid spatial resolutions, namely coarse, medium and fine (see Table 2). The simulations were performed at two test conditions with = 8 • and 12 • to assess the accuracy of the predicted aerodynamic loads. A quantitative analysis of the grid sensitivity of the solutions is reported in Fig. 2. In this figure, C L and C D values are compared with the data measured from wind tunnel [32]. The comparisons show that all computed C L values are slightly lower than the experimental data, small differences between CFD results can be appreciated for the C L values. However, the C D value computed by the coarse grid is significantly different from the results of the medium and fine grids. A weak sensitivity on the grid resolution can be observed for the computed C L and C D values, especially for the medium and fine grids. Therefore, the spatial discretization of the medium grid can be considered for further analysis to guarantee precision and computational efficiency simultaneously.

Effect on computational stability
Regarding the stability of each VC model, Fig. 3 shows the convergence histories of flow solution residual at = 0.002, 0.005 and 0.02, as well as the lift coefficient at = 0.02. The data of the non-VC model case is used for comparison. All flow solutions were simulated using the RANS equations, Table 1 Details of the background and vortex grids (Minimum spacing is outlined in terms of airfoil section chord c) Background 1545480  159  108  90  50  50  52  Vortex  2528800  160  145  109 37.5 7.5 7.5 Table 2 Details of the wing grid ( : chord-wise, : span-wise, : normal, minimum spacing is outlined in terms of airfoil section chord c.)   Fig. 3a-c, it is observed that the residuals of the OVC2 model go up as the value of the confinement parameter increases. It indicates that the stability of the OVC2 model case gradually worsens as the value is increasing. In contrast, the stability of the calculations carried out with two vortex feature-based VC (FVC2-Q and FVC2-L2) models could be better preserved for all tested values, although there are minor differences between them when = 0.02. More detailed stability behavior of VC models on the C L characteristic is plotted in Fig. 3d. In this example, an irregular oscillation can be observed with the OVC2 model, whereas the cases with the two vortex feature-based VC models provide periodic solutions. This discrepancy of the convergence history expresses a potential advantage of the vortex feature-based VC models in maintaining the robustness of the original ROSITA solver. From above analysis, the optimal value of the confinement parameter could be determined: o = 0.002 for the OVC2 model, o = 0.02 for both the FVC2-Q and FVC2-L2 models.

Influence on aerodynamic loads prediction
Different VC models have been tested to check the consistency of the aerodynamic loads ( C L , C D ), when the optimal parameters were applied. Due to the slightly oscillation behavior of the solutions, an average procedure was activated over 250 pseudo time steps at the end of the calculations. In Fig. 4, the computed sectional distribution of the aerodynamic loads ( C L , C D ) along the wing span are compared against the measured wind-tunnel data. It appears that, for the optimal value, although a good prediction of the sectional C L distribution is obtained by the OVC2 model case, the C D distribution is overestimated due to the application of confinement inside the boundary layer. Nevertheless, two vortex featurebased VC models provide more acceptable and robust results with minor discrepancies between them and with the simulation without VC. These improvements may be explained as the over-confinement in the boundary layer region is avoided by introducing the vortex feature-detection methods.

Effect on vorticity preservation
The contours of vorticity are plotted at five stream-wise stations (x/c = 0.05, 0.25, 0.45, 0.65, 0.85), as illustrated in Fig. 5. Figure 6 presents the formation procedure of wing-tip vortex structures at the onset phase calculated by different VC models with the optimal value and non-VC model with = 0.0. Two strong vortex systems, one from the suction side of the wing-tip and the other associated with the pressure side of the wing-tip, are shown. The improvement of vorticity for both two vortex systems could be seen from the solutions of the OVC2, FVC2-Q, and FVC2-L2 models. Meanwhile, compared with the OVC2 model case, the vortex feature-based VC models present better behavior of vorticity preservation. Figure 7 shows iso-surface of Q criterion at Q = 0.8 with the optiaml used. It is clear that the streamwise vortex shedding from the wing-tip is preserved much better by the vortex feature-based VC models. In addition, the FVC2-L2 model provides an almost identical distance of the downstream tip vortices with the FVC2-Q model. As for the OVC2 model case, although the downstream tip vortices are preserved better than the non-VC model, the vortices diffused much faster than the vortex feature-based VC models.
The normalized z-velocity through the wingtip vortex core is one of the indices to evaluate the performance of the VC models. The velocity profile is extracted by a line passing through the vortex core parallel to the span-wise direction, as indicated by Fig. 8. Figure 9 further quantify the improvement of the VC models on vorticity preservation by comparing the swirl velocity profile recorded at two-and four-chord lengths downstream of the trailing edge with measurements. As state previously, the OVC2 model presents the improved velocity profile but a fast diffusion downstream. In contrast, the vortex feature-based VC models are able to offer more desirable values. In addition, the improvement in the predicted maximum z-velocity for the vortex feature-based VC models are apparent with respect to the non-VC model. For example, at two chords downstream, the solutions of the FVC2-Q and FVC2-L2 models have increased the maximum z-velocity value by 43.9% and 44.3% , respectively, compared with the non-VC model. At location of four chords downstream, the increase becomes 32.4% and 33.1% , respectively. It can be said that the vortex feature-based VC models have somehow more robust values than the OVC2 model. Furthermore, these results demonstrate that

Caradonna-Tung rotor in hover
In this section, the flow around the Caradonna-Tung rotor blades, in hover, is used to demonstrate the performance of the implemented VC models on a threedimensional rotor flow. Due to the public availability of the experimental datasets, this rotor represents a benchmark to validate the VC models for hovering helicopter rotors.
The 2-bladed Caradonna-Tung model rotor has untwisted planform of 1.143 m radius. The blades are comprised of symmetric NACA0012 airfoil of 0.1905 m chord length. Experiments were carried out in the Army Aeromechanics Laboratory's hover test facility, where a large chamber with special ducting was designated to eliminate room recirculation [35]. The hover conditions considered here employ a blade collective pitch angle of 8 • , blade-tip Mach number of 0.877, and precone angle of 0.5 • .
A moving embedded grid system is designed for the present study, composed of background grid, and two identical body-fitted blade meshes. The grid system is presented in Fig. 10. The background mesh consists of two different cylindrical grids. A coarse grid (far-field grid) is generated to represent the flow domain far from the rotor, where outer boundaries are located 4R (above), 8R (below), and 8R (radial) away from the blade hub. A finer grid (near-field grid) is created to model the flow region close to the blades. The body-fitted blade grid is modelled with a C-H topology. The wall distance of the first layer of blade surfaces is set to 1 × 10 −5 c so that the y + value is less than 1. A non-slip boundary condition is applied on the blade surface.
To account for the grid sensitivity effect for the solutions around the blade, three structured blade grid, increasing the mesh size from 0.8 to 2.4 million cells, are generated to perform the grid sensitivity analysis. Furthermore, three near-field meshes are considered for the assessment of VC models on vorticity preservation and computational efficiency. The summary of the far-and near-field grids is shown in Table 3. The blade grid details are reported in Table 4 3.

Grid sensitivity study
A grid sensitivity study is carried out first to obtain the grid convergence solutions when no VC models are enabled. Computations are performed by employing the NG1 near-field mesh and three different blade meshes with coarse, medium and fine grid densities. Fig. 11 shows the comparison of the calculated C P distributions with experimental data [34] at three radial positions, r/R = 0.5, 0.89 and 0.96. At inboard station (r/R = 0.5), the sectional C p distributions are almost identical for three grid density cases. At two outboard positions (r/R = 0.89 and 0.96 ), the grid resolution effect becomes more notable on the prediction of transonic shock wave. A more smeared shock is observed for the coarsest grid case. However, the results of the medium and fine grid cases are in good agreement with the measurements with negligible difference between them. Therefore, it is reasonable to believe that sufficiently converged results are obtained with the medium blade grid, which can be used for further study.

Effect on computational stability
In terms of the computational stability of the implemented VC models, Fig. 12 shows the flow solution residuals of the Caradonna-Tung rotor blade ( NG1 near-field mesh for all the test approaches) under seven confinement parameter values, where = 0.0 stands for the solution without any VC models enabled. All flow solutions were simulated by solving the RANS equations, coupled with the Spalart-Allmaras turbulence model. The governing equations were integrated with the implicit dual-time stepping method of ROSITA, using a pseudo-time CFL number equal to 3.0. Typically, 4000 iterations are required to reduce the residuals by four levels of flow solutions in most cases. However, it is observed that, for the OVC2 solution with = 0.02 , some robustness problems occur after about 500 pseudo-time steps; therefore values greater than = 0.02 were not tested for the OVC2 model. The result for the FVC2-Q model with = 0.05 starts diverging after 6000 pseudo-time steps. The FVC2-L2 model is able to maintain the stability of the entire calculation process until = 0.07. This supports the idea that with the introduction of vortex feature-detection method, the robustness of the standard VC2 model is maintained with higher confinement parameters, especially for the FVC2-L2 model. Therefore, the optimal confinement parameter value can be roughly obtained: o = 0.01 for OVC2 model, o = 0.04 for FVC2-Q model, and o = 0.06 for FVC2-L2 model.

Influence on aerodynamic loads prediction
In this part, the medium blade mesh and NG1 nearfield mesh are employed for all tested approaches. Furthermore, the parameter is set at its optimal value for each VC model. The sectional lift coefficient C L are employed to investigate the effect of different VC models on aerodynamic loads prediction. An averaging process is activated over 500 pseudo-time steps at the end of each simulation. Figure 13 shows the variation of sectional C L along the rotor radial direction. It is observed that the FVC2-L2 model case provides notable improvements, although overpredictions are obtained for all calculations. As state previously, this improved result could be attributed to the use of the vortex feature detection methods, which avoid the over-confinement in the boundary layer.

Effect on vorticity preservation
As shown in Fig. 14, the contour plots of vorticity magnitude are extracted from five chord-wise sections (y/c = 0.5, 0.25, 0.45, 0.65, 0.85). The development of the tip vortex system over the Caradonna-Tung rotor blade is presented in Fig. 15, where the results with three different VC models are compared against the non-VC model solution. From the section x = 0.25c , it can be seen that the application of VC models has improved the coherent vortex structures. However, the improvements in the vortex formation region are not so evident.
Visualization of the Caradonna-Tung rotor wake using the iso-surface of Q-criteria (Q = 0.15) is shown in Fig. 16, in which the blade-tip vortex sections at five selected wake ages are extracted to better assess the downstream vorticity preservation capability of employed VC models. It should be reminded that, the VC models are employed with the optimal values; the non-VC model is used for comparison on the same near-field mesh. It is observed that the helical vortex filaments that are shed from the blade-tip are better preserved with the use of the implemented VC models when compared with the solution of non-VC model case. Furthermore, the capacities of the vortex feature-based VC models to preserve downstream vorticity are well demonstrated through the comparison with the solution of the OVC2 model. Quantitatively, Fig. 17 illustrated the vorticity magnitude of blade-tip vortex core as function of the wake age in degrees for the solutions with and without VC models on NG1 mesh. Although a dissipation of vorticity at the vortex core is observed along with the wake ages, the simulations performed with the VC models result in an improved vorticity preservation if compared with the solution without VC model. In detail, at wake ages of ∕6 , an improvement of vorticity by 16% , 20.8% and 25.1% , with respect to the non-VC case, is reported for the OVC2 model, FVC2-Q model and FVC2-L2 model, respectively. As the wake age increases, a more notable discrepancy of vorticity is witnessed for the solutions of three VC models. After 5 /3, the improvement of the core vorticity shown by the FVC2-L2 model is about twice as large as found for the FVC2-Q model if compared with the OVC2 model. These studies confirm that the introduction of vortex feature-detection methods can help to effectively enhance the capability of vorticity preservation of the OVC2 model. Moreover, the 2 -based VC2 model shows lower dissipation than the Q-based VC2 model.
The wake structures of the non-VC model case performed over two finer near-field meshes (NG2 and NG3) are also shown in Fig. 18. It is observed that, with the use of finer near-field meshes, the blade-tip vortex of the case without VC model can be traced up to about 11 ∕6 and 2 , respectively. However, the coherent structures are better preserved up to further distance by introducing the FVC2-Q and FVC2-L2 model on NG1 mesh, as previously illustrated in Fig. 16c and d.
This conclusion is further confirmed looking at the variations of the blade-tip vortex strength in Fig. 19. The vortex strength captured by the two finer nearfield meshes without VC model has a faster decay rate than those of the FVC2-Q and FVC2-L2 solutions. At the onset phase ( = ∕6 ), the core vorticity of the non-VC model with NG2 and NG3 meshes is greater than the results of the vortex feature-based VC  models with NG1 mesh. At the wake age of = ∕3 , the core vorticity value of the NG2 case turns to be lower than those of the NG1 cases. When it is close to the wake age of = 2 ∕3 , the calculated vorticity at the vortex core center for the NG3 case starts to become smaller than the vorticity value of the NG1 cases. In other words, it indicates that the vortex feature-based VC models express an excellent vorticity preservation capability. Table 5 reports the computational time cost for the Caradonna-Tung rotor blade simulations with the FVC2-Q and FVC2-L2 models on the NG1 mesh as well as the non-VC model on the NG1, NG2, and NG3 meshes. Solutions were computed on 64 cores of the high-performance cluster Galileo100 of CINECA, comprised of Intel CascadeLake 8260 2.4 GHz processors, nodes interconnected by a Mellanox Infiniband HDR100 high-performance network. As can be observed, cases with the vortex feature-based VC models on the coarse near-field mesh present a much higher computational efficiency than the non-VC model case performed on the finer near-field meshes.

Conclusions
In the present work, the implementation of the vortex feature-based second vorticity confinement (FVC2-Q and FVC2-L2) models in the ROSITA CFD solver has been presented to obtain the high resolution of vortical structure for three-dimensional vortex-dominated flows. To assess the performance of the vortex  non-VC model solutions in terms of computational stability, aerodynamic prediction, vortex resolution, and computational time consumption. It was found that the vortex feature-based VC models express better performance on above aspects. Furthermore, the present methods allowed higher confinement parameters on a coarse grid with a relatively high computational efficiency to achieve better solutions than those obtained without VC on a finer grid. In particular, the 2 -based VC2 model showed a higher resolution of the vortical structure, more robust computational procedure , and more accurate aerodynamic loads if compared with the Q-based VC2 model.
On the basis of the conclusions mentioned above, the future work will apply the vortex feature-based VC models for enhancing the prediction of the helicopter rotor blade-vortex interaction phenomenon.   Vorticity of the blade-tip vortex core, obtained with vortex feature-based VC models over NG1 mesh as well as the non-VC model over NG2 and NG3 meshes