This study explores the dynamics and statistical patterns of coherent long-lived vortices spontaneously forming in bluff body wakes. The analysis is based on a series of two-dimensional direct numerical simulations performed for a wide range of Reynolds numbers. We demonstrate that the majority of coherent vortices beyond the recirculation zone are well represented by the canonical Lamb–Oseen solution. This observation is used to develop a low-order census of long-lived eddies in terms of their core sizes and vorticity magnitudes. We demonstrate that the increase in the Reynolds number (Re) leads to the systematic reduction in the initial core radii (*r*_{0}), whereas the core vorticity (*ζ*_{0}) increases. These dependencies exhibit singular behavior in the inviscid limit (Re → ∞), which is captured by the proposed explicit relations for *r*_{0}(Re) and *ζ*_{0}(Re).

## I. INTRODUCTION

One of the most fundamental and ubiquitous properties of large-scale flows impinging on rigid obstacles is their tendency to generate a row of coherent vortices, commonly referred to as the von Kármán vortex street. von Kármán vortices have been the subject of strong and persistent interest due to the role they play in numerous engineering (Green and Unruh, 2006; Platzer *et al.*, 2008; and Demori *et al.*, 2017), geophysical (Dong *et al.*, 2007; Liu and Chang, 2018), and biological (Ortega-Jimenez *et al.*, 2013; 2014; and Wieskotten *et al.*, 2011) applications. For instance, the aerial (Muller *et al.*, 2015) and satellite-based (Young and Zawislak, 2006; Ito and Niino, 2015) observations present some of the most spectacular examples of vortex streets—in both the ocean and the atmosphere—downstream of small islands. Another series of studies focus on vortices generated by propagating solid objects (Voropayev and Smirnov, 2003; Afanasyev, 2004; Chongsiripinyo *et al.*, 2017; and Radko and Lorfeld, 2018). These efforts are motivated by the growing interest in developing hydrodynamically-based methods for the tracking and categorization of wakes. Our current understanding of the physical properties of von Kármán vortices is largely gained from laboratory (Kwon *et al.*, 2016; Afanasyev and Korabel, 2006) and numerical (Dietrich *et al.*, 1996; Ponta, 2010; and Nunalee and Basu, 2014) studies. This body of knowledge consistently indicates that the dynamics of von Kármán vortices is controlled by the Reynolds number

where *U* represents the reference velocity of the large-scale flow, *D* is the effective diameter of the obstacle, and *ν* is the molecular viscosity (see the schematic in Fig. 1). The shedding of von Kármán vortices occurs for Re ≥ 50, and as the Reynolds number increases, the process becomes progressively more irregular and turbulent (Williamson, 1996). The present study attempts to bring insight into the behavior of von Kármán vortices in the effectively inviscid limit of high Reynolds numbers (Re → ∞). To achieve such high values of Reynolds numbers, we perform two-dimensional simulations of the governing (Navier–Stokes) equations. The two-dimensional model is particularly relevant in the context of typical large-scale geophysical flows. Due to their large aspect ratios (with width-to-height ratios on the order of 100:1 to 1000:1 in the ocean), vertical displacements of fluid elements in geophysical flows are much reduced relative to their horizontal motions. As a result, such geophysical mesoscale phenomena are often conceptualized using two-dimensional models—see, for example, the reviews of geophysical vortex theories in Carton (2001) or Sokolovskiy and Verron (2014). Elements of two-dimensional dynamics are also realized in the evolution of stratified late wakes, where buoyancy forces tend to suppress the vertical displacement of fluid. Thus, a clear understanding of two-dimensional effects can yield important insights into the dynamics of wakes and vortices.

This investigation is based on a series of high-resolution direct numerical simulations (DNS) in which Re is systematically increased. Modeling high-Re systems is ultimately constrained by the requirement to resolve a wide range of vastly dissimilar spatial scales—from the scale of energy injection into the flow field to its viscous dissipation scale. At present, DNS modeling of vortex streets in most vessel generated (Re ∼ 10^{8}) and geophysical (Re ∼ 10^{10}–10^{13}) systems still remains largely inaccessible. However, due to continuous advancements in high-performance computing, the current generation of models at least starts to approach levels where it is possible to offer meaningful predictions and extrapolations aimed to represent the elusive “hard turbulence” limit.

