Stratified flow past a sphere at moderate Reynolds numbers

A numerical study of stably stratified flows past spheres at Reynolds numbers $Re=200$ and $Re=300$ is reported. In these flow regimes, a neutrally stratified laminar flow induces distinctly different near-wake features. However, the flow behaviour changes significantly as the stratification increases and suppresses the scale of vertical displacements of fluid parcels. Computations for a range of Froude numbers $Fr\in [0.1,\infty]$ show that as Froude number decreases, the flow patterns for both Reynolds numbers become similar. The representative simulations of the lee-wave instability at $Fr=0.625$ and the two-dimensional vortex shedding at $Fr=0.25$ regimes are illustrated for flows past single and tandem spheres, thereby providing further insight into the dynamics of stratified flows past bluff bodies. In particular, the reported study examines the relative influence of viscosity and stratification on the dividing streamline elevation, wake structure and flow separation. The solutions of the Navier-Stokes equations in the incompressible Boussinesq limit are obtained on unstructured meshes suitable for simulations involving multiple bodies. Computations are accomplished using the finite volume, non-oscillatory forward-in-time (NFT) Multidimensional Positive Definite Transport Algorithm (MPDATA) based solver. The impact and validity of the numerical approximations, especially for the cases exhibiting strong stratification, are also discussed. Qualitative and quantitative comparisons with available laboratory experiments and prior numerical studies confirm the validity of the numerical approach.


