Application of a Combinatorial Vortex Detection Algorithm on 2 Component 2 Dimensional Particle Image Velocimetry Data to Characterize the Wake of an Oscillating Wing

: To investigate the vortical wake pattern generated by water flow past an oscillating symmetric airfoil, using experimental velocity fields from particle image velocimetry (PIV), a novel combinatorial vortex detection (CVD) algorithm is developed. The primary goal is to identify and characterize vortices within the wake. Experimental flows introduce complexities not present in numerical simulations, posing challenges for vortex detection. The proposed CVD approach offers a more robust alternative, excelling in both vortex detection and quantification of essential parameters, unlike widely-used methods such as Q -criterion, λ 2 -criterion


Introduction
Vortices constitute fundamental flow features with widespread relevance in various fluid dynamic applications.Their formation in the wake of airfoils and streamlined bodies holds a paramount significance in both research and practical applications.In turbulent flow regimes, temporally evolving and spatially coherent structures, commonly identified as vortices [1], exert substantial influence.For instance, studying the airflow around a pitching airfoil not only aids in comprehending undesirable phenomena such as wing flutter [2] but also finds applications in understanding biologically inspired aquatic propulsion [3,4].Researchers have extensively compared the simulation of pitching airfoils using full Navier-Stokes and vortex models [5][6][7].In this context, precise detection and characterization of coherent wake structures emerge as central research objectives [2,8,9].Notably, at a given Reynolds number, the generation of various wake schemes hinges on the positioning and organization of vortices in the wake [3].
The classification of wakes characterized by the presence of vortices can be effectively depicted in a phase diagram that illustrates significant wake transitions as functions of two independent nondimensional variables [3,10,11].These parameters are the Strouhal number (St = fC/U), which is based on the flapping frequency (f ), the chord length (C), and the free stream velocity (U), and the dimensionless oscillation amplitude derived from the chord length.To construct accurate phase diagrams, it is crucial to identify vortices in the wake with high precision and reliability.Furthermore, beyond the fundamental task of vortex detection, it is important to compute several key parameters that are essential for describing the flow, especially in the context of aerodynamics.These parameters encompass crucial Fluids 2024, 9, 53 2 of 30 metrics such as the coordinates of vortex cores, their drift velocities, peak vorticity values, overall circulation, the radius defining vortex boundaries, and the highest circumferential velocity [12].
The existing methods for detecting vortices, such as the widely used Q-criterion [13], λ 2 -criterion [1], ∆-criterion [14][15][16], and Γ1-criterion [17], have inherent limitations [18][19][20].A primary challenge arises from their reliance on threshold values to discern vortices.Selecting an appropriate threshold is often subjective, and there may not be a universally optimal value across flow types.This is particularly pronounced in turbulent flow studies, marked by significant, multiscale, intermittent variations in both vorticity and strain rates.More significantly, the threshold sensitivity of these methods can lead to unreliable and inaccurate vortex identification [21,22].Excessively low thresholds cause frequent false positives, incorrectly labeling non-vortex regions as vortices.Conversely, high thresholds result in many missed detections, failing to identify real vortices in the flow.This tradeoff hinders finding a robust threshold that generalizes across the wide variability in real-world flow dynamics [23,24].
One core limitation identified in the Γ1-criterion, developed by Laurent Graftieaux et al. (2001) [17], is that the Γ1 function used to locate vortex centers lacks Galilean invariance, rendering it susceptible to the constant advection of the flow field.This inherent lack of robustness becomes particularly evident in scenarios where the overall flow undergoes consistent shifts.Additionally, Graftieaux's theoretical relations linking the Γ1 and Γ2 functions to vortex properties assume axisymmetry in the swirling flow structure.While this assumption simplifies the theoretical framework, real unsteady flows often deviate from the idealized axisymmetric case, introducing complexities that challenge the method's applicability to diverse flow conditions.Moreover, directly applying the Γ1-criterion to flows with multiple, interacting vortices poses challenges.Ambiguity can arise in attributing vortex indices, especially in scenarios where vortices dynamically interact, merge, or evolve over time.The method's reliance on theoretical relations may lead to difficulties in precisely characterizing complex vortex structures within the flow.Consequently, the method's effectiveness diminishes in situations where accurate identification and differentiation of interacting vortices are crucial [17].
To provide a robust approach for vortex detection, researchers often use combinations of criteria or data-driven techniques such as machine learning to adaptively determine optimal thresholds for a given flow [25][26][27].However, subjectivity and extensive parameter tuning may still be required.Furthermore, while the Q-criterion, λ 2 -criterion, and ∆-criterion are effective in binary classification, categorizing regions as vortical or non-vortical, they do not inherently provide direct quantification of vortex strength or intensity [28].Evaluating vortex strength is often crucial for in-depth analysis of vortex dynamics, comparative studies between vortices, and various related objectives [29,30].Researchers typically need to perform supplementary analyses to calculate circulation, maximum vorticity, core radius, or other relevant parameters to obtain strength information [29,30].
Recently, advanced vortex identification techniques such as the Ω-method [24] and 'Rortex' method [31,32] have shown promise in addressing some of the limitations of traditional criteria.The Ω-method defines vortices as regions where vorticity overtakes deformation, providing a clear physical interpretation.It does not require arbitrary threshold selection and can simultaneously capture strong and weak vortex structures [24].The 'Rortex' method decomposes vorticity into rotational and non-rotational parts, excluding shear contamination.As a vector quantity, 'Rortex' identifies both vortex axis orientation and strength [31,32].These methods offer greater objectivity and physical insight compared to conventional techniques.However, the Ω-method and 'Rortex' method have thus far been predominantly applied to direct numerical simulations (DNSs) data [24,28,[31][32][33][34].Their effectiveness in processing experimental measurements with inherent noise remains to be rigorously validated.
Substantial effort has been dedicated to the development and validation of advanced techniques using DNS data [1,35,36].DNS provides high-resolution flow field information for testing novel vortex identification strategies, especially in three-dimensional turbulent flows such as boundary layers [37,38].However, DNS is limited to relatively simplified flow conditions and small domains due to computational constraints.Therefore, while beneficial for initial development, additional rigorous validation of new methods with experimental data remains essential before broader utilization across complex flows.DNS aids technique refinement, but ultimate validation still requires experiments.
Detecting vortices in particle image velocimetry (PIV) data poses unique challenges due to inherent measurement noise and experimental uncertainties [39].Typically, mitigating measurement noise involves prior filtering before applying a detection algorithm [12].However, selecting appropriate thresholds becomes crucial to avoid overlooking weak or small vortices [39].Additionally, the flow field near the vortex core exhibits substantial gradients, potentially causing difficulties in particle seeding and correlating PIV data [12].
In contrast to DNS data that can utilize local mesh refinement for better spatial detail in areas like vortex centers, PIV analysis methods must seek alternative strategies to effectively manage the pronounced gradients encountered around the vortex core [40].
Due to camera resolution limits, PIV analysis often requires integrating images from several cameras, necessitating image de-warping and stitching, which might lead to discontinuities in the flow-field representation, potentially leading to the erroneous identification of elevated local vorticity as a vortex [1].
Given the increasing availability of PIV data in experimental fluid mechanics [12], the need for developing dedicated and effective algorithms for vortex detection and characterization in PIV datasets has emerged as a pivotal area of research.The precise definition of a vortex remains elusive [1], making their detection in practical flows a complex task.The most common and intuitive conception of a vortex revolves around the perception of swirling fluid motion centered around a focal point [35].However, translating this intuitive notion, rooted in human perception, into a rigorous numerical characterization detectable by algorithmic methods poses considerable challenges.Lugt (1979) [41] proposed that "A vortex is the rotating motion of a multitude of material particles around a common center", aligning with the fundamental notion of swirling motion.While this concept is intuitively appealing, it lacks specific criteria for describing vortical structures and does not easily translate into a workable vortex identification algorithm.Building upon Lugt's concept, Robinson (1991) [36] offered a more geometrically precise definition: "A vortex exists when instantaneous streamlines, projected onto a plane perpendicular to the vortex core, exhibit a roughly circular or spiral pattern when observed from a reference frame moving with the center of the vortex core".While this definition offers clarity in geometric characterization, it necessitates a priori information about the vortex core's location, posing challenges for systematic vortex detection.In contrast, Jeong and Hussain (1995) [1] put forth the notion that "A vortex core must possess net vorticity, hence net circulation".This concept offers a direct method to pinpoint potential vortex regions without requiring prior knowledge of the vortex core's location or drift velocity.By computing scalar vorticity fields from velocity vectors, researchers can generate and label regions of interest.Combining Jeong and Hussain's definition [1] with Robinson's geometric criteria [36] forms a solid groundwork for creating an effective vortex detection algorithm.
Vortex identification techniques are often designed for specific data types and unique flow fields, aiming to identify well-defined vortex structures based on precise definitions [35].These algorithms exhibit a high degree of sensitivity to various parameters, including vortex spacing, particularly when multiple vortices coexist within a given flow field, vortex size, angular velocity, and the presence of shear flow [12].Sadarjoen and Post (2000) [42] proposed a classification of vortex detection algorithms found in the literature, distinguishing them into two primary categories.The first category encompasses traditional vortex detection algorithms, which rely on the evaluation of physical flow-field properties at specific points.Commonly assessed properties include velocity, vorticity, pressure, and the velocity gradient tensor.However, one notable limitation of these algorithms, as emphasized by Sadarjoen and Post (2000) [42], is their reduced sensitivity to large, slowly rotating vortices.For instance, the maximum vorticity method defines a vortex core based on local vorticity maxima [35].Inherently, this method disregards regions where vorticity falls below a predefined threshold, potentially leading to the failure of detecting weak, slowly rotating vortices.
The second category comprises geometric methods, representing a more recent set of vortex detection algorithms.These methods prioritize the examination of flow patterns and trajectories, such as streamlines and pathlines, over the analysis of scalar attributes measured at distinct grid locations [35].Although geometric methods demand more computational resources, they offer a significant advantage over traditional techniques by effectively differentiating between actual vortices and false positives, addressing a common limitation in earlier approaches [35,42].
The limitations of the most widely adopted vortex identification criteria, including their sensitivity to threshold values and their inability to quantify vortex strength directly, underscore the pressing need for complementary techniques, modifications to existing criteria, combinations of criteria, and the exploration of emerging methods to enable accurate and comprehensive vortex detection across a wide range of fluid mechanic applications.
This work introduces a vortex detection algorithm that combines selected features from three individual algorithms to interrogate PIV data.The primary aim of this algorithm is to efficiently process data, detecting and locating vortices while calculating characteristic vortex parameters.Given the potential magnitude of the datasets associated with the PIV results, the algorithm also incorporates automated procedures for identifying false positive vortices without the need for human intervention.The algorithm's performance is evaluated by analyzing vortex behavior in the extensively studied wake behind a symmetrically pitching airfoil [2,10].

Background on Existing Vortex Detection Algorithms
The development of an effective vortex detection algorithm relies on the fundamental principles and methodologies established in the existing literature.This background reviews three pivotal methods that serve as the basis for the proposed combinatorial vortex detection (CVD) framework tailored for the analysis of PIV data.These methods encompass the maximum vorticity (MV) method [43], the cross-sectional lines (CSL) method [12], and the winding angle (WA) method [44].The following sections delve into the core principles of each method and discuss key aspects leveraged in developing the CVD algorithm.

Maximum Vorticity (MV) Method
The MV method, proposed by Strawn et al. (1999) [43], assumes that a vortex core exists at locations characterized by a local maximum in vorticity magnitude.This approach proves effective in detecting vortices in situations where their cores may overlap [35].However, a notable limitation of the MV method is its propensity to not only identify vortices but also regions exhibiting shearing activity [1].Consequently, devising a vortex detection algorithm solely based on vorticity fields becomes more challenging, particularly in non-free shear flow [1,35].An alternative utility of vorticity fields lies in outlining regions of interest (ROIs) where potential vortices may be located, enabling the use of other algorithms to evaluate these ROIs individually.
To enhance the performance of the MV method within the framework of the CVD method, a Gaussian low-pass spatial filter is applied to the velocity data to mitigate measurement noise and experimental uncertainties [39].This filtering process produces a smoother field, effectively eliminating small-scale fluctuations that might otherwise lead to localized spikes in the vorticity field.Scalar vorticity, ω, calculated numerically at each grid point, is defined within the two-dimensional velocity field, In the pursuit of detecting weak vortices characterized by relatively low vorticity, it becomes necessary to apply a low vorticity threshold.However, this approach carries the inherent risk of grouping stronger, nearby vortices into a single ROI, effectively evaluating Fluids 2024, 9, 53 5 of 30 them as a singular vortex structure.To mitigate this, multiple threshold levels are used, and the indexed image is further analyzed using logical operations and image morphology (IM) techniques to distinguish between closely situated vortices [45].
The application of multilevel thresholds serves the dual purpose of allowing the identification of weak vortices while ensuring that even when strong vortices are close together or their cores overlap, they can still be distinctly identified and analyzed as separate entities.Figure 1 provides a visual representation of this concept, depicting two vorticity field test cases where a three-level thresholding approach is applied.The resulting indexed images, accompanied by the corresponding desired ROI maps, are presented for these sample vorticity fields.The left-hand side of Figure 1 illustrates the indexed images, where index levels are denoted as i = 1, 2, and 3. On the right, the desired ROI map is displayed, signifying the regions to be extracted from the corresponding indexed image.
and the indexed image is further analyzed using logical operations and image morphology (IM) techniques to distinguish between closely situated vortices [45].
The application of multilevel thresholds serves the dual purpose of allowing the identification of weak vortices while ensuring that even when strong vortices are close together or their cores overlap, they can still be distinctly identified and analyzed as separate entities.Figure 1 provides a visual representation of this concept, depicting two vorticity field test cases where a three-level thresholding approach is applied.The resulting indexed images, accompanied by the corresponding desired ROI maps, are presented for these sample vorticity fields.The left-hand side of Figure 1 illustrates the indexed images, where index levels are denoted as i = 1, 2, and 3. On the right, the desired ROI map is displayed, signifying the regions to be extracted from the corresponding indexed image.
Without the utilization of multilevel thresholds, a single threshold level would be chosen, leading to one of two scenarios.When solely threshold level i = 1 is applied, as seen in both examples, the broader vorticity region, clearly exhibiting two distinct peaks of vorticity, is incorrectly grouped into one single ROI, consequently treated as a single potential vortex.However, this area features two distinct peaks and should rightfully be considered as two individual ROIs.Conversely, when exclusively threshold level i = 3 is employed, the weaker vortices failing to meet the vorticity ωpeak level i = 3 threshold are excluded from consideration as ROIs and remain unidentified as potential vortices.Without the utilization of multilevel thresholds, a single threshold level would be chosen, leading to one of two scenarios.When solely threshold level i = 1 is applied, as seen in both examples, the broader vorticity region, clearly exhibiting two distinct peaks of vorticity, is incorrectly grouped into one single ROI, consequently treated as a single potential vortex.However, this area features two distinct peaks and should rightfully be considered as two individual ROIs.Conversely, when exclusively threshold level i = 3 is employed, the weaker vortices failing to meet the vorticity ω peak level i = 3 threshold are excluded from consideration as ROIs and remain unidentified as potential vortices.

Removed by erosion
To address the notable limitations associated with the use of a single-level threshold for identifying ROIs, the algorithm employs a three-level threshold to generate an ROI map.In this approach, the vorticity field is segmented, resulting in a classified image using threshold values specified by the threshold intensity vector (TIV): Subsequently, the algorithm executes morphological opening for each of the thresholds specified in Equation (1).Morphological opening combines two fundamental image morphology techniques, erosion and dilation, performed in a specific sequence.Initially, the image undergoes erosion using a diamond-shaped structuring element of size S e .This operation effectively eliminates clusters of pixels smaller than S e , while reducing the size of clusters larger than S e .
Following erosion, a dilation step is performed with the same structuring element S e to bring the pixel clusters, which remained post-erosion, back to their original size.This combination of erosion and dilation, performed in a well-defined sequence, plays a pivotal role in refining the ROI map and contributes to enhancing the algorithm's efficacy in vortex detection.
Connected pixel groups in the vorticity intensity indexed image are categorized by intensity levels (I 1 , I 2 , I 3 ).Groups associated with index I 1 , denoting the lowest vorticity intensity level, undergo an initial morphological opening procedure, employing an IM diamond-shaped structuring element of size S e1 .Subsequently, these groups are labeled and individually investigated.Structures corresponding to index I 2 also undergo morphological opening, but this time utilizing an IM diamond-shaped structuring element of size S e2 .They are then sub-labeled within each of the previously labeled groups from index I 1 .
The algorithm follows a conditional path: If none of the groups from the second threshold index I 2 are found within a group from the first level I 1 , the algorithm concludes, and the group of pixels from index level I 1 is designated as an ROI.Conversely, if at least one I 2 group is found within a group from index I 1 , the algorithm progresses further.Pixels corresponding to index I 3 undergo morphological opening using an IM diamondshaped structuring element of size S e3 and are further categorized within each of the index I 2 groups.
At this stage, the algorithm enumerates the group counts at each threshold level.If an I 2 cluster encloses several I 3 groups, each individual index I 3 group becomes a distinct ROI.However, if only one index I 3 group is counted, the algorithm reverts to a lower level, considering the number of index I 2 clusters.Once again, if multiple index I 2 groups are detected, each individual index I 2 group is recognized as an ROI.However, if only one index I 2 label is counted, the group from index I 1 becomes the designated ROI.This process ensures that the ROI encompasses the largest possible area while partitioning into multiple ROI's when distinct, higher vorticity peaks are detected.

Cross-Sectional Lines (CSL) Method
The CSL method is designed to identify potential vortex cores within a designated ROI by evaluating the velocity component perpendicular to parallel straight cross-sectional lines intersecting the ROI at arbitrary angles [12].In this context, the CSL method is integrated into the CVD framework for its proven ability to accurately pinpoint vortex cores within ROIs labeled using the MV algorithm.
It is essential to clarify that the CSL method does not serve as a vortex detection mechanism itself.Instead, its primary function revolves around locating potential vortex cores and defining boundary radii within the ROIs.The method provides two-dimensional coordinates for a core, whether or not an actual vortex exists within the current ROI.
Each labeled ROI undergoes a CSL algorithm akin to the one described by Vollmers (2001) [12].A schematic outlining the fundamental CSL procedure is presented in Figure 2, featuring a sample vortex for illustrative purposes.The schematic employs concentric ellipses to represent vorticity contours, with the vortex core identified by an '×'.
Each labeled ROI undergoes a CSL algorithm akin to the one described by Vollmers (2001) [12].A schematic outlining the fundamental CSL procedure is presented in Figure 2, featuring a sample vortex for illustrative purposes.The schematic employs concentric ellipses to represent vorticity contours, with the vortex core identified by an '×'.In the context of a vector velocity field with dimensions Α × Β, the y-component of velocity, represented as  (, ), is assessed for columns i = 1, 2, 3, …, Α and rows j = 1, 2, 3, …, Β along each row within the defined grid.During the discussion of the CSL method, this velocity component,  (, ), is referred to as the "perpendicular velocity," denoted as  (, ).
The determination of maximum and minimum perpendicular velocities for each row, designated as j, within a specified ROI is a pivotal aspect of the analysis.These extreme values are calculated as follows: , ( ) = min ( ( , ),  ( 1, ), … ,  ( , )) The indices As and Ae are used to denote the first and last velocity vectors within a ROI on a given row, respectively.Within the ROI, the row ψ that exhibits the difference  , ( ) −  , ( ) is identified as the critical cross-section line and is represented as CSLψ.This 'critical line' corresponds to the y-coordinate of a potential vortex core, should one exist within the ROI under evaluation.
The perpendicular velocity along this critical row is denoted as  (, ).The x-coordinate of the vortex core is determined by the location along CSLψ where the velocity satisfies the condition and the coordinates of the vortex core are, therefore, represented as (ς, ).
For the sake of simplicity, Figure 2 illustrates a hypothetical vortex with only three rows/lines of the  velocity field.In this example, the critical perpendicular line is denoted as CSLψ = CSL2.Among the three lines, CSL2 exhibits the most significant difference between its maximum and minimum perpendicular velocities,  , ( ) and  , ( ) , respectively.Consequently, the y-coordinate for the hypothetical vortex core is determined as ψ = 2.
Using Equation (4) along CSL2, the velocity  (ς, 2) is evaluated.The algorithm searches for the point along CSL2 where  (ς, 2) matches, thus finding the value of ς.In the context of a vector velocity field with dimensions A × B, the y-component of velocity, represented as v y (i, j), is assessed for columns i = 1, 2, 3, . . . A and rows j = 1, 2, 3, . .., B along each row within the defined grid.During the discussion of the CSL method, this velocity component, v y (i, j), is referred to as the "perpendicular velocity," denoted as v p (i, j).
The determination of maximum and minimum perpendicular velocities for each row, designated as j, within a specified ROI is a pivotal aspect of the analysis.These extreme values are calculated as follows: The indices A s and A e are used to denote the first and last velocity vectors within a ROI on a given row, respectively.Within the ROI, the row ψ that exhibits the difference v p,max (j) − v p,min (j) is identified as the critical cross-section line and is represented as CSL ψ .This 'critical line' corresponds to the y-coordinate of a potential vortex core, should one exist within the ROI under evaluation.
The perpendicular velocity along this critical row is denoted as v p (i, ψ).The x- coordinate of the vortex core is determined by the location along CSL ψ where the velocity satisfies the condition and the coordinates of the vortex core are, therefore, represented as (ς, ψ).
For the sake of simplicity, Figure 2 illustrates a hypothetical vortex with only three rows/lines of the v y velocity field.In this example, the critical perpendicular line is denoted as CSL ψ = CSL 2 .Among the three lines, CSL 2 exhibits the most significant difference between its maximum and minimum perpendicular velocities, v p,max (2) and v p,min (2) , respectively.Consequently, the y-coordinate for the hypothetical vortex core is determined as ψ = 2.
Using Equation (4) along CSL 2 , the velocity v p (ς, 2) is evaluated.The algorithm searches for the point along CSL 2 where v p (ς, 2) matches, thus finding the value of ς.These two computed coordinates effectively pinpoint the vortex core of a potential vortex within the ROI.However, it is important to note that this process alone does not confirm the presence of an actual vortex within the ROI.
The CSL method also offers an estimation of vortex size by fitting a circle with a radius r v centered at the vortex core, which is determined as Here, ς v,max (ψ) represents the column position on row ψ where the velocity vector v p,max (ψ) is located, and ς v,min (ψ) refers to the column position on row ψ where the velocity vector v p,min (ψ) is found.
The CSL method proves to be particularly effective for identifying vortex cores and estimating vortex size and drift velocity in flow fields where available information regarding the vortical structures present is scarce.Additionally, Vollmers (2001) [12] proposed that the transverse and streamwise components of the vortex drift velocity, represented as ⇀ v dri f t = v dri f t,x , v dri f t,y , correspond to the x and y components of velocity evaluated at the core coordinate (ς, ψ).This algorithm operates on a predefined ROI, which is de- termined using the MV method.The algorithm comprises two main nested loops: one iterating over individual ROIs and the second for individual cross-sectional lines within each ROI.The input data include the 2C2D velocity field and the labeled ROI map.
A drawback of the CSL method lies in its inherent inability to reject false positives.It may identify a vortex core in an ROI that does not actually contain a vortex.To address this challenge, a hybrid detection algorithm is proposed that combines physical quantitybased methods, such as MV and CSL, with geometric methods integrated into the CVD framework, as geometric methods are known for their effectiveness at identifying and rejecting false positives, a capability often overlooked by other methods.

Winding Angle (WA) Method
The WA method, initially proposed by Portela in 1999 [44], represents a geometric approach to vortex detection.This method involves the assessment of discretized streamlines to determine their affiliation with a vortex.According to the criteria suggested by Sadarjoen and Post in 2000 [46], a streamline must meet two conditions to be classified as part of a vortex.
Firstly, for a streamline denoted as S k , the winding angle α w,k must satisfy the equation and ⇀ V 2 on a given streamline is expressed as where and and Here, P x and P y are the respective (x,y) location on a given streamline and ⇀ V n is a unit vector normal to the flow-field of interest.Figure 3 demonstrates how the individual angles α w,i are calculated on a section from streamline S k .The winding angle α w,k for a streamline  A crucial characteristic of a vortex-associated streamline is the formation of a closed semi-elliptical path [1,47,48].However, relying solely on the winding angle criterion is insufficient to meet this requirement [46].It is possible for a streamline to exhibit a substantial winding angle without tracing a closed path.Therefore, an additional condition is necessary.
The second requirement dictates that a streamline deemed part of a vortex must have its initial and terminal points in close proximity [42].The distance between the starting point ( , ,  , ) and the endpoint ( , ,  , ) of a streamline Sk is calculated using The initiation point of a streamline is predetermined, and each streamline consists of a fixed number of points with uniform spacing between them.Consequently, for a streamline tracing a closed path, its starting and ending points can be anywhere along the path, as the streamline's length is pre-defined.To address this challenge in establishing a Dse threshold, consider a scenario where streamlines follow a closed circular path.In such cases, the maximum separation between the starting and ending points of any streamline following a closed circular trajectory would equal the circle's diameter, which outlines the path.This presents a significant challenge in setting a Dse threshold, as streamlines associated with larger vortices would require a larger Dse.However, an excessively large Dse threshold may incorrectly associate non-vortical streamlines with relatively small vortices.To mitigate this issue, streamlines are split whenever the cumulative angle sum reaches the highest multiple of 2π.This technique ensures that the beginning and end points of a streamline on a closed loop remain as close as possible.
While the WA threshold is defined as  , = 2π, selecting a suitable value for the threshold Dse is necessary.The Dse threshold depends on the length scale of the vortical structures present, but the previously defined ROIs can aid in determining a sensible Dse threshold.Within each ROI, circles with a radius  are used to identify potential vortices.This parameter serves as a reasonable length scale for coherent structures in the flow.Therefore, an appropriate choice for the Dse threshold would be a value greater than  2 ⁄ .
This criterion ensures that only streamlines with minimal endpoint separation are considered as part of a vortex.Streamlines that extend significantly beyond the typical size of flow structures are excluded to maintain the accuracy of vortex identification.The final step of the algorithm involves verifying whether individual streamlines adhere to the specified thresholds.Streamlines meeting all required criteria are marked and assigned unique identifiers.The sign of the winding angle  , for each streamline determines the direction of rotation associated with the vortex corresponding to that streamline.A crucial characteristic of a vortex-associated streamline is the formation of a closed semi-elliptical path [1,47,48].However, relying solely on the winding angle criterion is insufficient to meet this requirement [46].It is possible for a streamline to exhibit a substantial winding angle without tracing a closed path.Therefore, an additional condition is necessary.
The second requirement dictates that a streamline deemed part of a vortex must have its initial and terminal points in close proximity [42].The distance between the starting point (P x,1 , P y,1 ) and the endpoint (P x,N , P y,N ) of a streamline S k is calculated using The initiation point of a streamline is predetermined, and each streamline consists of a fixed number of points with uniform spacing between them.Consequently, for a streamline tracing a closed path, its starting and ending points can be anywhere along the path, as the streamline's length is pre-defined.To address this challenge in establishing a D se threshold, consider a scenario where streamlines follow a closed circular path.In such cases, the maximum separation between the starting and ending points of any streamline following a closed circular trajectory would equal the circle's diameter, which outlines the path.This presents a significant challenge in setting a D se threshold, as streamlines associated with larger vortices would require a larger D se .However, an excessively large D se threshold may incorrectly associate non-vortical streamlines with relatively small vortices.To mitigate this issue, streamlines are split whenever the cumulative angle sum reaches the highest multiple of 2π.This technique ensures that the beginning and end points of a streamline on a closed loop remain as close as possible.
While the WA threshold is defined as α w,k = n2π, selecting a suitable value for the threshold D se is necessary.The D se threshold depends on the length scale of the vortical structures present, but the previously defined ROIs can aid in determining a sensible D se threshold.Within each ROI, circles with a radius r v are used to identify potential vortices.This parameter serves as a reasonable length scale for coherent structures in the flow.Therefore, an appropriate choice for the D se threshold would be a value greater than r v /2.This criterion ensures that only streamlines with minimal endpoint separation are considered as part of a vortex.Streamlines that extend significantly beyond the typical size of flow structures are excluded to maintain the accuracy of vortex identification.
The final step of the algorithm involves verifying whether individual streamlines adhere to the specified thresholds.Streamlines meeting all required criteria are marked and assigned unique identifiers.The sign of the winding angle α w,k for each streamline determines the direction of rotation associated with the vortex corresponding to that streamline.
In situations where an ROI encompasses multiple vortices, it becomes essential to appropriately label all the closed streamlines associated with these vortices.This labeling process involves mapping each of the closed streamlines to a specific point and subsequently identifying clusters of closely situated points [42].The mapping of the streamlines S k belonging to a particular vortex is achieved by associating them with corresponding points MP k , which are determined by and To ascertain the distance between two such mapped points, MP k and MP k' , the following computation is applied: The labeling process begins with the first point, MP 1 , serving as the foundation for the initial cluster group.When considering point MP 2 , a crucial criterion is whether the distance between MP 1 and MP 2 falls within a specified tolerance.If it does, MP 2 is assigned to the existing group labeled as group 1.Conversely, if the distance exceeds the defined tolerance, MP 2 initiates the formation of a new cluster group.
Moving on to point MP 3 , the algorithm evaluates two distances, namely D MP 3,1 ′ and D MP 3,2 .If neither of these distances meets the tolerance criteria, MP 3 takes the role of establishing yet another new cluster group.However, if either D MP 3,1 ′ or D MP 3,2 falls within the tolerance range, MP 3 is assigned to the nearest existing group, ensuring optimal grouping based on proximity.
This meticulous labeling process is applied consistently to all other points, ensuring that each point is appropriately categorized within the designated cluster groups.Finally, after all the streamlines have been meticulously labeled, the algorithm proceeds to compute the arithmetic average for each cluster group.
The WA algorithm implemented in this study utilizes a three-tiered nested loop structure: the first loop investigates the ROIs, the second loop delves into the instantaneous streamlines within each ROI, and the third loop focuses on the individual grid points comprising each streamline.This intricate process operates on the input data, which consist of streamlines computed from local velocity fields with a reference velocity equal to the drift velocity of the corresponding ROI.
It is important to note that the WA method's accuracy is closely tied to the precise selection of the reference velocity for streamline computation.Deviations from the actual vortex core's drift velocity can lead to inaccuracies in estimating the vortex boundary.While the WA method may not excel in providing accurate estimates of vortex size and shape, it exhibits remarkable robustness in effectively eliminating false positives, as demonstrated by previous research [42].Additionally, it serves as a valuable means of validating the accuracy of other vortex detection techniques.
Typically, verifying the presence of a vortex in a flow field involves visual inspection, but this approach becomes impractical when dealing with large datasets.As such, the WA method assumes the role of a swirling flow verification tool once potential vortex core locations and drift velocities have been computed through other detection methods.

Comparison of MV to other Detection Methods
Vorticity magnitude alone cannot distinguish between shear generated and free shear vorticity.However, in the current study, MV serves only as an initial means of identifying ROIs.Detection methods such as Q-criterion, λ 2 -criterion, or ∆-criterion are far less likely to identify a shear region, as a vortex, than the MV method and ultimately produce less false positives but are more complex [12,49].The simplicity of the MV method makes it a good candidate for initially identifying ROIs in the flow and to effectively study the WA's ability to reject shear generated ROIs by deliberately identifying some of these shear regions as potential vortices.It is also important to consider that the vortices in the pitching airfoil flow-field have vastly different vorticity magnitudes.The Q, λ 2 , and ∆ methods have difficulty distinguishing between individual vortices with potentially overlapping cores and selecting the right threshold in order to distinguish individual vortices can be complex [49].For example, in incompressible flows, the Q-criterion identifies vortices as regions where the vorticity magnitude prevails over the strain-rate magnitude [13,18].It accomplishes this by finding connected regions in the flow field having a second invariant of the velocity gradient tensor that is less than a negative threshold value II E [13] shown as The second invariant of the velocity gradient tensor (II) is expressed as As a result, vortices with overlapping cores will blend together and appear as a single vortex without the proper threshold value −II E .
When the entire vortex identification procedure is carried out in a single method, and it is not guaranteed that the flow is free shear, then it is favorable to use a more elaborate method such as the Q-criterion, ∆-criterion, or λ 2 -criterion.
In summary, MV is favored as an ROI identifier for several reasons: • MV is not employed as a standalone detection method; instead, it is integrated with multilevel thresholds, CSL, and WA methods.• The flow field comprises vortices with a wide range of magnitudes and potentially overlapping cores.Therefore, it is recommended that the use of multilevel thresholds and adjusting WA thresholds is more intuitive, given that the threshold unit is expressed as a vorticity magnitude.

•
Unlike other methods such as Q-criterion, ∆-criterion, or λ 2 -criterion, MV provides information about the direction of rotation.This additional feature aids in distinguishing closely spaced vortices by symmetrically extending the multilevel threshold.

•
MV's simplicity is valuable in the development of a combinatorial method as it facilitates the isolation of each method's contribution.

•
For the purpose of evaluating WA's effectiveness in identifying and eliminating false positives, it is useful to include some shear-generated ROIs (false positives).• The computation of derivatives required by the Q-criterion, ∆-criterion, or λ 2 -criterion methods is computationally demanding and sensitive to noise, which is often present in experimental data.

New Combinatorial Vortex Detection (CVD) Method
The CVD method has been developed through the integration of three distinct detection algorithms, MV, CSL, and WA, alongside straightforward image morphology techniques.The primary objective in crafting this method was to ensure its consistent and reliable capability to detect and characterize multiple vortices within velocity vector maps generated by PIV.The CVD method's functionalities encompass labeling each vortex within the flow field, pinpointing vortex cores, determining drift velocities, calculating circulation, assessing peak vorticity, and estimating boundary radii for individual vortical structures.Above all, the method serves the crucial purpose of reducing the dimensions of the original dataset while preserving essential vortex parameters accurately.
The CVD algorithm is succinctly outlined in the flow chart presented in Figure 4. Within the flowchart, the three individual detection algorithms are represented by green boxes, and the input and outputs are indicated in grey boxes, encompassing the 2C2D velocity vector field and the characterized vortex field, respectively.To begin, a 2D scalar vorticity field is derived from the global velocity map.This vorticity field undergoes indexing through the multilevel threshold algorithm, leading to the generation and labeling of sensible ROIs via the MV technique.

Verification on a Simulated Flow Field
To determine the minimum required number of velocity vectors, obtained through PIV, for the CVD algorithm to accurately identify and describe a vortex, an investigation was conducted using a simulated flow field measuring 15 mm × 15 mm.This simulated field contained an axisymmetric Burgers vortex positioned at the center of the region.Various levels of discretization of the velocity field were employed, enabling an examination of vector density per vortex.
The Burgers vortex solution serves as a well-known model for illustrating key aspects of modern turbulence theory [50].This vortex model offers an exact solution to the cylindrical Navier-Stokes equations, depicting the flow on a cylindrical vortex core inducing a circulation denoted as  at large distances [51].For an axisymmetric Burgers vortex in a fluid with kinematic viscosity ν and radius  , the distribution of circumferential velocity  and vorticity ω are [52] and Subsequently, the CSL method is employed to evaluate each ROI individually, providing coordinates for the vortex core, boundary radius, and drift velocity.One of the CVD method's standout features is its reliance on Galilean invariant indicators for vortex core detection.The CSL method, which searches for critical points along perpendicular velocity profiles, ensures consistent and reliable identification of vortex cores, independent of uniform flow shifts.This robust approach significantly broadens the method's applicability in dynamic flow scenarios.However, it is essential to note that the CSL algorithm yields these vortex parameters regardless of the presence of an actual vortex, underscoring its inability to ascertain vortex existence within a given ROI.
The WA method serves the purpose of validating the presence of a vortex within an ROI.This is achieved by the WA method searching each individual ROI for closed streamlines and automatically confirming or denying the presence of one or multiple vortices.By identifying vortex-affiliated streamlines based on their inherent geometry, independently of any symmetry assumptions, this approach allows for the accurate detection of vortices in complex, asymmetric flows, even in scenarios involving significant vortex interactions.
It is noteworthy that, for consistency with Robinson's definition [36], the WA method necessitates the computation of instantaneous streamlines from a reference frame moving with the vortex core [12].Consequently, prior knowledge of the vortex, specifically its core coordinates and drift velocity, is required.This poses a unique challenge as, unlike CSL and MV, WA cannot be conducted within the global coordinate system.Instead, local velocity vector maps must be generated for each ROI, with the drift velocity, calculated from CSL, subtracted before computing the streamlines.Ultimately, the WA method either accepts an ROI and labels it as a vortex or rejects it.
The CVD algorithm has been implemented in MATLAB.The code is freely available on GitHub (link in Supplementary Material).This includes an example vector field generated from a commercial PIV software (DaVis 8.1, LaVision GmbH Göttingen Germany) and examples in the form of videos of the oscillating wing, determined vorticity field, and detected vortices.

Verification on a Simulated Flow Field
To determine the minimum required number of velocity vectors, obtained through PIV, for the CVD algorithm to accurately identify and describe a vortex, an investigation was conducted using a simulated flow field measuring 15 mm × 15 mm.This simulated field contained an axisymmetric Burgers vortex positioned at the center of the region.Various levels of discretization of the velocity field were employed, enabling an examination of vector density per vortex.
The Burgers vortex solution serves as a well-known model for illustrating key aspects of modern turbulence theory [50].This vortex model offers an exact solution to the cylindrical Navier-Stokes equations, depicting the flow on a cylindrical vortex core inducing a circulation denoted as Γ ∞ at large distances [51].For an axisymmetric Burgers vortex in a fluid with kinematic viscosity ν and radius r v , the distribution of circumferential velocity v θ and vorticity ω are [52] and where the parameter γ represents the axial strain (∂w/∂z) within a velocity field ⇀ u described by an irrotational pure strain component denoted as ⇀ u s = (αx, βy, γz) and a rotational component confined to the x-y plane denoted as ⇀ u w = u x , u y , 0 .Specifically, for the case of an axisymmetric Burgers vortex, γ > 0 and β = α = −γ/2.
The simulated flow field consists of (n,m) velocity vectors, denoted as ⇀ v (i, j) with indices i in the x-direction and j in the y-direction, expressed by In this simulation, an axisymmetric Burgers vortex is modeled, which travels in the negative y-direction at a freestream velocity of U ∞ = −2.25 mm/s.The vortex exhibits a peak vorticity of ω peak = −3.96s −1 , a boundary radius of r v =3.18 mm, and induces a circulation of Γ ∞ =100 mm 2 /s at far distances.
The effects of changing the grid resolution of the flow field and introducing Gaussian white noise to the velocity vectors are investigated.This simulation aims to establish the minimum grid resolution and the acceptable level of velocity field noise that ensures the accurate application of the CVD algorithm to a universal flow field.Figure 5 illustrates the simulated vector field, with the background color map representing velocity magnitude.Additionally, Figure 6 displays the simulated velocity field after subtracting the freestream velocity, accompanied by a background color map representing vorticity.
minimum grid resolution and the acceptable level of velocity field noise that ensures the accurate application of the CVD algorithm to a universal flow field.Figure 5 illustrates the simulated vector field, with the background color map representing velocity magnitude.Additionally, Figure 6 displays the simulated velocity field after subtracting the freestream velocity, accompanied by a background color map representing vorticity.minimum grid resolution and the acceptable level of velocity field noise that ensures the accurate application of the CVD algorithm to a universal flow field.Figure 5 illustrates the simulated vector field, with the background color map representing velocity magnitude.Additionally, Figure 6 displays the simulated velocity field after subtracting the freestream velocity, accompanied by a background color map representing vorticity.The resolution is quantified by the number of velocity vectors (n,m) spanning the diameter of the vortex, which corresponds to ∅ = 2 .This ratio,  ∅ ⁄ , defines the minimum necessary number of velocity vectors spanning a vortex's diameter for accurate The resolution is quantified by the number of velocity vectors (n,m) spanning the diameter of the vortex, which corresponds to ∅ v = 2r v .This ratio, nm/∅ v , defines the minimum necessary number of velocity vectors spanning a vortex's diameter for accurate identification and characterization.Figure 7 illustrates the vortex radius calculated by the CVD algorithm across different nm/∅ v values, alongside the known vortex radius employed in the Burgers model.Vortices with nm/∅ v < 5 are categorized as Type II errors (as described later).Considering the expected radius of r v = 3.18 mm, the computation of vortex radius for r v = 3.18 mm ± 5% is achieved when nm/∅ v > 17. Figure 7 also exhibits discretization artifacts originating from the simulation.
Figure 8 depicts the vortex circulation (Γ) calculated by the CVD algorithm across various nm/∅ v values.Circulation is not computed for cases with nm/∅ v < 5, as these vortices are rejected by the WA method.The circulation derived using the CVD algorithm pertains specifically to the core region circulation, representing the circulation induced within the vortex boundaries within −r v ≤ r ≤ r v .This value is expected to be smaller than Γ ∞ = 100 mm 2 /s, which corresponds to the circulation induced by the vortex at far distances (as r v → ∞ ).Using Γ ∞ = 100 mm 2 /s is impractical in this study due to interference from neighboring vortices, making it unobtainable experimentally.Measuring core region circulation is a more feasible approach to gauge vortex strength, as it can be calculated within a finite circle with a radius of r v , within which most of the vorticity is concentrated [53].Throughout this section, the term 'circulation Γ core ' refers to the core region circulation, and the expected core region circulation is determined as By substituting r v = 3.18 mm into Equation ( 19), the expected core circulation of the simulated vortex is determined to be Γ core = 71 mm 2 /s.Computation of core circulation within Γ core = 71 mm 2 /s ± 5% range is achieved when nm/∅ v > 25, as shown in Figure 8.For lower values of nm/∅ v , the circulation is consistently underestimated, primarily due to discretization errors in the summation of vorticity pixels within the vortex boundaries.The grid resolutions used for calculating vorticity radius (r v ) and core circulation (Γ core ) are summarized in Table 1 for reference., as these vortices are rejected by the WA method.The circulation derived using the CVD algorithm pertains specifically to the core region circulation, representing the circulation induced within the vortex boundaries within − ≤  ≤  .This value is expected to be smaller than Γ = 100 mm s ⁄ , which corresponds to the circulation induced by the vortex at far distances (as  → ∞).Using Γ = 100 mm s ⁄ is impractical in this study due to interference from neighboring vortices, making it unobtainable experimentally.Measuring core region circulation is a more feasible approach to gauge vortex strength, as it can be calculated within a finite circle with a radius of  , within which most of the vorticity is concentrated [53].Throughout this section, the term 'circulation Γ ' refers to the core region circulation, and the expected core region circulation is determined as By substituting  = 3.18 mm into Equation ( 19), the expected core circulation of the simulated vortex is determined to be  = 71 mm /s.Computation of core circulation within  = 71 mm /s ± 5% range is achieved when  ∅ > 25 ⁄ , as shown in Figure 8.For lower values of  ∅ ⁄ , the circulation is consistently underestimated, primarily due to discretization errors in the summation of vorticity pixels within the vortex boundaries.The grid resolutions used for calculating vorticity radius ( ) and core circulation (Γ ) are summarized in Table 1 for reference.The computed circulation pertains to the core region circulation, representing circulation induced within the vortex boundaries (− ≤  ≤  ).
Table 1.Grid resolution thresholds for simulated flow field.

Computation of 𝑟 within 5% of true value 𝑛𝑚 ∅ > 17 ⁄
Computation of Γ within 5% of true value  ∅ > 25 ⁄ Having established the minimum grid resolution requirements using an ideal simulated flow field, the robustness of the CVD algorithm to noise is now examined.Specifically, Gaussian white noise is deliberately introduced to the velocity vector components to determine the impact on vortex characterization.A simulated vector field with n = m = 200 was chosen, yielding a resolution of  ∅ = 43 ⁄ , which satisfies the minimum requirements for accurate computation of circulation and radius, as determined previously.Gaussian white noise is generated by adding random real number vectors to the x and y Table 1.Grid resolution thresholds for simulated flow field.

Threshold nm/∅ v
Type II error nm/∅ v < 5 Computation of r v within 5% of true value nm/∅ v > 17 Computation of Γ within 5% of true value nm/∅ v > 25 Having established the minimum grid resolution requirements using an ideal simulated flow field, the robustness of the CVD algorithm to noise is now examined.Specifically, Gaussian white noise is deliberately introduced to the velocity vector components to determine the impact on vortex characterization.A simulated vector field with n = m = 200 was chosen, yielding a resolution of nm/∅ v = 43, which satisfies the minimum requirements for accurate computation of circulation and radius, as determined previously.Gaussian white noise is generated by adding random real number vectors to the x and y components of each velocity vector, where both components follow a probability distribution with a zero mean and a predetermined variance.The noise vectors are statistically independent and have a standard deviation of ε σ times that of the vector fields.
The impact of adding white noise on vortex characterization is assessed by calculating the standard deviation and mean values of both the core region circulation and the detected vortex radii across a range of noise levels.For this study, 10 noise levels ranging from 0 ≤ ε σ ≤ 0.2 with 100 samples per noise level were investigated.Beyond ε σ = 0.2, the algorithm starts to fail in detecting the vortex.Figure 9 presents the mean circulation results, along with standard deviation error bars, and the expected circulation value for each of the 10 noise levels examined.Similarly, Figure 10 displays the mean vortex radius results for each of the noise levels.
Type II error  ∅ < 5 ⁄ Computation of  within 5% of true value  ∅ > 17 ⁄ Computation of Γ within 5% of true value  ∅ > 25 ⁄ Having established the minimum grid resolution requirements using an ideal simulated flow field, the robustness of the CVD algorithm to noise is now examined.Specifically, Gaussian white noise is deliberately introduced to the velocity vector components to determine the impact on vortex characterization.A simulated vector field with n = m = 200 was chosen, yielding a resolution of  ∅ = 43 ⁄ , which satisfies the minimum requirements for accurate computation of circulation and radius, as determined previously.Gaussian white noise is generated by adding random real number vectors to the x and y components of each velocity vector, where both components follow a probability distribution with a zero mean and a predetermined variance.The noise vectors are statistically independent and have a standard deviation of  times that of the vector fields.
The impact of adding white noise on vortex characterization is assessed by calculating the standard deviation and mean values of both the core region circulation and the detected vortex radii across a range of noise levels.For this study, 10 noise levels ranging from 0 ≤  ≤ 0.2 with 100 samples per noise level were investigated.Beyond  = 0.2, the algorithm starts to fail in detecting the vortex.Figure 9 presents the mean circulation results, along with standard deviation error bars, and the expected circulation value for each of the 10 noise levels examined.Similarly, Figure 10 displays the mean vortex radius results for each of the noise levels.The velocity vectors near the center of the vortex exhibit larger magnitudes compared to those near the boundary.Consequently, the impact of adding white noise is less pronounced near the core.However, increasing  eventually results in the entire vector field appearing incoherent.Notably, this incoherence starts from the outer regions and progresses inward due to the increasing velocity magnitude in the negative radial direction.Consequently, as  increases, the detection algorithm tends to underestimate the size of the vortex.It primarily identifies the center of the vortex, as this region is relatively less affected by the added noise (or is at least affected to a lesser extent, still resembling a vortex).Figure 11   The velocity vectors near the center of the vortex exhibit larger magnitudes compared to those near the boundary.Consequently, the impact of adding white noise is less pronounced near the core.However, increasing ε σ eventually results in the entire vector field appearing incoherent.Notably, this incoherence starts from the outer regions and progresses inward due to the increasing velocity magnitude in the negative radial direction.
Consequently, as ε σ increases, the detection algorithm tends to underestimate the size of the vortex.It primarily identifies the center of the vortex, as this region is relatively less affected by the added noise (or is at least affected to a lesser extent, still resembling a vortex).Figure 11 illustrates the probability density function (PDF) of the boundary radius of a detected vortex based on 100 simulations conducted with ε σ values of 0.05 and 0.1.For the higher noise level (ε σ = 0.1), a distribution with a negative bias in its mean relative to the expected value of 3.18 mm is observed, indicating that the algorithm tends to underestimate the size of detected vortices at higher noise levels, as anticipated.This highlights the robustness of the CVD algorithm under controlled conditions.The CVD algorithm was able to detect and importantly return the characteristics of a vortex defined by Burger's model.The introduction of noise to the velocity field did affect that detection and determination of the characteristic circulation but was only detrimental at high noise levels.
The velocity vectors near the center of the vortex exhibit larger magnitudes compared to those near the boundary.Consequently, the impact of adding white noise is less pronounced near the core.However, increasing  eventually results in the entire vector field appearing incoherent.Notably, this incoherence starts from the outer regions and progresses inward due to the increasing velocity magnitude in the negative radial direction.Consequently, as  increases, the detection algorithm tends to underestimate the size of the vortex.It primarily identifies the center of the vortex, as this region is relatively less affected by the added noise (or is at least affected to a lesser extent, still resembling a vortex).Figure 11 illustrates the probability density function (PDF) of the boundary radius of a detected vortex based on 100 simulations conducted with  values of 0.05 and 0.1.For the higher noise level ( = 0.1), a distribution with a negative bias in its mean relative to the expected value of 3.18 mm is observed, indicating that the algorithm tends to underestimate the size of detected vortices at higher noise levels, as anticipated.This highlights the robustness of the CVD algorithm under controlled conditions.The CVD algorithm was able to detect and importantly return the characteristics of a vortex defined by Burger's model.The introduction of noise to the velocity field did affect that detection and determination of the characteristic circulation but was only detrimental at high noise levels.

Verification on an Experimental Flow Field
The velocity flow field immediately downstream from an oscillating NACA 0012 airfoil is captured using a PIV system with the experimental setup described below.The CVD approach is used to identify and characterize vortical structures.

Experimental Facility and PIV Setup
The experimental facility was a water flume measuring 0.7 × 0.4 m (27.5 ′′ × 16 ′′ ), designed to exhibit low turbulence characteristics as detailed by Hilderman (2004) [54].Within this setup, an airfoil is vertically suspended, allowing it to extend through the free surface into the water channel, positioned perpendicular to the upstream flow direction.The experimental setup and subsequent analysis were conducted in a two-dimensional domain, focusing on the planar cross-section of the wing's wake.
The experimental PIV setup, depicted in Figure 12, was equipped with four 2112 × 2072 pixel resolution, 14-bit, dual-frame CCD cameras (Imager Pro X 4M, LaVision).These cameras view the investigation plane through four independent local coordinate systems.Calibration of the cameras was carried out using a 300 mm × 800 mm calibration target featuring 1.3 mm diameter markers spaced at 3 mm intervals.Subsequently, the captured images were de-warped to account for variations in camera viewing angles and were stitched together with specified offsets to create a unified global field.
To mitigate surface refraction effects, an acrylic sheet was positioned at the free surface, with the cameras capturing images through it.In the illumination process, a double-pulse Nd:YAG laser (Solo III-15Hz, New Wave Research Fremont, California, United States) was employed to illuminate the flow.This flow was seeded with 18µm hollow glass spheres (SPHERICEL ® , Potters Industries, Malvern, PA, USA).The laser beam was focused into a thin sheet and directed upstream by a mirror and a submerged periscope.
Furthermore, the measurement error associated with determining the precise location of the correlation peak for an 8-pixel particle displacement is approximately 1-2% of the full-scale velocity for similar planar PIV systems [56].Lastly, in terms of event timing accuracy, the system achieves a resolution of 10 ns with a jitter of less than 1 ns.
In summary, the key measurement accuracy specifications of the present PIV system include a particle slip velocity error of 0.47% relative to the freestream velocity, 0.3% uncertainty in image magnification, 1-2% error in determining particle displacement correlation peaks, and 10 ns timing resolution with sub-ns jitter.Before performing cross-correlation on image pairs, the raw images underwent preprocessing using commercial software (LaVision GmbH, DaVis 8.05).This preprocessing step was essential as it enhances particle intensity and shape, ultimately resulting in an improved correlation peak [55].The strength of correlation can be affected by disparities in image intensity caused by factors such as non-uniform light sheets, shadows, reflections, or variations in particle size [55].To enhance correlation strength, several image preprocessing methods were employed, including background intensity subtraction, sliding minimum subtraction, and min-max filtering for intensity normalization.In this study, the generation of the vector field involved three passes of a cross-correlation algorithm with window shifting.The first pass used a 64×64 pixel interrogation window, while the subsequent two passes employed a 32 ×32 pixel interrogation window with a 50% overlap.
An important consideration in this study is the error that arises from the disparity between the motion of seed particles (r) and the actual fluid motion.This error is primarily attributed to particle slip, wherein the seed particles lag behind the fluid motion by a finite quantity.To compute the slip velocity to first order, the approach outlined by Adrian and Westerweel (2011) [56] was utilized: In this equation, v p represents the velocity of the seed particle, v f is the fluid velocity, and g is the gravitational constant with a value of 9.81 m/s 2 .The time constant, τ 0 , is defined as where d p is the diameter of the seed particles, set to 18 µm, and the densities of the seed particle and water are ρ p = 600 kg/m 3 and ρ f = 998 kg/m 3 , respectively.Additionally, the density ratio, ρ , is defined as The slip velocity error, which amounts to 0.47% relative to the freestream velocity (U ∞ = 0.017 mm/s), is initially approximated using Equation (20).To address variations in image magnification across the image domain, calibration with a target is employed.However, it is important to note that image magnification also varies along the thickness of the light sheet, introducing a magnification uncertainty typically around 0.3% in most PIV setups, as discussed by Adrian and Westerweel (2011) [56].
Furthermore, the measurement error associated with determining the precise location of the correlation peak for an 8-pixel particle displacement is approximately 1-2% of the full-scale velocity for similar planar PIV systems [56].Lastly, in terms of event timing accuracy, the system achieves a resolution of 10 ns with a jitter of less than 1 ns.
In summary, the key measurement accuracy specifications of the present PIV system include a particle slip velocity error of 0.47% relative to the freestream velocity, 0.3% uncertainty in image magnification, 1-2% error in determining particle displacement correlation peaks, and 10 ns timing resolution with sub-ns jitter.
The aluminum airfoil used was a continuous extrusion featuring a cross-sectional shape conforming to the NACA 0012 airfoil profile, with a chord length (C) of 75 mm.The airfoil's cross-section and the parameters governing its motion are visually represented in Figure 13.The airfoil is suspended on a sturdy shaft, which is driven by a stepper motor (PK258-02Dl, ORIENTAL MOTOR CO.LTD., Tokyo, Japan).This shaft passes through the airfoil's aerodynamic center, enabling precise pitch oscillations of any desired waveform.The aluminum airfoil used was a continuous extrusion featuring a cross-sectional shape conforming to the NACA 0012 airfoil profile, with a chord length (C) of 75 mm.The airfoil's cross-section and the parameters governing its motion are visually represented in Figure 13.The airfoil is suspended on a sturdy shaft, which is driven by a stepper motor (PK258-02Dl, ORIENTAL MOTOR CO.LTD., Tokyo, Japan).This shaft passes through the airfoil's aerodynamic center, enabling precise pitch oscillations of any desired waveform.
To orchestrate these oscillations and maintain control, a commercial hardware system (DS1104, dSPACE Inc., Paderborn, Germany) was programmed.This system governs the motion of the stepper motor and generates output trigger signals when the airfoil's pitch reaches the desired angle.Consequently, this setup synchronizes the imaging system to capture data at predefined positions, facilitating the acquisition of both phase-averaged and time-averaged datasets.To maintain a well-controlled flow and avoid the formation of disruptive leadingedge vortices in the wake, the airfoil, characterized by its chord length (C), is limited to small pitch oscillation amplitudes ( ≤ 8°).The flow can be characterized by the Reynolds number (Re), defined in terms of the airfoil's chord thickness (D) as  To orchestrate these oscillations and maintain control, a commercial hardware system (DS1104, dSPACE Inc., Paderborn, Germany) was programmed.This system governs the motion of the stepper motor and generates output trigger signals when the airfoil's pitch reaches the desired angle.Consequently, this setup synchronizes the imaging system to capture data at predefined positions, facilitating the acquisition of both phase-averaged and time-averaged datasets.
To maintain a well-controlled flow and avoid the formation of disruptive leading-edge vortices in the wake, the airfoil, characterized by its chord length (C), is limited to small pitch oscillation amplitudes (θ A ≤ 8 • ).The flow can be characterized by the Reynolds number (Re), defined in terms of the airfoil's chord thickness (D) as In the experiments conducted in this study, a constant freestream velocity (U ∞ ) of 17 mm/s was maintained, using an airfoil thickness (D) of 8.6 mm and a kinematic viscosity of water (ν) equal to 1 × 10 −6 m 2 /s, resulting in a Reynolds number (Re) of 146.
The oscillation waveform of the airfoil is given by By adjusting the oscillation amplitude (θ A ) and frequency (f ), various wake conditions can be achieved.

Analysis of Oscillating Airfoil Wake
A schematic depicting the oscillating airfoil and the flow features observed in its wake is illustrated in Figure 14, highlighting the relationship between the airfoil and the generated vortices.The use of small airfoil oscillation amplitudes results in an orderly wake pattern, characterized by the shedding of precisely two counter-rotating vortices per oscillation cycle, as detailed by Bohl and Koochesfahani (2009) [2].The airfoil, with a chord length C, is shown in the context of the uniform flow with a velocity U ∞ .
These vortices are organized into two distinct rows separated by a distance S y and aligned with the flow direction.Within each row, the vortices are spaced apart by S x .The core coordinates of the vortices are labeled, and a core boundary radius (r v ) is defined.
Additionally, the drift velocity, ⇀ v dri f t = v dri f t,x , v dri f t,y , is represented by an instantaneous velocity vector at the grid point coinciding with the vortex core.Furthermore, the vortices are characterized by their peak vorticity (ω peak ) and circulation (Γ).The circulation is approximated by summing the vorticity at each velocity measurement location determined by PIV within a radius r v of a vortex as follows: Here, s refers to the surface integral over area ds, and A pixel is the surface area of a rectangle formed from the coordinates of 4 adjacent velocity vectors.
By manipulating the parameters θ A and f, various wake configurations can be achieved.Common wake patterns for sinusoidal pitching symmetric airfoils include the following [57][58][59][60]: 1.
von Karman Wake: This configuration, illustrated in Figure 15a, resembles the wake pattern typically associated with vortex shedding from a cylinder.

2.
Aligned Wake: In this arrangement, the transverse separation distance S y between vortices approaches zero.

3.
Inverted von Karman Wake: This wake pattern, illustrated in Figure 15b, differs from the von Karman wake by reversing the rotation direction of each row of vortices, as described by Godoy-Diana et al. ( 2009) [10].
At an airfoil pitch angle of 0 • , instantaneous velocity vector fields are computed for 100 wake datasets using PIV.These fields are then averaged to generate a phase-averaged representation since all the data are collected at a fixed oscillation phase.Subsequently, WA and CSL methods are applied to the phase-averaged fields, and the results are depicted in Figure 15.Within Figure 15, the CSL method identifies vortex core locations, represented as white crosses, and vortex boundary radii, denoted by black circles, superimposed on the vorticity field.
oscillation cycle, as detailed by Bohl and Koochesfahani (2009) [2].The airfoil, with a chord length C, is shown in the context of the uniform flow with a velocity  .
These vortices are organized into two distinct rows separated by a distance Sy and aligned with the flow direction.Within each row, the vortices are spaced apart by Sx.The core coordinates of the vortices are labeled, and a core boundary radius ( ) is defined.
Additionally, the drift velocity,  ⃑ =  , ,  , , is represented by an instantaneous velocity vector at the grid point coinciding with the vortex core.Furthermore, the vortices are characterized by their peak vorticity ( ) and circulation (Γ).The circulation is approximated by summing the vorticity at each velocity measurement location determined by PIV within a radius  of a vortex as follows: Here, s refers to the surface integral over area ds, and  is the surface area of a rectangle formed from the coordinates of 4 adjacent velocity vectors.By manipulating the parameters  and f, various wake configurations can be achieved.Common wake patterns for sinusoidal pitching symmetric airfoils include the following [57][58][59][60]: Figure 15a illustrates the von Karman wake pattern, resulting from sinusoidal pitching at f = 1.6 rad/s and θ A = 8 • , while Figure 15b showcases the inverted von Karman wake pattern generated by sinusoidal pitching at f = 3.4 rad/s and θ A = 8 • .The von Karman wake typically displays weak vortices characterized by low peak vorticity and circulation.It exhibits minimal vortex decay over downstream distances and features larger streamwise spacing (S x ).Conversely, the inverted von Karman wake shown in Figure 15b exhibits vortices with a higher peak vorticity and circulation.These vortices also experience the most rapid decay as they progress downstream and display relatively smaller streamwise spacing (S x ).These two wake configurations at Re = 146 exemplify the extreme cases of wake vorticity for θ A = 8 • .circulation.It exhibits minimal vortex decay over downstream distances and features larger streamwise spacing (Sx).Conversely, the inverted von Karman wake shown in Figure 15b exhibits vortices with a higher peak vorticity and circulation.These vortices also experience the most rapid decay as they progress downstream and display relatively smaller streamwise spacing (Sx).These two wake configurations at  = 146 exemplify the extreme cases of wake vorticity for  = 8°.A false positive vortex diagnosis, or a Type I error, occurs when a vortex is identified where none exists [35], potentially due to local shear, boundary influences, or image stitching inaccuracies.Figure 15 highlights this issue with the CSL method, which may incorrectly confirm a ROI containing no vortex.However, when the CSL data are processed through the WA algorithm, these false positives are eliminated, as it checks for streamlines forming semi-closed, semi-elliptical paths, which are absent in the false detections.
The Supplementary Material includes animations (Supplemental Files S1 and S2) featuring instantaneous velocity vector fields for the von Karman Wake and inverted von Karman Wake cases.These animations dynamically illustrate vortex core locations and vortex boundary radii superimposed on the vorticity field for each wake configuration.

Sample Experimental Vortex for WA and CSL Verification
Figure 16 provides a detailed view of a vortex field sample from Figure 15a, highlighted within a red square, and includes a locally calculated velocity field centered on the vortex core.This calculation employs a moving reference frame that matches the drift velocity of the sample vortex, as described in the pre-WA section of the combinatorial algorithm.
WA algorithm does not correct these errors directly.Instead, it acts as an indicator, and addressing Type II errors requires additional processing, such as using clustering algorithms.The primary vortex is determined by the largest number of complying streamlines, while other streamline clusters are marked as potential Type II errors.This underscores that while the WA algorithm detects deviations, further analyses and processing are necessary to correct and characterize potentially missed vortices within the ROI.In Figure 16a, the streamlines generated from this velocity field are plotted over the local vorticity field.Figure 16b highlights five streamlines that meet the criteria set by the WA method.Each of these compliant streamlines is represented by a point (depicted as black dots), and the vortex core, previously determined using the CSL method, is indicated with a black '+' symbol.The WA algorithm serves to confirm the presence of a vortex within the sample ROI, consistent with the vortex definition of Robinson (1991) [36].It also verifies the existence of a single vortex within the sample ROI and corroborates that the core location, computed by CSL, closely aligns with the geometric center of all the complying streamlines.Should the WA algorithm identify multiple vortices within a single ROI, it signals the possibility of Type II errors.However, it is crucial to note that the WA algorithm does not correct these errors directly.Instead, it acts as an indicator, and addressing Type II errors requires additional processing, such as using clustering algorithms.The primary vortex is determined by the largest number of complying streamlines, while other streamline clusters are marked as potential Type II errors.This underscores that while the WA algorithm detects deviations, further analyses and processing are necessary to correct and characterize potentially missed vortices within the ROI.
A false negative, or a Type II error, occurs when a vortex exists in the flow but remains undetected by the algorithm.Unfortunately, due to the a priori knowledge required by the WA algorithm, it is inevitable to encounter Type II errors for vortices that fall outside the predefined ROI.However, when multiple vortices coexist within a single ROI, the WA algorithm has the capability to identify a Type II error, but only if the drift velocity of the primary vortex closely matches that of any secondary vortices represented by separate streamline clusters.
The proposed CVD algorithm offers ways to minimize Type II errors through two key mechanisms: • Wider Threshold Bounds: By expanding the threshold bounds in the threshold intensity vector (TIV) defined in Equation (1) for the MV algorithm, the ROI can encompass weaker vortices, including those with lower vorticity.This adjustment ensures that such vortices are considered by both the CSL and WA algorithms.

•
Erosion Process Optimization: Employing two smaller sizes of image morphology (IM) structuring elements in the erosion process preserves smaller vortices, those with small values of (r v ), allowing them to remain within the ROI.These vortices can then be assessed by the CSL and WA algorithms.
It is important to note that implementing these modifications significantly increases the computational time.Therefore, it is recommended to define a minimum vortex size r v,min and strength ω peak,min before running the algorithm.Weak and/or small vortices falling below the specified threshold are intentionally disregarded, reducing the algorithm's computational demands.The selection of both the threshold defined in Equation ( 1) and the sizes of IM structuring elements (S e1 , S e2 , and S e3 ) are based on the goal of minimizing Type II errors while maintaining efficient computation.

Vorticity Distribution Profiles
With the vortex identified and core location determined, the distribution of vorticity within the vortex can be examined.The vorticity profiles of these vortices play a crucial role in assessing various aspects of vortex characteristics, including the accuracy of the vortex radius (r v ), vortex shape and symmetry and their suitability for fitting analytical models.
In Figure 17, the vorticity distribution of the sample vortex depicted in Figure 16 is plotted against the dimensionless radius (y − y c )/r v along a constant y line traversing the vortex core.The vorticity profile indicates a Gaussian distribution with slight asymmetry, which may be attributed to interactions with nearby counter-rotating vortices.
(IM) structuring elements in the erosion process preserves smaller vortices, those with small values of ( ), allowing them to remain within the ROI.These vortices can then be assessed by the CSL and WA algorithms.
It is important to note that implementing these modifications significantly increases the computational time.Therefore, it is recommended to define a minimum vortex size  , and strength  , before running the algorithm.Weak and/or small vortices falling below the specified threshold are intentionally disregarded, reducing the algorithm's computational demands.The selection of both the threshold defined in Equation (1) and the sizes of IM structuring elements (Se1, Se2, and Se3) are based on the goal of minimizing Type II errors while maintaining efficient computation.

Vorticity Distribution Profiles
With the vortex identified and core location determined, the distribution of vorticity within the vortex can be examined.The vorticity profiles of these vortices play a crucial role in assessing various aspects of vortex characteristics, including the accuracy of the vortex radius ( ), vortex shape and symmetry and their suitability for fitting analytical models.
In Figure 17, the vorticity distribution of the sample vortex depicted in Figure 16 is plotted against the dimensionless radius ( −  )  ⁄ along a constant  line traversing the vortex core.The vorticity profile indicates a Gaussian distribution with slight asymmetry, which may be attributed to interactions with nearby counter-rotating vortices.Figure 17 also includes vorticity profiles for ideal Burgers and Rankine vortices.These profiles have the same peak vorticity and radii as those computed by the CSL algorithm.The Burgers vortex was described in Equations ( 17) and ( 18) and the Rankine vortex represents a simplified model that attempts to mimic real vortices by dividing them into two regions: the inner core with uniform vorticity, resembling a forced vortex, and the outer core, which lacks vorticity, simulating an irrotational or free vortex [61].In contrast, the Burgers vortex solution is a more intricate model used to illustrate fundamental elements of modern turbulence theory, providing an exact solution to the cylindrical Navier-Stokes equations that accounts for flow on a cylindrical vortex core inducing circulation (Γ ∞ ) at large distances [51].
The vorticity distribution plot in Figure 17 reveals that the sample vortex profile closely resembles a Burgers vortex, sharing the same radius and peak vorticity values computed by the CSL algorithm.This correspondence validates the CSL algorithm's effectiveness in predicting vortex radius and core coordinates from a PIV-generated velocity vector field.A slight spatial offset between the experimental vortex and the Burgers vortex is attributed to the spatial resolution of the experimental PIV data.This discrepancy arises from the CSL algorithm's sequential reading of v p (ς, ψ) in Equation ( 4), from left to right.When the true core coordinate falls between two velocity vectors, the algorithm consistently selects the leftmost value.

Circumferential Velocity Profiles
In addition to vorticity, circumferential velocity profiles are valuable for characterizing vortical flow structures.The circumferential, or azimuthal, velocity represents the velocity component perpendicular to any straight line passing through the vortex core.
Figure 18 provides a comparison of the circumferential velocity profile for the sample vortex from Figure 16 against theoretical profiles for a Rankine vortex and Burgers vortex.The sample vortex profile is plotted alongside the analytical curves.As observed in Figure 17, defining the vortex boundary radius based solely on vorticity profiles requires an arbitrary cutoff [1].This highlights the limitations of relying solely on the MV method for radius determination.However, circumferential velocity profiles offer a distinct advantage for defining the boundary.Within the vortex radius, the circumferential velocity magnitude generally increases with radial distance, reaching a maximum value precisely at the boundary (r v ) [1].This characteristic peak at r v is clearly illustrated for the sample vortex in Figure 18.Using the peak velocity avoids ambiguity and arbitrary cutoffs.Additionally, circumferential velocity can be measured farther from the high shear at the core, where PIV more accurately represents curvature.Thus, circumferential profiles provide an unambiguous vortex boundary definition without relying on vorticity gradients alone.

Circulation Analysis in the Von Karman Wake
With the CVD approach now in place, it can be used to analyze a large number of PIV datasets without the need for setting limits to individual vector sets to detect vortices.The flow field in Figure 15a has some unique characteristics.In the near wake, vortices exhibit more consistent spatial locations compared to those further downstream.In other words, vortex cores close to the trailing edge of the flapping wing consistently appear in the same positions across oscillation periods.However, structures farther downstream demonstrate greater spatial variation.Phase averaging multiple vector fields reduces the magnitudes of spatially varying vortices while preserving more consistent structures.The CVD method enables quantifying key wake parameters after computing vorticity distributions and defining vortex boundaries and cores.As highlighted before, these parameters include circulation (Γ), peak vorticity ( ), vortex radii ( ), vortex drift velocity ( ), streamwise spacing (Sx), and transverse spacing (Sy).Figure 19 plots the circulation of 25 instantaneous vector fields from the case illustrated in Figure 15a against dimensionless downstream distance when the airfoil angle is  = 0°.Individual field circulations from the instantaneous images appear as black dots, while phase-averaged circulation is shown in red.This illustrates the increasing variability of the vortex core streamwise coordinate  with downstream distance.In the near field  =   ⁄ < 2, tight vertical data groups indicate consistent vortex locations.However, groups become indistinguishable farther downstream as locations vary.Consequently, phase averaging underestimates circulation versus individual fields.
Circulation provides critical insights into vortex strength and influence within the wake.Examining circulation across instantaneous fields reveals vortex evolution and interactions over time, granting a deeper understanding of von Karman wake dynamics.Circumferential velocity profiles effectively showcase the CSL algorithm's ability to distinctly identify the boundary of a vortex core by pinpointing the radial location where the absolute circumferential velocity is maximized.This method of defining the core boundary radius is less prone to ambiguity compared to one derived from vorticity profiles, making it a superior approach for computing vortex boundary radii.

Circulation Analysis in the Von Karman Wake
With the CVD approach now in place, it can be used to analyze a large number of PIV datasets without the need for setting limits to individual vector sets to detect vortices.The flow field in Figure 15a has some unique characteristics.In the near wake, vortices exhibit more consistent spatial locations compared to those further downstream.In other words, vortex cores close to the trailing edge of the flapping wing consistently appear in the same positions across oscillation periods.However, structures farther downstream demonstrate greater spatial variation.Phase averaging multiple vector fields reduces the magnitudes of spatially varying vortices while preserving more consistent structures.The CVD method enables quantifying key wake parameters after computing vorticity distributions and defining vortex boundaries and cores.As highlighted before, these parameters include circulation (Γ), peak vorticity (ω peak ), vortex radii (r v ), vortex drift velocity (v dri f t ), streamwise spacing (S x ), and transverse spacing (S y ).
Figure 19 plots the circulation of 25 instantaneous vector fields from the case illustrated in Figure 15a against dimensionless downstream distance when the airfoil angle is θ a f = 0 • .Individual field circulations from the instantaneous images appear as black dots, while phase-averaged circulation is shown in red.This illustrates the increasing variability of the vortex core streamwise coordinate x C with downstream distance.In the near field x C = x/C < 2, tight vertical data groups indicate consistent vortex locations.However, groups become indistinguishable farther downstream as locations vary.Consequently, phase averaging underestimates circulation versus individual fields.

Considerations for Higher Reynolds Number Applications
At the demonstrated Reynolds number of 146 based on airfoil thickness and freestream velocity, the flow remains laminar and ordered, facilitating structured vortex shedding amenable to characterization with the combinatorial vortex detection (CVD) technique.However, applying the algorithm to scenarios with significantly higher Re introduces additional complexities that must be considered.
As Re increases in the transitional flow regime, turbulence levels intensify.This places greater demand on spatial resolution to fully resolve steeper velocity gradients and smaller-scale flow structures.If spatial resolution becomes insufficient, discretization errors may prevent accurate calculation of vorticity and vortex parameters.Additionally, heightened turbulence increases variability in vortex properties such as circulation, boundary radius, and peak vorticity.To achieve a statistically reliable description of vortex characteristics, an increased number of velocity field samples is necessary.Lastly, elevated turbulence can cause vortex core boundaries to become more diffuse [62,63], complicating boundary identification.Addressing these challenges would necessitate modifications such as improving localized velocity calculations and using an adaptive ROI threshold approach.
Proper validation of the CVD algorithm at higher Re requires a more controlled experimental facility and PIV setup meeting minimum resolution requirements.While not explored presently with the available resources, the proposed method shows promise for characterization of transitional flows.Future efforts should focus on method refinements targeting higher Re applications.

Conclusions
In this study, a robust vortex detection and characterization algorithm has been developed and rigorously tested on 2C2D velocity vector fields obtained from the wake of an oscillating NACA 0012 airfoil within a uniform flow.Existing vortex detection methods face major limitations in handling experimental data, including sensitivity to measurement noise and substantial velocity gradients near vortex cores that complicate particle seeding and the performance of the PIV algorithm [12,39].
The proposed combinatorial vortex detection (CVD) algorithm aimed to address these challenges and accurately identify wake vortices and their attributes.It harnesses Circulation provides critical insights into vortex strength and influence within the wake.Examining circulation across instantaneous fields reveals vortex evolution and interactions over time, granting a deeper understanding of von Karman wake dynamics.

Considerations for Higher Reynolds Number Applications
At the demonstrated Reynolds number of 146 based on airfoil thickness and freestream velocity, the flow remains laminar and ordered, facilitating structured vortex shedding amenable to characterization with the combinatorial vortex detection (CVD) technique.However, applying the algorithm to scenarios with significantly higher Re introduces additional complexities that must be considered.
As Re increases in the transitional flow regime, turbulence levels intensify.This places greater demand on spatial resolution to fully resolve steeper velocity gradients and smallerscale flow structures.If spatial resolution becomes insufficient, discretization errors may prevent accurate calculation of vorticity and vortex parameters.Additionally, heightened turbulence increases variability in vortex properties such as circulation, boundary radius, and peak vorticity.To achieve a statistically reliable description of vortex characteristics, an increased number of velocity field samples is necessary.Lastly, elevated turbulence can cause vortex core boundaries to become more diffuse [62,63], complicating boundary identification.Addressing these challenges would necessitate modifications such as improving localized velocity calculations and using an adaptive ROI threshold approach.
Proper validation of the CVD algorithm at higher Re requires a more controlled experimental facility and PIV setup meeting minimum resolution requirements.While not explored presently with the available resources, the proposed method shows promise for characterization of transitional flows.Future efforts should focus on method refinements targeting higher Re applications.

Figure 1 .
Figure 1.Schematic illustrating the application of a multilevel threshold technique for identifying vortices in two vorticity fields.(left): indexed images with levels i = 1, 2, and 3. (right): corresponding desired ROI maps for vortex extraction.

Figure 2 .
Figure 2. Schematic representation of the cross-sectional lines (CSL) procedure with four concentric ellipses illustrating vorticity contours and the vortex core indicated by an '×'.

Figure 2 .
Figure 2. Schematic representation of the cross-sectional lines (CSL) procedure with four concentric ellipses illustrating vorticity contours and the vortex core indicated by an '×'.

Figure 3 .
Figure 3. Computation of angle  , on streamline Sk for the winding angle algorithm.

Figure 3 .
Figure 3. Computation of angle α w,i on streamline S k for the winding angle algorithm.

Fluids 2024, 9 , 31 Figure 4 .
Figure 4. Flowchart of the combinatorial vortex detection (CVD) algorithm, utilizing three detection methods, MV, CSL and WA (green boxes), to process a 2C2D velocity vector field and generate characterized vortex data.

Figure 4 .
Figure 4. Flowchart of the combinatorial vortex detection (CVD) algorithm, utilizing three detection methods, MV, CSL and WA (green boxes), to process a 2C2D velocity vector field and generate characterized vortex data.

Figure 5 .
Figure 5. Simulated velocity vector field (200 × 200 resolution) with velocity magnitude represented by a color map, showing every 6th vector.

Figure 6 .Figure 5 .
Figure 6.Simulated velocity vector field (200 × 200 resolution) with the freestream velocity ( ) subtracted, displaying every 6th vector, and featuring vorticity as a color map in the background.The resolution is quantified by the number of velocity vectors (n,m) spanning the diameter of the vortex, which corresponds to ∅ = 2 .This ratio,  ∅ ⁄ , defines the minimum necessary number of velocity vectors spanning a vortex's diameter for accurate

Figure 5 .
Figure 5. Simulated velocity vector field (200 × 200 resolution) with velocity magnitude represented by a color map, showing every 6th vector.

Figure 6 .
Figure 6.Simulated velocity vector field (200 × 200 resolution) with the freestream velocity ( ) subtracted, displaying every 6th vector, and featuring vorticity as a color map in the background.

Figure 6 .
Figure 6.Simulated velocity vector field (200 × 200 resolution) with the freestream velocity (U ∞ ) subtracted, displaying every 6th vector, and featuring vorticity as a color map in the background.

Fluids 2024, 9 ,
x FOR PEER REVIEW 15 of 31 identification and characterization.Figure 7 illustrates the vortex radius calculated by the CVD algorithm across different  ∅ ⁄ values, alongside the known vortex radius employed in the Burgers model.Vortices with  ∅ < 5 ⁄ are categorized as Type II errors (as described later).Considering the expected radius of  = 3.18 mm, the computation of vortex radius for  = 3.18 mm ± 5% is achieved when  ∅ > 17 ⁄ .Figure 7 also exhibits discretization artifacts originating from the simulation.

Figure 7 .
Figure 7. Vortex radius variations calculated by the CVD algorithm for different  ∅ ⁄ values, compared with the expected radius ( = 3.18 mm ± 5% ) in the Burgers model.Vortices with  ∅ ⁄ < 5 are disregarded.Discretization artifacts from the simulation are also evident.

Figure 8
Figure 8 depicts the vortex circulation (Γ) calculated by the CVD algorithm across various  ∅ ⁄ values.Circulation is not computed for cases with  ∅ < 5 ⁄, as these vortices are rejected by the WA method.The circulation derived using the CVD algorithm pertains specifically to the core region circulation, representing the circulation induced within the vortex boundaries within − ≤  ≤  .This value is expected to be smaller than Γ = 100 mm s ⁄ , which corresponds to the circulation induced by the vortex at far

Figure 7 .
Figure 7. Vortex radius variations calculated by the CVD algorithm for different nm/∅ v values, compared with the expected radius (r v = 3.18 mm ± 5%) in the Burgers model.Vortices with nm/∅ v < 5 are disregarded.Discretization artifacts from the simulation are also evident.Fluids 2024, 9, x FOR PEER REVIEW 16 of 31

Figure 8 .
Figure 8. Circulation data are excluded for cases with  ∅ ⁄ < 5, as rejected by the WA method.

Figure 8 .
Figure 8. Circulation data are excluded for cases with nm/∅ v < 5, as rejected by the WA method.The computed circulation pertains to the core region circulation, representing circulation induced within the vortex boundaries (−r v ≤ r ≤ r v ).

Figure 9 .
Figure 9. Mean circulation computed by the CVD for 100 simulations at various noise levels, along with standard deviation error bars.

Fluids 2024, 9 , 31 Figure 9 .
Figure 9. Mean circulation computed by the CVD for 100 simulations at various noise levels, along with standard deviation error bars.

Figure 10 .
Figure 10.Mean radius computed by the CVD for 100 simulations at various noise levels, along with standard deviation error bars.
illustrates the probability density function (PDF) of the boundary radius of a detected vortex based on 100 simulations conducted with  values of 0.05 and 0.1.For r (mm)

Figure 10 .
Figure 10.Mean radius computed by the CVD for 100 simulations at various noise levels, along with standard deviation error bars.

Figure 12 .
Figure 12.Experimental setup with a 0.7 × 0.4 m water channel designed for low turbulence.It includes a vertically suspended airfoil extending into the channel and a PIV setup with four CCD cameras capturing the measurement plane.

Figure 12 .
Figure 12.Experimental setup with a 0.7 × 0.4 m water channel designed for low turbulence.It includes a vertically suspended airfoil extending into the channel and a PIV setup with four CCD cameras capturing the measurement plane.

Figure 13 .
Figure 13.Representation of the aluminum airfoil, a continuous extrusion with an NACA 0012 cross-sectional shape and a chord length (C) of 75 mm.The figure also illustrates the motion parameters, including pitch oscillations, enabled by a sturdy shaft driven by a stepper motor passing through the airfoil's aerodynamic center.

Figure 13 .
Figure 13.Representation of the aluminum airfoil, a continuous extrusion with an NACA 0012 crosssectional shape and a chord length (C) of 75 mm.The figure also illustrates the motion parameters, including pitch oscillations, enabled by a sturdy shaft driven by a stepper motor passing through the airfoil's aerodynamic center.

Figure 14 .
Figure 14.Schematic representation of the oscillating airfoil and the associated wake flow features.The airfoil, with chord length C, interacts with a uniform flow ( ), generating organized vortices.These vortices form two distinct rows, separated by Sy and spaced apart by Sx, aligned with the flow direction.Core coordinates and the core boundary radius ( ) are identified.The drift velocity is represented by an instantaneous vector at the vortex core grid point.

Figure 14 .
Figure 14.Schematic representation of the oscillating airfoil and the associated wake flow features.The airfoil, with chord length C, interacts with a uniform flow (U ∞ ), generating organized vortices.These vortices form two distinct rows, separated by S y and spaced apart by S x , aligned with the flow direction.Core coordinates and the core boundary radius (r v ) are identified.The drift velocity is represented by an instantaneous vector at the vortex core grid point.

Figure 15 .
Figure 15.Phase-averaged results obtained at an airfoil pitch angle of 0 degrees using PIV data.The cross-section lines (CSL) method identifies vortex core locations (white crosses) and vortex boundary radii (black circles), superimposed on the vorticity field; (a) showcases the von Karman wake pattern generated by sinusoidal pitching at f = 1.6 rad/s and θ A = 8 • , while (b) illustrates the inverted von Karman wake pattern resulting from sinusoidal pitching at f = 3.4 rad/s and θ A = 8 • .

Figure 16 .
Figure 16.Illustration of a sample vortex along with a locally computed velocity field centered on the vortex core.The velocity field is calculated in a moving reference frame matching the drift velocity of the sample vortex.In (a), streamlines generated from the velocity field overlay the local vorticity field, and (b) highlights five streamlines meeting WA criteria, along with cluster points (black), determined via the CSL method, denoted by a black '+' symbol.

Figure 17 .
Figure 17.Vorticity distribution profiles of an identified vortex, showing the vorticity plotted against the dimensionless radius along a constant y line traversing the vortex core.Additionally, vorticity profiles for ideal Burgers and Rankine vortices are included, sharing the same peak vorticity and radii as computed by the CSL algorithm.

Fluids 2024, 9 , 31 Figure 18 .
Figure 18.Comparison of circumferential velocity profiles for the sample vortex, a Rankine vortex, and a Burgers vortex.

Figure 18 .
Figure 18.Comparison of circumferential velocity profiles for the sample vortex, a Rankine vortex, and a Burgers vortex.

Figure 19 .
Figure 19.Circulation vs. dimensionless downstream distance for 25 instantaneous fields (black dots) and phase-averaged data (red circles) when the airfoil is at θ a f = 0 • .Near field (x C < 2) shows tight vertical data grouping, indicating consistent vortex positions.