The goal of the current analysis is the identification of the prevailing spatial, temporal, and statistical patterns of nearly inviscid von Kármán vortices. The presented DNS of two-dimensional flows impinging on a circular obstacle are performed with Reynolds numbers in the range of 5000 ≤ Re ≤ 200 000. The resulting flow fields are characterized by the ubiquitous presence of persistent coherent eddies downstream of the obstacle that are highly suggestive of von Kármán vortices. These structures frequently take the form of narrow high-vorticity cores surrounded by much broader, nearly irrotational annuli. The effective outer scale of such vortices remains finite with increasing Re. However, the high-vorticity cores exhibit singular behavior associated with the systematic decrease in their sizes (*r*_{c} → 0) and the corresponding unbounded increase in the peak vorticity values (*ζ* → ∞) in the inviscid limit (Re → ∞).

The singular distribution of vorticity in this regime raises a number of intriguing questions and opens interesting opportunities for theoretical development. Of particular significance are the proposed (Sec. III) explicit analytical models for the size and intensity of high-vorticity cores as functions of Reynolds number in the asymptotic limit of Re → ∞. The systematic analysis of DNS in the present investigation indicates that the spatial and evolutionary patterns of individual von Kármán vortices generally become well represented by the exact Lamb–Oseen solution of two-dimensional equations of motion (Oseen, 1911; Lamb, 1932),

where *r* is the distance to the vortex center, $rc(t)=rc2(0)+4\nu t$ is the effective radius of the vorticity core, *V*_{θ} is the mean azimuthal velocity, and Γ is the net far-field circulation. This finding supports the suggestion of Ponta (2010), who argued that the vortex cores in wake of a circular cylinder at lower Re exhibit vorticity profiles and decay rates consistent with the Lamb–Oseen model, yet extends the analysis to higher Re flows. While being of interest in their own right, these observations also afford the opportunity to design a low-order census of vortices in terms of two structure variables *r*_{c}(0) and Γ and explore the variation in their generation and evolution as a function of the Reynolds number. The current investigation draws attention to the dissimilar dynamics of two distinct wake regions: (i) the accretion zone, characterized by a series of vortex mergers indicative of the upscale energy cascade, and (ii) the dissipative downstream region at which Lamb–Oseen dynamics controls the evolution of von Kármán vortices.

## II. FORMULATION

In this study, we explore the dynamics of von Kármán vortices using the configuration illustrated in Fig. 1, which shows a uniform flow impinging on a circular rigid obstacle. The investigation is based on a series of two-dimensional DNS. Effects of planetary rotation and compressibility are ignored, and therefore, the governing equations are represented by the 2D momentum and continuity equations. This system is non-dimensionalized using the speed of the impinging flow (*U*_{0}) as the unit of velocity, the diameter of the circular obstacle (*D*) as the unit of length, $DU0$ as the unit of time, and *ρU*^{2} as the unit of pressure, where *ρ* is the fluid density. The non-dimensionalization reduces the governing equations to

where *p* is the non-dimensional dynamic pressure and **v** is the non-dimensional velocity field. Systems (3) and (4) can be also be conveniently expressed in terms of vorticity (*ζ*) and streamfunction (*ψ*), which reduces it to

where $\zeta =\u2202v\u2202x\u2212\u2202u\u2202y$, $(u,v)=\u2212\u2202\psi \u2202x,\u2202\psi \u2202y$, and J represents the Jacobian.

### A. Model configuration

The DNS of the governing system [(3) and (4)] is performed using OpenFOAM, version 1806—an open source finite volume numerical model initially described in Weller *et al.* (1998). To model an effectively unbounded flow field, we used the computational domain with size (*L*_{x}, *L*_{y}) = (17, 16) and (*L*_{x}, *L*_{y}) = (34, 16) greatly exceeding the diameter of the solid cylindrical object. The origin of the coordinate system is placed at the center of the object. In this frame of reference, the inflow (outflow) boundaries are located at *x* = −*L*_{in}(*x* = *L*_{x} − *L*_{in}), where *L*_{in} = 4. The calculations were performed on an unstructured grid aligned to polar coordinates near the cylinder and evolving to a Cartesian grid in the far field near the boundaries. In the x-direction, a uniform velocity (*u*_{0} = 1) is prescribed at the inflow, and an open boundary condition $(\u2202u\u2202x=0,p=0)$ is used at the outflow while restricting backflow from exiting vortices from entering back into the domain. We apply symmetry plane boundary conditions in the y-direction, which are not expected to influence the dynamics of von Kármán vortices, given the large cross-flow extent of the computational domain ($LyD\u226b1$). At the boundary of the bluff body, we apply no-slip boundary conditions. A second order central differencing scheme was utilized for spatial discretization and a backward scheme (also second order) for temporal advancement. The Pressure Implicit with Splitting of Operators (PISO) method was used for solving the pressure–velocity coupling [see Versteeg and Malalasekera (1995)]. A series of experiments (not shown) has also been performed using MITgcm (Marshall *et al.*, 1997a; 1997b) and the general consistency of both models indicates that the results are robust and not model-dependent.