Introduction
Stably stratified flows past a sphere display an array of intricate physical phenomena epitomising flow responses to localised three-dimensional obstacles. The benchmark emulating a spherical body moving through a stratified fluid provides a convenient framework for investigating both the flow behaviour and the numerical methods employed to simulate flows at a range of Reynolds and Froude numbers. Viscous stratified flows have important applications in diverse areas of engineering and geoscience, where examples include wind turbines and marine appliances submersed in thermoclines [9], flows of ocean currents past natural bathymetry [56] as well as atmospheric flows past mountains and islands [21,34]. Such practical interests motivate numerical simulations complementing idealised experimental studies of flows past a sphere. However, despite a large body of experimental and numerical work devoted to stratified wakes behind spheres, the numerical simulations actually resolving the viscous boundary layer past a sphere are rare due to their technical complexity; see [54,30] for recent succinct reviews.
Experimental studies, e.g. [26,15,23,3], demonstrate the importance of the relative ratio of the characteristic time scales of the viscous and buoyant forces to the time scale of inertia. These studies suggest that the structure of a flow depends largely on the fluid stratification and its interaction with viscous effects-characterised, respectively, by the Froude F r = V o /N h and Reynolds Re = V o L/ν numbers, with V o , N , h, L, and ν corresponding to free stream velocity, buoyancy frequency, obstacle height, characteristic length scale, and kinematic viscosity. As F r ∞, weakly stratified laminar flows demonstrate three-dimensional structures reminiscent of neutrally stratified flows, whereas for F r 0 they become more planar and two-dimensional vortex shedding can be observed. The latter matches satellite observations of mesoscale eddies in the lee of mountainous islands [22]. As reported by [23,3], a range of different flow patterns occurs between the two asymptotic flow regimes. The main responses such as lee waves, obstacle-scale eddies, and upwind blocking are commonly observed in experiments. However, in terms of delimiting the transitions between the patterns, experimental data has been found to be in poor agreement [13].
The pioneering numerical study of a flow past a physical sphere towed in a linearly density-stratified incompressible fluid at Re = 200 was documented by Hanazaki [14]. Hanazaki's focus was on the spatial evolution of the nearfield flow response to the sphere. Subsequent investigations spawned a range of distinct numerical strategies that were driven by diverse interests. Studies focused on temporal evolution of turbulent wakes formed behind underwater vehicles and the parameter-dependent far-field wake flow dynamics employed direct numerical simulation (DNS) [11,2,5], large eddy simulation (LES) [7], and implicit large eddy simulation (ILES) [5]. These studies played an important role in exploring turbulent wakes and parameter-dependent flow dynamics for stratified flows. However, they either completely excluded the sphere from their simulation domains while relying on suitable initial conditions, or hybridised the spatial/temporal evolution [54]-in the spirit of secondary application models where an application of interest is driven by an output of a different model [42]. Alternative models which include a sphere immersed in the moving stratified fluid and utilise high resolution grids required to accurately represent evolution of boundary layers, separation, and the dynamics of fluids close to the sphere are computationally demanding; their application is therefore still relatively uncommon.
An explicit representation of a sphere's geometry in a computational domain can adopt distinctly different approaches based on immersed boundaries or body-fitted meshes. Seminal DNS studies at Re = 3700 [31,32] solved Navier-Stokes equations in cylindrical coordinates on a staggered grid, with a sphere represented by an immersed boundary. A broad range of stratifications were considered with F r ∈ [1, ∞] [32] and F r ∈ [0.025, 1] [31], respectively. These studies provided systematic investigations of buoyancy effects on the behaviour of wakes. They confirmed that buoyancy suppresses turbulence as the stratification increases, however, for low Froude numbers below approximately F r = 0.25, a regeneration of turbulence occurs due to a new regime of vortex shedding.
Here, we adopt an alternative approach that employs unstructured meshes fitted to a sphere. It enables simulations fully resolving viscous boundary layers, important for our particular interest with quantitative drag predictions. Similar to our study, recent contributions employing body fitted meshes solved incompressible Navier Stokes equations in the Boussinesq approximation. Among them, simulations in [30] were performed for F r ∈ [0.1, ∞] at Re = 200 and F r ∈ [1, ∞] at Re = 1000. These solutions used generalised curvilinear coordinates and were obtained on a regular non-staggered grid. A structured spherical grid conforming to the sphere was also employed in [13] with a staggered arrangement for velocity and pressure. The solutions in [13] used a hybrid finite-difference discretisation in spherical coordinates, and investigated flow regimes at 0.005 < F r < 100, for 1 < Re < 500. Examples of earlier contributions utilising body fitted grids include [22], reporting simulations of flows at Re = 200 and F r ∈ [0.02, 10]. Their simulations where performed on a structured mesh generated in a cylindrical domain composed of brick and prismatic shaped elements with a mixed finite element discretisation for velocity and pressure. The ancestral work of Hanazaki [14] also employed a structured body fitted grid, in a study of flows at Re = 200 and F r ∈ [0. 25,200]. In [14] the structured grid was generated by rotating a 2D mesh created with a form of conformal mapping. Computations used cylindrical coordinates and finite difference discretisations of different orders for pressure and velocity.
The present study builds on our earlier work [51] that documented density stratified flow calculations for Re = 200 and for high Reynolds numbers in a turbulent regime appropriate for ILES. Numerical capabilities employed here extend the non-oscillatory forward-in-time finite volume (NFT-FV) approach documented in [51] through MPI parallelisation of the NFT-FV model. This extension allows for the refinement of earlier Re = 200 results, in addition to investigating stratified flow past a sphere at Re = 300. Interestingly, for weak stratifications, both cases have distinctly different wake patterns: where for high Froude numbers steady state is achieved for Re = 200, whereas a characteristic structure of the hairpin vortex shedding is observed for Re = 300. However, as stratification increases the flow behaviour of the two cases starts to show strong similarities. In addition to reducing the time to solution, the parallelisation benefits the quality of simulation results. Namely, it allows for the enlargement of the computational domain that reduces the impact of the external boundaries. Moreover, it enables the refinement of the mesh resolution in proximity to the solid boundaries and in the flow wake region. Last but not least, the simulations can be run for a longer time; so flow patterns which exhibit themselves on a longer time-scale can be better detected.
Previous numerical efforts reported in the literature have been conducted exclusively on structured meshes, whereas the current implementation permits the use of more flexible unstructured meshes and overcomes the constraint of the regular grids historically adopted in simulation of stratified flows. Structured meshes (including topologically regular grids) have merits of simplicity and computational economy for homogeneous flows in simple domains. However, geophysical flows have a large degree of heterogeneity, complex geometry of bounding domains, and multiplicity of scales. In the Earth's atmosphere and oceans internal gravity waves are both ubiquitous and intricate, as their occurrence and form depend on the relative magnitude and structure of ambient flows/currents, entropy/density stratification, and forcings. The latter in particular must account for complicated variability of shorelines and topography/bathymetry. There is also an abundance of natural multiscale phenomena relevant to weather and climate for which high variation in mesh resolution is desirable. Factors limiting the flexibility of structured meshes in construction of differential operators are usually related to rigid data structures, dictated by the topological uniformity of computational cells. Furthermore, generation of body-fitted structured meshes is restricted to moderately complex geometries. For geophysical applications, terrain following coordinates are commonly used to account for orography or bathymetry. However, this approach cannot be used for shapes involving overhanging orography, underwater caves and reservoirs, or submerged objects. The generation of structured body-fitted meshes for simple geometries such as a sphere is usually accomplished via conformal mapping. The use of multiblock techniques can add some flexibility to generate high quality structured meshes for more complex and multibody configurations, cf. [48], however, the control of the point distribution, required to capture key flow features, with conformal meshes is difficult and limited. In contrast, fully unstructured meshes adopted in this work have flexibility allowing for easy accommodation of optimal meshes and mesh adaptivity. Even more importantly, unstructured meshes allow for simulations involving arbitrary complex geometries.
For the considered flows past a sphere, unstructured meshes enable local refinement and the optimisation of mesh resolution in regions where the flow is most complex, such as the boundary layer and wake. Additionally, the capability of unstructured meshes to account for multibody configurations is demonstrated by extending our study to calculations of stratified flows past two spheres, for Re = 300, where mesh refinement in the region between the two spheres can also be used.
The reminder of the paper is organised as follows. The next section outlines the governing equations and the numerical approach employed. Section 3 documents a systematic numerical study for a range of Froude numbers 1/F r ∈ [0, 10] for flows past a single sphere at Re = 200 and 300 and two spheres in a tandem configuration at Re = 300. Remarks in section 4 conclude the paper.

Governing equations
This study is motivated by atmospheric applications, therefore the actual PDE system solved in our NFT-FV model is the Lipps-Hemler [24,25] anelastic system suitable for simulating a variety of mesoscale atmospheric flows [43,44,46]. For an ideal atmosphere, the anelastic conservation laws of mass, momentum and entropy fluctuations (1)-(3) are compactly written as where V I (=1,2,3) indicate the velocity components in the x I Cartesian coordinate directions (x, y, z), with time t, the Kronecker delta δ IJ , and potential temperature θ related to specific entropy via ds = c p d ln θ (with c p denoting specific heat at constant pressure). Density and pressure are symbolised by ρ and p, while an overline and subscript "e" mark the static reference state and a stably stratified inertial ambient state, respectively; cf. [43] for a discussion. A density normalised perturbation pressure in the momentum equation (2) is given by ϕ = (p − p e )/ρ, while the gravitational acceleration g enters the buoyancy term on the RHS of (2), where θ = θ − θ e is the potential temperature perturbation about the ambient profile θ e (x 3 ) = θ o exp(Sx 3 ); hereafter stratification S is assumed constant and "o" subscript refers to constant reference values. The deviatoric stress tensor τ in the last term on the RHS of (2) is defined by where µ = ρν is the dynamic viscosity.
The entropy equation (3) amounts to the adiabaticity statement Dθ/Dt = 0, thanks to the assumption of vanishing thermal diffusivity-the approximation consistent with laboratory studies conducted in saline solutions with large Schmidt numbers. The governing equations are deliberately written in terms of potential temperature and pressure perturbations to facilitate initial and boundary conditions as well as to enable a semi-implicit amplitude-error-free representation of the buoyant modes. To closely link the anelastic equations for low Mach number atmospheric flows with studies addressing marine applications, the governing set (1)-(3) is taken in the incompressible Boussinesq limit with ρ = ρ o , θ = θ o , and θ e (x 3 ) = θ o + Sx 3 . In this limit-valid for displacements of fluid parcels that are small compared to scale height S −1 of the stratified environment-the equations written in terms of density stratification or potential temperature stratification are mathematically equivalent.

Numerical integration
The governing equations are solved numerically using a semi-implicit NFT-FV solver, presented in comprehensive detail in [51]. Its key components consist of the unstructured and hybrid mesh based MPDATA [40,41] as the nonlinear advection operator, and a robust non-symmetric Krylov-subspace elliptic solver [43] for the Poisson equation for pressure. At the highest level of abstraction, the NFT model algorithm for the prognostic PDEs (2)-(3) can be viewed as an analogue to the trapezoidal integral of the Lagrangian ODEs underlying (2)-(3) [38]. This leads to a 4x4 implicit linear problem of the form in each node of a collocated mesh. In (5), Ψ symbolises the unknown vector of four dependent variables-namely, three components of V and θ -at the future t + δt time level, Ψ marks the explicit part of the solution including advection of the past Ψ values at t combined with their corresponding RHS forcings, and R(Ψ) symbolises the future RHS forcings comprising the pressure gradient, buoyancy and the convective derivative of the ambient θ e . 1 The linear problem (5) is inverted algebraically in each node, leading to the closed form discrete expression where V is a modified explicit part of the solution V, C is a 3x3 matrix of known coefficients, and ϕ is the unknown solution for the pressure variable. Substituting (6) into the discrete form of (1) forms the discrete Poisson problem for ϕ . In the spirit of the exact projection, it ensures that the updated velocity V can satisfy the discrete mass continuity to a round-off (as opposed to a truncation) error. This is particularly important for simulation of stratified flows, as it assures the compatibility of conservative MPDATA advection with the mass continuity and eliminates fictitious buoyancy forces in regions of locally undisturbed isentropes. The updated velocity is subject to Dirichlet boundary conditions, which implies consistent Neumann boundary conditions on the resulting pressure field. The reader interested in further details of semi-implicit NFT integrators and the associated elliptic solvers is referred to [51,49,39] and references therein.
Spatial discretisation adopts a median-dual finite-volume approach with edgebased connectivity outlined in [50]. The primary meshes employed in this work, illustrated in Fig. 14, combine prismatic elements near solid boundaries and tetrahedral elements elsewhere. The median-dual cell faces are constructed from the primary meshes, where the midpoint of each edge is connected both to the barycentres of the surrounding polyhedra and to the centres of the polygonal faces sharing the edge. Noteworthy is the synergism of the collocated arrangement of all dependent flow variables and the proven ILES property of MPDATA. On one hand, the data collocations enables algebraic inversion of the linear problem (5) at the heart of the semi-implicit NFT integrators, which generally can be quite complex [33,47]. On the other hand, collocated meshes admit non-trivial null spaces of the discrete differential operators, including in particular the generalised Laplacian that forms the viscous terms. The latter can manifest as a noise at fine scales near the limit of mesh resolution, even for low Reynolds number flows that are otherwise well resolved. The solution regularisation at the limit of mesh resolution characteristic of ILES remedies this problem. Importantly, as the MPDATA ILES regularisation is adaptive (as opposed to additive) to explicit dissipation, it has a negligible impact on the resolved scales [27,28,6]. In effect, our NFT-FV schemes benefit from the collocated data arrangement while being inherently stable even with lownumerical-diffusion [46,51,20].
While NFT codes operating on Cartesian meshes and exploiting the advantages of the distributed memory paradigm are well established, the successful implementations of NFT-FV codes for unstructured meshes on multi-node architectures are more recent. To date these efforts have focused on global models for the Earth's atmosphere with meshes that are quasi-uniform in the horizontal and prismatic in the vertical [45,20]. In this paper, we complement our earlier work with new NFT-FV model developments that provide capabilities for parallel simulations on fully unstructured, irregular meshes. We utilise the multilevel graph partitioning scheme supplied by the MeTis package [19] to assign partitions of the computational spatial domain. This library effectively balances the sub-domain workloads for irregular meshes, while minimising the number of edges crossing partition boundaries.
The NFT-FV solver utilises a localised stencil on collocated data, and an auxiliary bespoke OpenMP parallelised code configures partition halos, with a halo depth of one. The second order operators found in (1-3) and the Poisson equation therefore require an internal halo exchange after applying first order operators, but benefit from a reduced memory footprint of internal variables. Halos operate on a partition-to-partition exchange basis and halo points are chunked in the data-structure. Most MPI communications exchange information only between immediate neighbours, with some limited global communications as required by the Krylov-subspace solver-all of which utilise the standard MPI library.