The structure of governing equations (3) and (4) implies that the system dynamics is largely controlled by the Reynolds number (1), and therefore, our investigation is based on a series of simulations in which Re is systematically increased from Re = 5 · 10^{3} to Re = 2 · 10^{5} (see Table I). The increase in Reynolds number is associated with the reduction in both the scale of viscous dissipation in the turbulent wake and the width of frictional boundary layers forming around the solid object. It is of critical importance to fully resolve these scales. Therefore, the simulations were performed with fine resolutions in the vicinity of the cylinder, ranging from Δ*x* = Δ*y* = 0.0025 in the lower Re runs to Δ*x* = Δ*y* = 0.0008 for the higher Re runs. Initial selection of these resolutions was guided by the need to resolve the flow field within the boundary layer. An expansion of the grid spacing away from the cylinder wall to boundaries is limited (Δ*x*_{far} ∼ 2Δ*x*_{cyl}), and adequate resolution of all components of the flow field was validated by inspecting the spectra of enstrophy and perturbation energy and by performing selected runs with increased resolution for limited periods of time. An example is shown in Fig. 2 for our most challenging simulation with Re = 200 000. To test our results, a high-resolution simulation (Δ*x*_{fs} ∼ 0.2Δ*x*) was run for 10 units of time, and then, the results compared with the base case. KE perturbation spectra in Fig. 2 confirm that the perturbation energy spectra are grid independent and the simulations are adequately resolved. The perturbation KE values between the baseline experiment and the one with refined grid mesh differed by 6%.

Run Nos. . | Re . | No. of grid points . | Time . | N_{av}
. | N_{t}
. | λ (%)
. |
---|---|---|---|---|---|---|

1 | 5 000 | 13600 × 6400 | 75 | 19.95 | 0.665 | 94.45 |

2 | 10 000 | 13600 × 6400 | 75 | 21.56 | 0.719 | 94.69 |

3 | 25 000 | 13600 × 6400 | 75 | 37.19 | 1.240 | 56.58 |

4 | 50 000 | 8500 × 8000 | 75 | 33.59 | 2.239 | 93.29 |

5 | 100 000 | 8500 × 8000 | 75 | 50.22 | 3.348 | 76.77 |

6 | 200 000 | 17000 × 16000 | 45 | 57.98 | 3.865 | 78.10 |

Run Nos. . | Re . | No. of grid points . | Time . | N_{av}
. | N_{t}
. | λ (%)
. |
---|---|---|---|---|---|---|

1 | 5 000 | 13600 × 6400 | 75 | 19.95 | 0.665 | 94.45 |

2 | 10 000 | 13600 × 6400 | 75 | 21.56 | 0.719 | 94.69 |

3 | 25 000 | 13600 × 6400 | 75 | 37.19 | 1.240 | 56.58 |

4 | 50 000 | 8500 × 8000 | 75 | 33.59 | 2.239 | 93.29 |

5 | 100 000 | 8500 × 8000 | 75 | 50.22 | 3.348 | 76.77 |

6 | 200 000 | 17000 × 16000 | 45 | 57.98 | 3.865 | 78.10 |

For each case, the model was run for sufficient time to ensure that solutions had statistically stabilized, at which point sampling was conducted on 0.25 intervals for the final 15 time units, with the exception of runs 1 and 2, where the last 30 were collected to ensure that a sufficient number of vortices were captured for analysis. Total run times are shown in Table I. A typical simulation is shown in Fig. 3, which was performed with Re = 50 000. The instantaneous patterns of vorticity and absolute velocity are presented for t = 55. The vorticity distribution is characterized by a series of sharp isolated features, surrounded by broader annuli of perturbation velocity.

### B. Vortex classification