Problem formulation
Simulation results for neutrally and stably stratified flows past a single sphere are discussed first, due to its extensive existing research history and the availability of pertinent data. Following [18,53], In the quantitative analysis, the drag C d , side force C s and lift C l coefficients, together with the Strouhal number St, are used to evaluate the results. For a flow impinging on a sphere, the force coefficients are defined as where A = π(D/2) 2 is the reference area, and ζ indicates indices d, s or l. The forces F κ -with κ = 1, 2, 3 corresponding to drag, side and lift forces acting on the sphere-are computed according to where for any index ι, n ι is the ι th component of the unit vector normal to the surface of the face dS of the finite-volume element and repeating indices imply the summation.
The Strouhal number is defined as St = f D/V o , where f is the frequency of the vortex shedding. Since V o and D are equal to one, the Strouhal number can be directly derived from the frequency spectrum of the relevant coefficient. We also note that for the sphere the Froude and Reynolds numbers become

Flow past a sphere at Re = 200
The steady state behaviour of a neutrally stratified flow past a sphere at the Reynolds numbers between 20 and 210 is characterised by the axisymmetric flow in the wake region. It is well studied by means of direct numerical simulations [1] and verifiable by experimental observations [52]. Computations with the NFT-FV scheme for Re = 200 and F r = ∞ are documented in [51]. They show that the resulting drag coefficient C d = 0.774, the length L r = 1.42D of the steady axisymmetric recirculation bubble behind the sphere, and the boundary layer separation angle φ s = 116.6 • are in agreement with values obtained from the experiment [29] and other numerical studies, e.g. [10,53,1]; see section 4.2 in [51] for further comparisons and references.
In line with [23], the NFT-FV computations indicate that for Froude numbers higher then F r ∼ 10, the flow resembles the neutrally stratified case; and, as the Froude number decreases the density stratification starts to significantly impact the flow behaviour [26,23,3]. For Re = 200, three dominant flow patterns have been identified in [23] throughout the Froude number regimes. However, the transitions between these flow patterns are not clearly defined and exhibit intermediate flow features leading to alternative classifications. For example, [13] distinguishes further different flow characteristics for low and moderate Reynolds numbers. Our simulations for a range of F r were carried out using the NFT-FV scheme until flow patterns have been established. The non-dimensional simulated time of T = t V o /D = 30 was used in [51]. The parallel implementation of the scheme enabled fivefold T = 150 to tenfold T = 300 extensions of the simulated time, used in this work. Numerical simulations of stably stratified flows past a sphere for moderate Reynolds numbers, especially for Re = 200, have been reported before (e.g. [30,13,31,32,22,14]) but, to the best of our knowledge, this study together with our earlier investigation in [51] are the first to document such simulations  This presentation is known as the Q-method for vortex identification [17]. The left column in Fig. 1 illustrates the 3D structures of the wake flow obtained for Re = 200 in horizontal and vertical planes at the reciprocal Froude numbers 1/F r = {0.1, 1.6, 4}, 2 representative of the three flow patterns identified in [23]. Figure 2 illustrates the corresponding vorticity patterns. For the consecutively increasing values of 1/F r the three patterns represent, respectively, a non-axisymmetric attached vortex, lee-wave instability, and two-dimensional vortex shedding. They are in agreement with the experimentally-based schematics provided by Figs. 2 and 3 in [23]. Due to the buoyancy-induced gravity waves radiating in the lee, even for weak stratification 1/F r = 0.1, the flow past a sphere loses its axisymmetric character even though viscous effects are still strong. Planar symmetry of the wake can be observed separately in the horizontal and vertical central planes in both the 3D images and vorticity plots, representing a steady recirculation bubble formed behind the sphere. Noteworthy, the symmetry of the solution is not affected by the use of irregular unstructured meshes. As the Froude number decreases, the separation point on the vertical plane moves towards the rear of the sphere [23] and the breaking of axial symmetry becomes more pronounced.
In particular, at 1/F r = 0.4 the vertical deformation results in a characteristic doughnut-shaped attached vortex, illustrated in the streamwise velocity plot at the x = 0.51 plane in the left image of Fig. 3. At 1/F r = 0.8 the increased stratification causes further deformation, such that the doughnut-shaped vortex becomes doubly connected, leaving two attached vortices as shown in the right image of Fig. 3. These results confirm the three-dimensional interpretative experimentally-based sketches in Figs. 20a,b in [23]. Similarly to the experimental observation (Fig. 19) and corresponding sketch (Fig. 20c) in [23], at 1/F r = 0.8 the NFT-FV computations indicate the formation of the bulge-like particle streak pattern (not shown).
The bulge-like particle streak pattern also occurs at 1/F r = 1, when the leewave amplitude approaches its maximum. Figure 4 presents the streamwise velocity plots and the associated streamlines at the planes z = 0 and y = 0. The bulge-like particle streak pattern can be clearly seen in the lower panel and the streamlines show that the flow remains attached to the sphere, while the separation angle moves towards the aft stagnation point. The patterns represent a transition between the non-axisymmetric and the lee waves instability regimes. They are in good agreement with those observed experimentally; cf.  horizontal and vertical planes are now distinctly different. The special case of F r = 1 is further discussed in [51] in the context of the dividing streamline [16] and the linear theories. Consistently with results reported in [3], at 1/F r ≈ 1.2, the wake collapses resulting in a "bow-tie"-shaped separation line, with a bulge at a distance of 1.5D from the rear stagnation point.
The propagation of the lee waves becomes profound between 1/F r = 1.6 and 1/F r = 2. The 3D plots in Fig. 1, vorticity patterns in Fig. 2, and the streamwise velocity components plots with associated streamlines in Fig. 5 were obtained for 1/F r = 1.6. They are representative of the lee-wave instability regime and show similar features to experimental patterns in Fig. 14 (Re = 500, 1/F r = 2.63) in [23], Figs. 6ii,iii (Re = 329, 1/F r = 1.66) in [3], and Fig. 1a (Re = 144, 1/F r = 2.27) in [15]. Here, the concept of the dividing streamline becomes relevant to the resulting flow pattern. The fluid parcels above the height of the dividing streamline have sufficient kinetic energy to flow over the top of the sphere and, due to the presence of lee waves, have a down-slope velocity greater than the mean stream; cf. [3]. The attached recirculation becomes elongated, Fig. 5, while the lee waves in the vertical centre-plane are clearly visible, although dissipated further downstream. The characteristic overturning motion of fluid parcels, due to the combined effects of buoyancy and shear [3,23], are represented by tiny recirculation regions that can be seen in the lower image of Fig. 5. They give the name to the "lee-wave instability regime" [23]. 3 Further increasing the stratification results in a predominantly two-dimensional 3 Note that the streamline plots are constructed from data stored on irregular meshes, which affects their symmetry.  that the time from the model initialisation to the onset of flow periodicity depends on the stratification. For 1/F r = 2.25 the periodicity of the flow is hardly visible, as the amplitude in the C d plot is almost one order smaller than for 1/F r = 2.5. However, examination of the corresponding streamlines and streamwise velocity at the central horizontal plane, Fig. 7 top, reveals that the periodic two-dimensional vortex shedding already takes place at 1/F r = 2.25, while for 1/F r = 2, Fig. 7 bottom, the flow is still steady but its symmetry in the horizontal plane has been already broken. The NFT-FV results confirm the experimental [23] and numerical [12] findings reporting that the flow remains stationary until 1/F r ≈ 2.5, while computations in [22] indicate that the flow remains stationary for 1/F r < 2. At 1/F r = 4, illustrated in Figs. 1 and 2, the flow follows a distinct twodimensional vortex shedding pattern in the horizontal central plane. The fluid parcels in the central part of the flow, between the dividing streamlines, do not have sufficient kinetic energy to travel over and under the sphere. Strong stratification forces them to flow around, resulting in a motion similar to a flow around a cylinder. In contrast, the parcels travelling above and below the dividing streamlines appreciate the obstacle as a hill of a moderate height which gives rise to the vertical pattern of linear gravity waves characteristic of flow past isolated mountains. Both horizontal and vertical patterns are clearly visible in Figs. 2 and 1. For 1/F r = 4 the vorticity plots illustrate a significant decrease of the gravity waves' amplitudes and lengths by comparison with results obtained for 1/F r = 1.6. Similar patterns to those shown in Fig. 1 for this flow regime were obtained both computationally (Figs. 5c-d in [13]) and experimentally (Fig. 10e in [23], and Fig. 5a in [3] with Re = 1961). The computations for 1/F r = 4 in [51] used a serial code and, following [14], adopted the non-dimensional simulation time T = 30. The current results show that for stratifications with 1/F r > 2.5, a longer simulation is required to trigger a transition from a steady to an unsteady pattern. The four panels in Fig. 8 illustrate the evolution of the flow pattern for 1/F r = 4 that gradually transitions from the symmetric steady state reported in [14] to a quasi-twodimensional vortex shedding. In experimental and computational studies, the transition to unsteady flow is typically triggered by arbitrary perturbations. Indeed, the addition of a white noise with amplitude 0.2V o introduced into the potential flow initialisation substantially accelerates the transition, and the onset of shedding is sensitive to the magnitude of the perturbations, Fig. 9. However, after long simulation times, all computed flows developed the same vortex shedding pattern as that shown for T = 150 in Fig. 8. Furthermore, a deliberate test replacing the second-order-accurate MPDATA in the NFT-FV solver with an over-diffusive first-order upwind scheme, shows that excessive numerical diffusion can entirely prevent the transition to an unsteady flow regime. Physically this is not surprising as the onset of shedding is related to the viscous boundary layer separation whereas the established flow depends robustly on the Reynolds and Froude numbers. Numerically, however, these results attest to the solution quality and robustness on the unstructured nonsymmetric meshes.
Following [14], the influence of stratification on the drag is quantified in terms of the departure of the drag coefficient ∆C d = C d (Re, 1/F r) − C d (Re, 0) from its value for the neutrally stratified flow. Figure 10 compares our calculations for 1/F r ∈ [0, 10] with the results obtained experimentally in [26], and computationally in [14] and [30]. For unsteady flows the drag coefficient in our calculations was averaged over the time interval T ∈ [100, 150]. Since the experimental data in [26] were obtained at varying Reynolds numbers, Fig. 10 shows  [14] show a significant departure from other data. For 1/F r ∈ [1,2], the remaining computations are in excellent agreement. In the range 1/F r ∈ [2,3], the NFT-FV results match the experimental values closely and show an increasing slope. At 1/F r = 3, the results from experiment, computation of [30] and the current simulation all show a sharp decrease in slope. This change of regime is consistent with the linear theory of [35,36] and the numerical/laboratory results in [37,55]. Between 1/F r ∈ [3.5, 5], experimental data is only available for different Reynolds numbers, however, both [30] and the current NFT-FV computations show similar decreasing slope tendencies as the experiments in [26]. For 1/F r ∈ [5, 10], both computations give comparable results but their departure from the experiments is significant. The authors in [30] attribute this to an insufficient mesh resolution for a strong stratification which could also be a factor affecting the NFT-FV results. We shall return to this issue in section 4. Further quantitative analysis compares the Strouhal numbers evaluated from the periodic evolution of the drag and side force coefficients at 1/F r = 3, 5, 10. Table 2 displays the NFT-FV results together with the corresponding outcome obtained in [22]. The results are in good agreement overall, with the magnitude of St evaluated from the drag coefficient being typically twice as large as that reached by St computed from the side force coefficient. The latter is consistent with an interpretation of the 2D vortex shedding from the sphere being akin to the vortex shedding in a flow past a cylinder [22]. The Strouhal numbers evaluated from the lift coefficients (not shown) obtained by the NFT-FV scheme are close to those obtained from the side force coefficients. Table 2 Strouhal numbers based on C d and C s coefficients (cf. [22]).