Our approach to analyzing the vortex field takes advantage of the general similarity of the fully developed coherent vortices in the wake of the rigid body to the corresponding Lamb–Oseen solutions (2). The vortex identification algorithm is based on the proper orthogonal decomposition of the velocity field, which was developed by Michard *et al.* (1997) and utilized by Graftieaux *et al.* (2001). Once all vortices at a given instant of time were identified, the velocity fields around each vortex were transformed into the polar coordinate system referenced to the diagnosed vortex center. The azimuthal velocities were then averaged in the polar angle, and the resulting radial patterns of mean azimuthal velocity *v*_{θ}(*r*) for each vortex were analyzed as follows: First, the patterns of *v*_{θ}(*r*) were fitted by the canonical Lamb–Oseen solution (2) with Γ and *r*_{c} as free parameters. Representative examples of the azimuthal velocity patterns and the corresponding Lamb–Oseen profiles are shown in Fig. 4.

Next, an attempt was made to ensure that the characterization of vortices in terms of Lamb–Oseen parameters was adequate by applying the quality-control criterion based on the confidence intervals (CI) of the resulting fits. If the 95% CI for the Lamb–Oseen fit was within 5% of the fitted parameter for both Γ and *r*_{c}, then the fit was considered sufficiently accurate, and those vortices were included in the subsequent analysis. Vortices that did not satisfy the CI criterion were discarded. In order to mitigate the impact of vortex interactions in the peripheral field, the comparisons between generated best fit curves and the numerical data were limited in range from the core center to 1.5 times the radius of maximum mean azimuthal velocity in the vortex.

## III. RESULTS

### A. Accretion and dissipation zones

At higher Re values, the applicability of Lamb–Oseen theory to vortex evolution only holds once a majority of strong vortices exit the immediate vicinity of the wake generating object, a region we term the vortex accretion zone. Figure 5 provides insight into the evolutionary differences between the accretion zone and the far-field regime, which we term the dissipative zone. It is clear that in the accretion zone near the cylinder, the general evolution of vortex core radii is not consistent with Lamb–Oseen theory. In this region (also see Fig. 1), newly generated vortices emerge from instabilities in the boundary layer with near singular patterns of vorticity. However, these vortices are characterized by significant variability, both in terms of their core sizes and also azimuthally. In many cases, they are not axisymmetric but elliptical, e.g., Ponta (2010), and in nearly all cases, they are subject to strong straining forces at nearly all radial ranges. Once separated from the bluff body, the majority of vortices undergo a spin-down and increase in core radius, generally accompanied by multiple interactions and mergers with other vortices. To a certain extent, the process results in a bimodal distribution of vortices entering the dissipative zone, with the resultant composition split into smaller vortices that escape interaction and larger vortices that undergo rapid accretion through interactions with other vortices or the turbulent, filamentary field. Nonetheless, the range of sizes is controlled by the relative dominance of the inertial forces generated by interaction with the bluff body in comparison with viscous dissipation.

In contrast, the dissipative zone is characterized by the predominance of axisymmetric vortices. In most cases, the vortices in this region are sufficiently separated such that mergers become much less frequent. In this region, dissipation plays the central role in the vortex core evolution, and hence, Lamb–Oseen theory becomes much more applicable. We caveat that the transition from accretion to dissipative zones is continuous in nature. Some vortices emerge into a dissipative regime almost immediately after being generated. Others undergo accretion and merger well after the majority have entered a diffusive regime. Nonetheless, we seek to determine, as accurately as possible, the approximate dividing line between the accretion zone and the dissipative zone.

To aid in this determination, we utilize the Okubo–Weiss (*OW*) parameter

where $Sn=\u2202u\u2202x\u2212\u2202v\u2202y$ is the normal strain field, $Ss=\u2202v\u2202x+\u2202u\u2202y$ is the shear strain field, and *ζ* is the vorticity. The *OW* parameter (Okubo, 1970; Weiss, 1991) and its associated components have traditionally been utilized as a diagnostic variable to identify coherent eddies in turbulent geophysical flows. Negative *OW* values are indicative of vorticity dominated flow and, therefore, generally associated with the presence of coherent eddies. In this experiment, the utility of the *OW* parameter and associated strain components lies in identifying regions where the strain field is sufficiently strong to limit the applicability of Lamb–Oseen dynamics. Figure 6 depicts a snapshot of the *OW* field for two runs (Re = 10 000 and Re = 100 000). In both panels, the presence of a large magnitude strain field (bright yellow) can be seen initially surrounding multiple vorticity dominated core regions (dark blue) that are generally not axisymmetric. The rapid reduction in strain field magnitude with X is clearly seen in both panels. This decline is evident in all other runs and occurs in conjunction with the merging and axisymmetrization of vortices.

To validate our observations in the *OW* intensity fields shown in Fig. 6, we inspected azimuthal variances of the strain (*σ*_{Sn,az}, *σ*_{Ss,az}) and vorticity (*σ*_{ζ,az}) fields surrounding each of the vortex center points identified in all our runs, defined, respectively, as

where *r* is a selected radius and $N=2\pi \Delta x$ is the number of samples taken along the radial curve. By examining each vortex in polar coordinates and plotting over their respective locations, we can identify features that could mark the transition into a more diffusive dominant dynamics regime. Fluctuations in the vortex peripheral regions outside the core could be expected to occur under the influence of external turbulent interactions. However, the presence of such fluctuations inside the core region would be an indicator that axisymmetry had not been established and that system adjustment could still be occurring. Conversely, a rapid reduction in fluctuation intensity would signal the onset of axisymmetry and a transition to the dissipative Lamb–Oseen regime. An example of this process is shown in Fig. 7, and it is clear that the radial variance rapidly decreases in the near field. This is the case for an entire range of radii both internal to the vortex core (values of *r* ≤ *v*_{max,azimuthal}) and outside the inner core region (*r* > *v*_{max,azimuthal}). All three *OW* components [Eqs. (7)–(9)] showed the substantial reductions in variance. As expected, vorticity magnitude becomes and remains much larger than the strain field.

The combination of these patterns provide strong evidence for the rapid establishment of axisymmetry and diffusive dynamics at a finite downstream location *X* = *L*_{B}, lending itself to Lamb–Oseen characterization and diagnostics. Given these observations, we set the downstream edge of the vortex accretion zone to be *L*_{B} = 2.5. Without loss of generality, we can then reference all vortex statistics to this initial point, beyond which Lamb–Oseen dynamics are the dominant process. Once beyond the accretion zone (*X* > *L*_{B}), this algorithm makes it possible to easily classify the entire set of von Kármán vortices, realized at various times and locations, in terms of two parameters: the net circulation Γ and the effective initial radius *r*_{0} ≡ *r*_{c}(0) at the exit from the vortex accretion zone.

### B. Low-order census of vortices in the dissipative zone

Since the majority of coherent vortices closely match the canonical Lamb–Oseen pattern, we now attempt to classify all eddies in terms of two parameters: (i) their effective radii at the boundary between accretion and dissipative zones (*r*_{0}) and (ii) their outer circulation, Γ. The temporal variation of the core radius in the Lamb–Oseen solution is represented, in non-dimensional units, by

where for each vortex, t = 0 corresponds to the crossing of the accretion/dissipative boundary. The presence of a temporal variable in the viscous expansion model (10) leads to a significant technical complication, typically requiring continuous tracking of individual vortices. However, this difficulty is bypassed by assuming—and subsequently verifying—that the displacement of von Kármán vortices in *x* can be largely attributed to passive advection by the large-scale flow. The passive advection model is not exact. For instance, one of the apparent limitations of this model is associated with occasional coalescence of vortices further downstream of the bluff body; such events are more prevalent in the higher Re flows. Nonetheless, the passive advection model has been validated and utilized in various systems, e.g., Lingevitch and Bernoff (1995), and in the present study, it proved to be representative of the general evolutionary trends of von Kármán vortices. Moreover, if individual vortex displacement magnitudes and re-directs due to the eddy–eddy interactions are randomly spread, the vortex motion over time will be statistically equivalent to the advection by the mean flow.

In this study, we considered two versions of the passive advection model, uniform and inhomogeneous. The uniform advection model assumes that the advective speed is close to unity for much of the flow field in the dissipative zone (*X* > *L*_{B}), and therefore, the effective age of each vortex at position (*x*_{c}, *y*_{c}) can be described by a simple relation

where *C* is a constant.

The inhomogeneous advection model makes an attempt to account for the relatively weak spatial variation in the wake velocity ($|u\xaf\u22121|\u226a1$), where $u\xaf$ is the temporal mean x-velocity of the flow. The wake defect introduces a correction to the foregoing uniform model, leading to

where $U\xaf$ is approximated by the spatial mean of $u\xaf$ taken at *y* = *y*_{c} over the x-interval [*L*_{B}*x*_{c}],

A comparison of these two methods yielded minimal differences in the statistics of the vortex field, less than 1% variation in mean radial values. This finding implies that the mean velocity defect in the wake has very limited influence on the statistics of von Kármán vortices in the dissipative zone.