Flow past a sphere at Re = 300
As observed experimentally (Figs. 39a,b in [18]) and numerically (e.g. Fig. 2c in [13]), the neutrally stratified flow past a sphere at Re = 300 results in periodic vortex shedding. The drag coefficient obtained from our NFT-FV simulation for Re = 300 and 1/F r = 0 is evaluated by averaging the drag force over the non-dimensional time interval between T = 250 and T = 300. Its mean value C d = 0.664 and oscillation amplitude C d = 0.002 are consistent with the results in [51] and the references therein. The Strouhal number derived from the frequency spectra of instantaneous C d is St = 0.138. It compares well with St = 0.136 and St = 0.137 reported in [18,4] and [53,57], respectively.
For weak stratification, 1/F r = 0.1, the character of the neutrally stratified flow vortex pattern is maintained; see Fig. 1 illustrating 3D features in the near wake and contrasting them with a distinctly different 1/F r = 0.1 flow at Re = 200. The units of the periodic wake structure consist of hairpin vortices with planar symmetry preserved in time and set by the incipient shedding. Three vortex-shedding structures can fit in the adopted computational domain. As the stratification increases the vortex shedding is suppressed and the flow becomes close to steady-state, Fig. 11, apart from tiny recirculation regions that form due to the combined effects of buoyancy and shear, just like for Re = 200. Figure 1 shows that for lower Froude numbers, 1/F r = 1.6 and 4, stratification dominates viscous effects and the flow patterns for Re = 200 and Re = 300 become close. The results for Froude numbers where the lee-wave instability regime occurs are also similar for both Reynolds numbers, with the lee-waves pattern being already well-developed at 1/F r = 1.6. However, for Re = 300 the change from the lee-waves instability into two-dimensional vortex shedding structure occurs around 1/F r = 1.8; i.e. earlier than for Re = 200. Furthermore, a closer inspection of Figs. 11 and 5 shows that at the same 1/F r = 1.6, the length of the steady recirculation formed behind the sphere in the horizontal plane increased by comparison with that measured for the Re = 200 case. Figure 12 compares the NFT-FV results with the experiment [26] for ∆C d as a function of 1/F r. As for the case of Re = 200, the drag coefficient in the NFT-FV calculations was averaged over the time interval T ∈ [100, 150]. Only the data for the Reynolds numbers close to 300, listed in Table 3, are shown in the graph. The results are well matched for 1/F r = 1.8 and display appropriate slope changes up to 1/F r ∼ 3. The results for 1/F r = 3 and 1/F r = 4 are close to the experimental values. Notably, the experimental result for 1/F r = 3.6 does not follow the slope defined by experiments at 1/F r = 3 and 1/F r = 4. This and the similar discrepancy for 1/F r = 4.8 are possible manifestations of the magnitude of experimental uncertainty. However when 1/F r > 5, the decline of the slope for the simulation results is system- atically shallower than in the experiment. A direct comparison of ∆C d for Re = 200 and Re = 300 (not shown), revealed differences for 1/F r ∈ [1.2, 1.8] which could be associated with an earlier change to the two-dimensional vortex shedding regime at Re = 300. Differences are also present for strongly stratified flows at 1/F r ∈ [5,10]. These are expected since Fig. 1 indicates that the amplitude of lee waves for Re = 300 in the vicinity of the sphere is larger than for Re = 200. Furthermore, a larger physical viscosity elicits a delayed development of a periodic vortex shedding. For example, for 1/F r = 4 the transition to the fully developed periodic flow takes place at T ≈ 100 for  The new parallel implementation of the NFT-FV code operating on fully unstructured meshes was extensively tested for the flow past a sphere at Re = 300. The computational mesh employed in the presented calculations (described in §2.3) is highly irregular, with the node spacing varying from 0.005D in proximity to the sphere surface to 1.00D at the outer boundaries of the model domain. Figure 13 highlights the code's speed-up, tested for Re = 300 and F r = 1. The speed-up ratio T 1 /T N , compares the computational time of the sequential and the N -processes runs needed to complete the simulation. The NFT-FV code shows good scaling until 160 cores, where the parallel efficiency is 0.88 on Hydra's Intel Xenon CPU E5-2680v2 HPC. 4 Beyond that, the parallel performance becomes affected by the limitations of the mesh size and partitioning variables of the MeTis library.