In all cases, core sizes of vortices systematically increase with the downstream distance from the bluff body, and this relation is well described by the model based on the Lamb–Oseen solution as

Support for the passive advection approximation and the diagnostic validity of Lamb–Oseen theory is provided by Fig. 8, which shows the distribution of *r*_{c}(*t*) in high-Re experiments (Table I) as a function of (*x*_{c} − *L*_{B}). While Fig. 8 is an aggregate plot, one can readily see many instances where the evolution of individual vortices can be traced. While some individual vortices displace at varying speeds due to interactions with other flow components, the majority propagate with the speed of the mean current. The agreement between the observed spreading rates of individual vortices and the Lamb–Oseen prediction is underscored by presenting the best fit theoretical slope of the $rc2(x)$ relations for each Re, indicated by the green lines in Fig. 8. Another striking feature of the diagnostics in Fig. 8 is the systematic decrease in vortex radii with increasing Reynolds numbers. The physics of the phenomenon can be traced to the processes involved in the generation of von Kármán vortices, as illustrated in the schematic diagram (Fig. 1). The flow field is irrotational (*ζ* = 0) in the upstream region, and therefore, the only source of vorticity in the system is the frictional boundary layer surrounding the bluff body. The continuous generation of vorticity in this region is a consequence of a no-slip boundary condition applied at the fluid–solid interface. The increase in the Reynolds number results in systemic reduction in the width of the boundary layer and in the corresponding increase in vorticity values. During the formation of von Kármán vortices, thin filaments peel-off from the boundary layer, detach, and reorganize into circular cores (Fig. 1). It is therefore natural to expect much higher values of vorticity in—and much smaller radii of—the rotational cores of vortices at high Reynolds numbers.

An inspection of the probability distribution functions (PDF) in Fig. 9 confirm and quantify our thesis: there is a systematic decrease in both the mean (*r*_{0,mean}) and the root mean square (*r*_{0,rms}) for initial radii with increasing Re. To further explore these dependencies, we plot *r*_{0,mean} and *r*_{0,rms} as functions of Re in logarithmic coordinates (Fig. 10). Both sets of data points in Fig. 10 are distributed along the approximately straight lines, which suggests that these relations can be adequately represented by the following power laws:

The prefactors and exponents of the power laws (15) were evaluated from their best linear fit to the data expressed in logarithmic coordinates (e.g., Fig. 10), which resulted in

The proposed explicit relations (16) provide a direct link between the effective Reynolds numbers and the scaling of coherent vortices in the late wake. While it remains to be determined whether the proposed laws remain relevant for even higher Reynolds numbers, our present calculations suggest the possibility of universal relations for the core radii in the effectively invicid limit (Re → ∞). In addition to the properties of nominal initial core radii, a related set of questions arises with regard to the behavior of the core vorticity values in the high-Re regime. This challenge was addressed by introducing the nominal initial vorticity (*ζ*_{0}) defined in the analogous manner as follows:

The foregoing evidence indicates the core vorticity systematically increases with increasing Reynolds number and the following diagnostics attempt to quantify this tendency. As in the analysis of the core radii, we focus on two measures of nominal initial vorticity: its mean (*ζ*_{0,mean}) and root mean square (*ζ*_{0,rms}) values. The results reveal monotonic increase in *ζ*_{0,mean} and *ζ*_{0,rms} with increasing Re. A predictive model for the core vorticity was designed by assuming the following power laws:

As with (15), the prefactors and exponents of the power laws (18) were evaluated from their best linear fit to the data in logarithmic coordinates, which resulted in

It is interesting to note that power laws in (16) and (19) imply the approximately inverse relation between the vorticity and the effective radii of von Kármán vortices ($\zeta 0\u223cr0\u22121$). The corollary of this observation is that the increase in Re leads to a systematic reduction in both net circulations (Γ) and initial radii (*r*_{0}) that occur at a comparable rate. This compensating reduction in Γ with Re makes the increase in *ζ*_{0} with decreasing *r*_{0} less rapid than one could have expected from (16). Nonetheless, the inverse relation between the vorticity and radii of individual vortices can be clearly seen in Fig. 11. Its consistency across all values of Re tested is remarkable considering the increasingly turbulent behavior. While there is a considerable scatter in the data, the distribution is consistent with the −1 power law realized for the entire range of Re and distances (*L*_{B} < *X* < 30) in this study.