Flow past two spheres
To illustrate flexibility of unstructured meshes, the NFT-FV solver was applied to simulate flows past two spheres at Re = 300. Two identical, D = 1, spheres were placed in a tandem arrangement aligned with a free stream flow. The neutrally stratified flow past the two spheres was simulated for T = 400. Figure 15 shows the streamline plot in the vertical central plane, y = 0. The streamlines reveal non-symmetric recirculations in the lee of both spheres. The shear layer generated by the windward sphere becomes reattached to the second sphere and interacts with the second sphere's boundary layer. This affects the wake behind the second sphere and decreases the size of the largest recirculation by a factor of 0.58 in comparison with the analogous single sphere case. The upper row of Fig. 16 displays the 3D flow pattern obtained for Re = 300 and 1/F r = 0. The Q-method visualisation shows that the hairpin vortices characteristic of the Reynolds number 300 flow past a single sphere (cf. [51]) are also developing in the present case. Furthermore, the right panel clearly shows how the wake generated by the first sphere overlaps the second sphere. Similarly, as with the single sphere in [51], the symmetry plane for the tandem configuration for the neutrally stratified flow is formed in the diagonal plane, rotated 45 • about the X-axis. Notably, even a weak stratification imposes the vertical position of the symmetry plane; cf. Re = 300, 1/F r = 0.1 case in Fig. 1.
Multiple vortex-shedding frequencies for each sphere were reported in [57]. The left panel of Fig. 17 shows the history of the lift coefficient obtained for both spheres. Frequencies from the NFT-FV simulation match closely those of Fig. 16 in [57]. The phase is shifted due to different initial conditions and the amplitude differs slightly due to the rotation of the symmetry plane in our calculations. and C d = 0.439 for the first and second sphere respectively.
Two dominant vortex-shedding frequencies f 1 and f 2 are found for each sphere. Table 4 provides the Strouhal numbers where f 1 and f 2 are derived from the frequency spectra of both drag and lift coefficients (not shown) in the time interval T ∈ [300, 400]. They are in a good agreement with those reported in [57]. Overall, these results reflect the influence of the first sphere's leeward flow reattachment in generating the vortices behind the second sphere. Table 4 Strouhal numbers for the Re = 300, 1/F r = 0 flow past two spheres.
First sphere Second sphere The flow past two spheres is dramatically influenced by the strength of stratification, as was previously observed for the single sphere. For 1/F r = 1.6, the symmetric steady state flow is established demonstrating the features characteristic to the lee-wave instability regime behind both spheres, Fig. 18. Consistent with the single sphere case, the vertical scales are inhibited by the stratification. Additionally, in comparison to Fig. 11 the length of steady recirculation formed behind the first sphere in the horizontal plane is approximately one diameter longer and stretches in the entire gap between the spheres. The recirculation after the second sphere is smaller and closer in length to the single sphere recirculation length. The drag coefficients computed for the first and Boundary layer separates behind both spheres, as shown in the lower panel of Fig. 18 for a vertical cross-section at y = 0. The separation angles, denoted by φ s , are measured clockwise from the upstream centerline on the vertical plane, where φ s = 164 • and φ s = 145 • for the first and second sphere respectively. In comparison, the separation angle for the single sphere at the same Reynolds and Froude numbers (Fig. 11) is φ s = 166 • , demonstrating that the separation angles behind the first sphere and the single sphere are similar. The first sphere blocks the incoming flow slowing down its motion between the spheres, and the flow separation limits the width of the wake induced by the first sphere. These two factors lead to a substantially reduced separation angle for the second sphere.
Similarly, for the strongly stratified flow at 1/F r = 4, the separation angle for the first sphere, φ s = 131 • , is higher than for the second sphere, φ s = 121 • , while the value for first sphere is lower than φ s = 138 • observed for the corresponding flow past a single sphere. The middle and bottom panels of Fig. 16 show 3D plots for flows with 1/F r = 1.6 and 4. The two-dimensional vortex-shedding regime at 1/F r = 4 also preserves the main features of the flow past a single sphere. As was observed for 1/F r = 0 and 1.6, the mean drag coefficients for 1/F r = 4 evaluated for the first sphere, C d = 1.477, is close to that obtained for the single sphere case. However, the mean drag coefficient computed for the second sphere, C d = 0.960, is significantly lower, reflecting the first sphere's blocking of the flow seen by the second sphere. The corresponding Strouhal numbers are St 1 (C d ) = 0.020 and St 2 (C d ) = 0.327 for the first sphere, and St 1 (C d ) = 0.020 and St 2 (C d ) = 0.337 for the second sphere. The Strouhal number for the second sphere is similar to that for the single sphere, indicating similarities in the vortex shedding in the far wake, visible also in Fig. 16. In contrast, the Strouhal number for the first sphere is strongly affected by the presence of the obstacle in its immediate wake.

Concluding remarks
A study exploring behaviour of stratified flows past spheres have been conducted using a newly developed parallel option of the nonoscillatory forwardin-time scheme operating on flexible unstructured meshes and employing finite volume spatial discretisation. The study repeated well documented benchmarks for neutrally-stratified flows at Re = 200 and 300, allowing for validation of the proposed numerical approach. Additionally, for the stably stratified flows at Re = 200, earlier numerical and experimental studies provided further means for qualitative and quantitative validation of the NFT-FV model.
Our study of stably stratified flows at Re = 300 provides further insights into stratified flows at moderate Reynolds numbers. Detailed documentation of numerical simulations and their analysis was supplemented with comparisons with experimental data when available. The presented results confirmed the general experimental characterisation for stratified flows past a single sphere, introduced in [23]. In particular, they showed that the overall flow patterns are dominated by stratification as the Froude number decreases. Although entirely different at high Froude numbers-where the flow at Re = 200 reaches the steady state and at Re = 300 results in a periodic hairpin vortex sheddingat 1/F r = 1.6 both flows become planar and are in the lee-wave instability regime. In this regime, both flows reach steady state except for tiny recirculations in the lee indicating the lee-waves breaking. Differences between the two flows are in the details, such as the lengths of the recirculations formed behind a sphere that can be observed in the central horizontal plane. This similarity between the flow patterns continues at lower Froude numbers where the flows enter a two-dimensional vortex shedding regime in the horizontal, while remaining planar in the vertical, e.g. at 1/F r = 4. Details of the stably stratified flow features, such as amplitudes and length of gravity waves, slightly differ between the two Reynolds numbers. For Re = 200, we also document the triggering of a transition to the 2D vortex shedding regime at 1/F r = 2.25. The exact transition between the regimes is dictated by viscosity and occurs at lower Froude number for lower Reynolds numbers.
In addition, the single sphere study was extended to akin simulations of flows past two spheres, facilitated in the framework of unstructured meshes. The results for the neutrally stratified flow at Re = 300 are in good agreement with those reported in the literature [57] and further confirm the validity of the NFT-FV model. The stratified flow past two spheres is, to our knowledge, documented for the first time and investigated for the chosen tandem configuration at Re = 300 with 1/F r = 1.6 and 1/F r = 4. The resulting flow patterns generally follow that recorded for the single sphere, however, a detailed analysis reveals particularities of the spheres' interactions that can be quite complex. These interactions can be associated with intricacies of atmospheric/oceanic flows past mountain ranges, islands and complex bathymetry [21,34,56], since they depend on the relative placement of the two spheres [57]. The tandem sphere simulations illustrate the potential of the NFT-FV solver operating on flexible meshes to provide a further insight into the class of stratified flows pertinent to engineering and geophysics, suggesting interactions between multiple spheres. Our preliminary simulations indicate that, similarly as for neutrally stratified flows, the flow patterns for multiple spheres are sensitive to the distance between the spheres and their relative location with respect to the flow direction. Identification of the flow regimes for configurations of multiple spheres will be studied in our future work.
Regardless of the overall satisfactory outcome of the study, it remains puzzling why numerical results for flows past a single sphere depart from experimental measurements for strong stratifications with 1/F r > 5, Figs. 10 and 12. Furthermore, the results using two distinct numerical models-ours and that of [30], available for Re = 200-give drag coefficients that agree well with each other, Fig. 10. The numerical results using yet another model [12], available for 1/F r = 5 and Re = 100, also appear to overestimate the drag coefficient. The design of numerical experiments is essentially the same for our and [30] calculations, where reduction of the Froude number is achieved by increasing the ambient stratification S (defined in §2.1). However in [12] the stratification is fixed, and the Froude number reduction is achieved by increasing the sphere's radius. The authors in [30] attribute the drag discrepancy to the loss of the resolution, which is plausible because as F r 0 the amplitude of the fluid parcel vertical displacements, η, scales as η = 0.5DF r, based on the dividing streamline arguments [51]. However, based on the same arguments, increasing the sphere radius effectively fixes the magnitude of the displacements and thus their effective resolution. Consequently, the result of [12] is at odds with the insufficient resolution hypothesis. Moreover, as the only common element in all three numerical designs is the Boussinesq approximation adopted in the governing equations, one might speculate invalidity of the Boussinesq model in the limit of small Froude numbers. Indeed, in our and [30] experiments the ratio of the parcel displacements and the density (or potential temperature) height scale S −1 increases as F r 0. Yet this ratio ηS is fixed in [12], which is again at odds with such speculations. Running out of ideas on how to blame the numerical results, we turn to experimental data. The design of the laboratory experiments in [26] was special, in that the size of the sphere and the ambient stratification were fixed with the buoyancy frequency N ∝ Re/F r (cf. §5.2 in [51] for discussion), effectively associating the Froude number reduction with the threading (viz. free stream) velocity. In such a design Re 0 as F r 0, with C d and ∆C d vanishing in the limit. On one hand, this shows an elementary incompatibility of these experimental data with the numerical experiments discussed, or vice versa. On the other hand it adds to the intricacy of the F r 0 limit, potentially rejuvenating interests in the asymptotic [8] theories.