Figure 12 present higher order statistics of the core radii, which include the variance, skewness, and kurtosis. These diagnostics are consistent with the non-normal probability distribution of *r*_{0}, which is apparent from the radii scatter plots in Fig. 8 and the corresponding distribution functions in Fig. 9. In particular, the higher Reynolds number distributions tend to be bi-modal, with a dominant peak at relatively small *r*_{0} and a secondary grouping around a larger radius. The tightening of the probability distribution is expressed through the reduction in *r*_{0} variance with increasing Re in Fig. 12(a). The corresponding presence of outlier radii (the second mode) in all runs is reflected in the plots of excess kurtosis remaining positive and the slight increase in right (positive) skewness, both shown in Fig. 12(b). The positive skewness in all runs indicates the stronger influence of a small number of larger vortices on the overall statistics. Our findings in this regard are similar to other two-dimensional studies of vortex decay and accretion in turbulence (Bracco *et al.*, 2000; Burgess *et al.*, 2017).

Another Re dependent feature observed is the frequency of generated vortices. This value systematically increases with increasing Reynolds number (Table I). The process has been noted before (Bracco *et al.*, 2000). However, the increase in spawned vortices would seem to contradict established experimental evidence relating to the Strouhal number, $St=fLU$ (where *f* is the vortex shedding frequency). Multiple experiments (Roshko, 1961; Williamson, 1996) have found this value ranging from 0.19 to ∼0.3 for a range of Re from less than a 1000 to ∼200 000. However, a closer visual inspection of the wake shows evidence that two distinct scales of motion are present (Fig. 13). The smaller formation scale is intuitively expected in the higher Re cases where small coherent vortices are generated in the vortex accretion zone. We utilize a modified logarithmic vorticity, defined as

to better visualize the smaller scales and signs of vorticity in this flow, as seen in the upper panel of Fig. 13. At the same time, inspection of the modified logarithmic vorticity plot in Fig. 13 also indicates that a broader organization of these small vortices on scales commensurate with the scale of the bluff body. At higher vorticity ranges, the process can be seen in groups of like spinning vortices rotating around each other before a similar patch of counter-rotating vortices generated from the opposite side of the cylinder cuts off the generation process. The second, larger formation scale can be discerned via the perturbation stream function defined by

displayed in the lower panel of Fig. 13. From this perspective, the organization of the complete vortex field can be seen to be governed by the dimensions of the cylinder. This large-scale pattern revealed by the stream function diagnostics could represent a closer analog of laboratory observations. It is likely that the variance in these larger motions would account for the majority of energy in the system, commensurate with the scale of the larger circulation fields (Burgess and Scott, 2017). This might explain the well-established *St* values, even with the increased localized shedding frequencies seen at larger Re.

## IV. DISCUSSION AND CONCLUSIONS

This study presents a systematic analysis of spatial and statistical patterns of von Kármán vortices forming in the wake of a circular obstacle in two-dimensional high Reynolds number flows. These vortices are frequently characterized by the presence of a compact high-vorticity core surrounded by a much wider irrotational annulus. Our analysis is based on the observation that the majority of von Kármán vortices are well represented by the canonical Lamb–Oseen solution (2). This feature makes it possible to classify the entire population of coherent eddies, identified at various instances of time and at various locations, in terms of just two parameters—the net circulation Γ and the nominal initial core radius *r*_{0} upon exiting the vortex accretion zone.

The average properties of von Kármán vortices, as well as their statistical patterns, are shown to be highly sensitive to the Reynolds number (1). Of particular interest is their singular behavior in the limit of Re → ∞, associated with the systematic reduction in the effective radii of vortex cores (*r*_{0} → 0) and the unbounded increase in the core vorticity ($\zeta 0\u2192\u221e$). Both features can be physically rationalized by focusing on the formation mechanism of high-vorticity cores, which is illustrated in Fig. 1. Since the upstream flow field is inherently irrotational, the main source of vorticity in the system is the frictional boundary layer surrounding the bluff body. For high Reynolds numbers, the boundary layer becomes increasingly narrow, and the peak vorticity in it systematically increases due to the no-slip boundary condition applied at the solid/fluid interface. The boundary layer intermittently sheds high-vorticity filaments, which subsequently reorganize into circular patches that ultimately form the cores of von Kármán vortices. A corollary of the proposed scenario is that the initial vorticity values in newly formed vortices should progressively increase and the core radii should decrease, with increasing Reynolds number. This effect is confirmed and quantified by formulating explicit empirical relations for *r*_{0}(Re) on the basis of high-resolution DNS performed for a wide range of Re.

While the approximate Re^{−0.5} scaling is an expected relationship for laminar processes, it is interesting that this scaling law remains relevant even for vortices formed after a complex and turbulent merging process within the accretion zone. With increasing Re, this result becomes less intuitive as the boundary layer near the cylinder wall becomes more turbulent, and therefore, its width is governed by very different shape-factor dependent empirical scaling relationships. See Schlichting and Gersten (2017) and references therein for in depth discussion of turbulent boundary layers. That the vortex field emerging from these mergers can still be represented with the scaling laws (*r*_{0} ∝ Re^{−0.5}) expected for laminar processes, and that the size of the accretion zone remains largely independent of Re, is intriguing and deserving of continued attention. The present study can be extended in a number of promising directions. One of its obvious limitations is that the accessible Reynolds numbers may not be representative of many applications. While the resolution of the current set of DNS is close to the present-day “state of the art,” it is still too crude to represent the entire spectrum of scales realized in typical geophysical systems. At geophysical scales, the vortex cores may be turbulent, and as such, the use of kinematic viscosity might need to be replaced with an eddy viscosity. In this regard, we pin much hope on the principle of ever-increasing computational capabilities (Moore’s law). The continuous modeling advancements would make it possible one day to test our qualitative inferences and the proposed power laws for more realistic parameters. Field observations of von Kármán vortices are relatively rare and more difficult to analyze. However, the identification of analogous trends in oceanic and atmospheric flow patterns would represent the most convincing validation of model predictions. Another interesting avenue of exploration, and a major step toward realism, would be the analogous census of coherent structures in three-dimensional stratified flows. While the relative vorticity in such systems is not conserved, the diagnostics could be based on the Ertel potential vorticity (PV)—the corresponding Lagrangian invariant of inviscid inhomogeneous 3D flows [e.g., Pedlosky (1987)]. Our preliminary 3D simulations (not shown) reveal the presence of localized PV patterns originating in the boundary layer. This, in turn, implies the possibility of singular PV distribution in the high Reynolds number limit. The analysis of the structure and statistics of von Kármán vortices in the presence of pre-existing upstream turbulence is yet another challenging and largely unresolved problem to be pursued in future studies.

It should be emphasized that while the present investigation is narrowly focused on von Kármán vortices, we believe that it can bring additional insight into a much broader field of two-dimensional turbulence (Kraichnan, 1967; Kraichnan and Montgomery, 1980). The spontaneous emergence of long-lived coherent structures is a common feature of turbulent 2D flows, explored and documented in various laboratory and numerical analogs [e.g., Bracco *et al.* (2000); Boffetta and Ecke (2012)]. These structures introduce spatial intermittency that is not represented in classical theories for the cascade of energy and enstrophy, which demands full understanding of their origin and dynamics. The formation mechanisms, patterns, and stability properties of persistent vortices illustrated here on the example of high Reynolds number bluff body wakes could prove to be sufficiently generic and relevant to other types of forced 2D turbulence.

Finally, it is important to mention some potential utilitarian ramifications of our findings. For instance, one of the interesting challenges in wake theory is the determination of the age of a wake and properties of the wake-generating body based on instantaneous measurements long after its passage [e.g., Radko and Lorfeld (2018); Radko and Lewis (2019)]. An effective solution to this problem has remained largely elusive, in part due to the “loss of memory” effect—the tendency of wake turbulence to evolve in time toward a universal pattern (Meunier and Spedding, 2004; Spedding, 2014). However, our study presents an interesting example demonstrating that some wake components (i) are highly sensitive to the system parameters (the Reynolds number in our case) and (ii) evolve over time in a simple and predictable manner. Of course, at this point, it is not clear whether these features are fully realized in more realistic three-dimensional settings. If so, our work could enhance the prospects of inferring operationally significant information about moving objects using late-wake tracking systems.

## ACKNOWLEDGMENTS

This work was funded by the Office of Naval Research (Grant No. N0001420WX00354). High performance computing resources were provided through the U.S. Department of Defense High Performance Computing Modernization Program. The authors would like to thank Dr. Justin Brown for his inputs.

## DATA AVAILABILITY

Raw data were generated at the Army Research Lab Department of Defense Supercomputing Resource Center (ARL DSRC). The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

^{*}