Tim Colonius

The Nonlinear One-Way Navier-Stokes (NOWNS) equations have been demonstrated to accurately evolve the early stages of boundary-layer transition, relative to direct numerical simulation (DNS), for low-speed flows. While NOWNS entails a higher computational cost than the nonlinear parabolized stability equations (NPSE), it is a more robust one-way approximation that accommodates forcing, converges for larger disturbance amplitudes, and accurately evolves non-modal and multi-modal disturbances. NOWNS is based on the OWNS projection (OWNS-P) approach, while the OWNS recursive (OWNS-R) approach is a more recent development that entails a reduced computational cost compared to OWNS-P. OWNS-P solves the projected (one-way) equations, while OWNS-R solves the two-way equations and projects the resulting solution. To reduce the computational cost of NOWNS and improving its scalability, we adapt the NOWNS method to employ a recursive approach similar to OWNS-R. We reformulate OWNS-R as a project-then-solve approach which improves its performance in linear cases and enables robust solution in the nonlinear case. NOWNS-R is validated against DNS for low-speed oblique-wave breakdown and K-type transition, where it accurately captures disturbance growth and the onset of transition. Finally, the method is demonstrated for oblique breakdown of Mack's second mode in a Mach 6 flat-plate boundary layer, illustrating its ability to model nonlinear high-speed transition mechanisms prior to breakdown. Consistent with DNS studies in the literature, we observe the formation of streaks before transition onset.

A shock-kinematic boundary condition (SKBC) formulation is developed for receptivity calculations of hypersonic boundary layers with attached shocks. The boundary condition is obtained from the linearized Rankine–Hugoniot equations by coupling the shock displacement to incoming and outgoing shock-normal characteristic amplitudes, and is incorporated into the linearized Navier–Stokes equations for analysis of Mach-6 flow over a 7 degree half-angle sharp tipped wedge. The formulation is used to study the response of the wedge boundary layer to two-dimensional slow acoustic, fast acoustic, vortical, and entropic disturbances incident on the shock. Parametric sweeps over incidence angle and forcing frequency show that the SKBC reproduces the expected branch-dependent shock transmission. Increasing frequency shifts the dominant amplified wall-pressure response upstream for all forcing classes. Vortical forcing generates a wall-pressure response indirectly, through shock processing and subsequent coupling to compressible post-shock disturbances, while entropic forcing primarily produces a convected thermodynamic disturbance with a weaker pressure signature except when the forcing frequency places the response in an unstable region. These results demonstrate that the SKBC provides a useful reduced framework for isolating shock-mediated receptivity mechanisms without resolving the full upstream disturbance field by direct numerical simulation (DNS).

Global resolvent analysis provides an efficient representation of coherent structures in turbulent jets, but predictive noise modeling remains limited by the need to specify forcing statistics. We explore a data-driven approach in which local RANS-derived quantities are mapped to resolvent forcing correlations within a localized input–output framework that connects near-field regions to a global far-field observer arc. The forcing correlation matrix is constructed from frequency-dependent latent vectors produced by a neural network, ensuring a positive semi-definite model by design. Trained on isothermal round jets across a range of Mach numbers and evaluated on both interpolated and extrapolated conditions, the approach captures dominant self-correlated structure and reproduces key spectral trends, but fails to recover cross-spectral phase relationships between spatially separated regions. These results indicate that while localized, RANS-driven parameterizations of resolvent forcing are a plausible path forward, additional physical constraints are needed to obtain predictive models.

Many problems of interest in engineering, medicine, and the fundamental sciences rely on high-fidelity flow simulation, making performant computational fluid dynamics solvers a mainstay of the open-source software community. Previous work MFC 3.0 was made a published, documented, and open-source solver via Bryngelson et al. Comp. Phys. Comm. (2021) with numerous physical features, numerical methods, and scalable infrastructure. MFC 5.0 is a significant update to MFC 3.0, featuring a broad set of well-established and novel physical models and numerical methods, as well as the introduction of GPU and APU (or superchip) acceleration. We exhibit state-of-the-art performance and ideal scaling on the first two exascale supercomputers, OLCF Frontier and LLNL El Capitan. Combined with MFC’s single-accelerator performance, MFC achieves exascale computation in practice, and achieved the largest-to-date public CFD simulation at 200 trillion grid points as a 2025 ACM Gordon Bell Prize finalist. New physical features include the immersed boundary method, N-fluid phase change, Euler–Euler and Euler–Lagrange sub-grid bubble models, fluid-structure interaction, hypo- and hyper-elastic materials, chemically reacting flow, two-material surface tension, magnetohydrodynamics (MHD), and more. Numerical techniques now represent the current state-of-the-art, including general relaxation characteristic boundary conditions, WENO variants, Strang splitting for stiff sub-grid flow features, and low Mach number treatments. Weak scaling to tens of thousands of GPUs on OLCF Summit and Frontier and LLNL El Capitan achieves efficiencies within 5% of ideal to over 90% of their respective system sizes. Strong scaling results for a 16-times increase in device count show parallel efficiencies over 90% on OLCF Frontier. MFC’s software stack has undergone further improvements, including continuous integration, which ensures code resilience and correctness through over 300 regression tests; metaprogramming, which reduces code length while maintaining performance portability; and code generation for computing chemical reactions.

This study presents band-ensemble Spectral Proper Orthogonal Decomposition (bSPOD). The approach is inspired by frequency smoothing, a method used to reduce estimator variance in power spectral density estimates, and is here extended to SPOD. The algorithm estimates SPOD modes from consecutive Fourier coefficients obtained from a single Fourier transform of the full time record and thus avoids time segmentation. In this study, bSPOD is applied to artificial test data and to a PIV data set of a broadband-tonal cavity flow. Compared to the more commonly used Welch-based SPOD formulation, bSPOD reduces spectral leakage, permits increased frequency resolution, and retains frequency information of tonal components at comparable computational cost. These features enable reduced estimator variance while maintaining low bias for tonal components, making bSPOD particularly effective for broadband-tonal flows.

This study presents band-ensemble Spectral Proper Orthogonal Decomposition (bSPOD). The approach is inspired by frequency smoothing, a method used to reduce estimator variance in power spectral density estimates, and is here extended to SPOD. The algorithm estimates SPOD modes from consecutive Fourier coefficients obtained from a single Fourier transform of the full time record and thus avoids time segmentation. In this study, bSPOD is applied to artificial test data and to a PIV data set of a broadband-tonal cavity flow. Compared to the more commonly used Welch-based SPOD formulation, bSPOD reduces spectral leakage, permits increased frequency resolution, and retains frequency information of tonal components at comparable computational cost. These features enable reduced estimator variance while maintaining low bias for tonal components, making bSPOD particularly effective for broadband-tonal flows.

SPOD identifies coherent structures in statistically stationary flows. The method estimates the cross-spectral density matrix from an ensemble of temporal Fourier modes obtained by segmenting time series data into blocks. This statistical approach enables the extraction of coherent structures even in broadband turbulence. However, SPOD faces a tradeoff: better statistical convergence requires more blocks, which reduces frequency resolution. This variance–bias tradeoff necessitates long time series for accurate mode extraction. In this talk we present a modified version of SPOD. Instead of block-segmenting the time series, the Fourier transform is computed using the full time history. SPOD modes are then obtained by averaging neighboring Fourier modes via Proper Orthogonal Decomposition. Importantly, information about the contribution of individual Fourier modes to each SPOD mode are stored in the expansion coefficients. Compared to standard SPOD, the new approach offers two key advantages: (1) the number of Fourier modes used to estimate SPOD modes can be varied with frequency, enabling adaptive tradeoff between spectral variance and spectral bias; and (2) the contribution of individual frequencies to each SPOD mode can be directly quantified, allowing for more accurate extraction of tonal features without sacrificing statistical convergence. This makes the new algorithm particularly effective for analyzing flows that exhibit both broadband and tonal components.

Symmetry-based analyses of multiscale velocity gradients highlight that strain self-amplification (SS) and vortex stretching (VS) drive forward energy transfer in turbulent flows. By contrast, a strain–vorticity covariance mechanism produces backscatter that contributes to the bottleneck effect in the subinertial range of the energy cascade. We extend these analyses by using a normality-based decomposition of filtered velocity gradients in forced isotropic turbulence to distinguish contributions from normal straining, pure shearing and rigid rotation at a given scale. Our analysis of direct numerical simulation (DNS) data illuminates the importance of shear layers in the inertial range and (especially) the subinertial range of the cascade. Shear layers contribute significantly to SS and VS and play a dominant role in the backscatter mechanism responsible for the bottleneck effect. Our concurrent analysis of large-eddy simulation (LES) data characterizes how different closure models affect the flow structure and energy transfer throughout the resolved scales. We thoroughly demonstrate that the multiscale flow features produced by a mixed model closely resemble those in a filtered DNS, whereas the features produced by an eddy viscosity model resemble those in an unfiltered DNS at a lower Reynolds number. This analysis helps explain how small-scale shear layers, whose imprint is mitigated upon filtering, amplify the artificial bottleneck effect produced by the eddy viscosity model in the inertial range of the cascade. Altogether, the present results provide a refined interpretation of the flow structures and mechanisms underlying the energy cascade and insight for designing and evaluating LES closure models.

This study investigates the onset of linear instabilities and their later nonlinear interactions in the shear layer of an initially laminar jet using high-fidelity simulations. We present a quantitative analysis of the vortex-pairing phenomenon by computing the spatial growth rates and energy budget of the dominant frequencies. Compared with a turbulent jet, the hydrodynamic instabilities and vortex pairing are enhanced in an initially laminar jet. Using local linear theory, we identify the fundamental as the frequency with the largest spatial growth rate, and its exponential growth causes the shear layer to roll up into vortices. Visualisations and conditional $x$ – $t$ plots reveal that fundamental vortices pair to form subharmonic vortices, which then merge to produce second subharmonic vortices. The energy transfer during this process is evaluated using the spectral turbulent kinetic energy equation, focusing on dominant coherent structures identified through spectral proper orthogonal decomposition. Spectral production and nonlinear transfer terms show that the fundamental frequency gains energy solely from the mean flow, while subharmonics gain energy both linearly from the mean flow and nonlinearly through backscatter from the fundamental frequency. Our results confirm Monkewitz’s theoretical model of a resonance mechanism between the fundamental and subharmonic, which supplies energy to the subharmonic. We highlight the energetic versus dynamical importance of tonal frequencies. The second subharmonic corresponds to the largest spectral peak, while the fundamental, though the fourth largest spectral peak, is dynamically dominant, as it determines all other spectral peaks and supplies energy to the subharmonics through a reverse energy cascade.

In this paper, we present resolvent4py , a parallel Python package for the analysis, model reduction and control of large-scale linear systems with millions or billions of degrees of freedom. This package provides the user with a friendly Python-like experience (akin to that of well-established libraries such as numpy and scipy ), while enabling MPI-based parallelism through mpi4py , petsc4py and slepc4py . In turn, this allows for the development of streamlined and efficient Python code that can be used to solve several problems in fluid mechanics, solid mechanics, graph theory, molecular dynamics and several other fields.

An axisymmetrically forced initially-laminar jet at a Reynolds number $Re = \rho U_{j}D/\mu = 50000$ is investigated to identify the presence of Widnall instability. The Widnall instabilty is a short-wave azimuthal instability that grows on isolated vortex rings. Using visualizations and Q-criterion, we identify a primary vortex rollup associated with the Kelvin-Helmholtz response of the jet to forcing and secondary azimuthal perturbations that grow on the rolled-up cores that have characteristics of the Widnall instability. The presence of this instability in forced transitional jets was previously identified by Monkewitz et al [1], and our numerical simulations provide direct evidence that supports and extends those findings. Conditional averaging highlights the dominance of higher azimuthal wavenumbers, particularly $m=16$, which contribute to the wavy deformation of the vortex rings in the initial part of shear layer. As the core of the ring thickens further downstream, lower azimuthal wavenumbers become dominant. Conditional space-time proper orthogonal decomposition is used to track the evolution of the perturbations. The leading mode captures both the growth of azimuthal distortions and the gradual spacing of the rings as they evolve downstream. Notably, the vortex rings develop a four-lobed structure within $ x \leq 5$ of the shear layer, which closely matches the most amplified wavenumber predicted by the Widnall and Tsai dispersion relation [2], providing strong evidence of Widnall instability. As this instability grows on the vortex ring, it generates counter-rotating streamwise vortices that induce a lift-up effect, resulting in enhanced entrainment and the formation of streaks.

In this paper, we present resolvent4py, a parallel Python package for the analysis, model reduction and control of large-scale linear systems with millions or billions of degrees of freedom. This package provides the user with a friendly Python-like experience (akin to that of well-established libraries such as numpy and scipy), while enabling MPI-based parallelism through mpi4py, petsc4py and slepc4py. In turn, this allows for the development of streamlined and efficient Python code that can be used to solve several problems in fluid mechanics, solid mechanics, graph theory, molecular dynamics and several other fields.

Solutions to hyperbolic systems comprise waves propagating at finite speeds. When wave propagation is predominantly unidirectional, one-way wave equations can be used to evolve only the right-going solution by removing support for left-going waves. The One-Way Navier-Stokes (OWNS) approach, which was originally developed for systems of first-order hyperbolic equations, constructs one-way approximations to the linearized Navier-Stokes equations using a recursive filter to remove left-going waves. The computational cost scales with the number of recursion parameters, which must be carefully chosen to ensure accuracy and stability of the resulting one-way equation. Previous work has chosen parameters based on heuristic estimates of key eigenvalues, which requires trial-and-error tuning while also yielding slow error convergence. We propose a greedy algorithm for automatic parameter selection, which we show yields faster convergence and a net decrease in computational cost for linear and nonlinear disturbance evolution in boundary-layer flows. We review the OWNS projection (OWNS-P) and recursive (OWNS-R) methods, comparing their convergence properties, and show through our numerical analysis and experiments that OWNS-P yields superior convergence and stability properties. Although we demonstrate the method for Navier-Stokes equations, we perform our analyses on systems of linear first-order hyperbolic equations and emphasize that the greedy algorithm is applicable to such systems.

Multiphase flows are ubiquitous in both nature and engineering. Over the past two to three decades, substantial progress has been made in developing numerical methods for simulating these complex flows. Yet, significant challenges persist in accurately capturing intricate interfacial dynamics and the multi-scale interactions inherent to multiphase systems. This review focuses on several key numerical approaches that have proven particularly relevant from both practical and theoretical perspectives. In particular, we discuss Volume-Of-Fluid techniques, level set methods, diffuse interface models, and front tracking methods, along with immersed boundary strategies designed for particle-laden flows. We also examine multi-fluid Eulerian frameworks, population balance models for reactive processes, and sub-grid scale techniques for handling unresolved dynamics. Furthermore, emerging hybrid strategies that integrate conventional numerical methods with data-driven machine learning techniques are highlighted as promising directions. In conclusion, while current methodologies offer valuable insights into multiphase flow behavior, continued interdisciplinary efforts are essential to enhance predictive accuracy, computational efficiency, and the overall applicability of these simulations to real-world challenges.

Boundary-layer instabilities for a finned cone at [Formula: see text], [Formula: see text], and zero incidence angle are examined using linear stability methods of varying fidelity and maturity. The geometry and laminar flow conditions correspond to experiments conducted at the Boeing Air Force Mach 6 Quiet Tunnel at Purdue University. Where possible, a common baseflow is utilized among the stability computations, and comparisons are made along the acreage of the cone where transition is first observed in the experiment. Stability results utilizing linear stability theory, planar parabolized stability equations, one-way Navier–Stokes, forced direct numerical simulation, and adaptive mesh refinement wavepacket tracking are presented. A dominant three-dimensional vortex instability occurring at [Formula: see text] is identified and correlates well with experimental measurements of transition onset. With the exception of linear stability theory, all of the higher-fidelity linear methods considered in this work were consistent in predicting the initial growth and general structure of the vortex instability as it evolved downstream. Some of the challenges, opportunities, and development needs of the stability methods considered are discussed.

We extend the one-way Navier Stokes (OWNS) approach to support nonlinear interactions between waves of different frequencies, which will enable nonlinear analysis of instability and transition. In OWNS, the linearized Navier–Stokes equations are parabolized and solved in the frequency domain as a spatial initial-value (marching) problem. OWNS yields a reduced computational cost compared to global solvers while also conferring numerous advantages over the parabolized stability equations (PSEs), despite its higher computational cost relative to PSE, that we seek to extend to nonlinear analysis. We validate the nonlinear OWNS (NOWNS) method by examining the nonlinear evolution of two- and three-dimensional disturbances in a low-speed Blasius boundary layer compared to nonlinear PSE (NPSE) and direct numerical simulation (DNS) results from the literature. We demonstrate that NOWNS can be used to simulate flows involving blowing/suction strips, is more robust to numerical noise, and converges for stronger nonlinearities, as compared to NPSE.

Modal decomposition methods are important for characterizing the low-dimensional dynamics of complex systems, including turbulent flows. Different methods have varying data requirements and produce modes with different properties. Spectral proper orthogonal decomposition (SPOD) produces orthogonal, energy-ranked spatial modes at discrete temporal frequencies for statistically stationary flows. However, SPOD requires long stretches of sequential, uniformly sampled, time-resolved data. These data requirements limit SPOD's use in experimental settings where the maximum capture rate of a camera is often slower than the Nyquist sampling rate required to resolve the highest turbulent frequencies. However, if two PIV systems operate in tandem, pairs of data can be acquired that are arbitrarily close in time. The dynamic mode decomposition (DMD) uses this pairwise data to resolve frequencies up to the Nyquist frequency associated with the small time step within a pair. However, these modes do not form a basis and have no set ranking. The present work attempts to compute SPOD modes from pairwise data with a small time step but with large gaps between pairs. We use DMD on pairwise data to estimate segment-wise, uniformly sampled series that can then be used to estimate the SPOD modes, intending to resolve frequencies between the gap and pair Nyquist limits. The method is tested on numerically obtained data of the linearized complex Ginzburg-Landau equation, as well as a Mach 0.4 isothermal turbulent jet. For the jet, pairwise SPOD can accurately de-alias the SPOD spectrum and estimate mode shapes at frequencies up to St = 1.0, while using over 90% less data.

Wall-modeled large-eddy simulation (WMLES) has attracted significant attention in engineering applications as a high-fidelity simulation technique not limited by near-wall resolution requirements. It significantly reduces computational costs by resolving the energy-containing and dynamically important scales of turbulence far from the wall while modeling the effects of near-wall eddies. However, the inherent grid dependence of WMLES results raises the question of how to determine an adequate yet minimal resolution for the outer region. To address this question, we propose an approach to assess grid resolution that does not rely on a priori knowledge of domain-specific flow statistics. This approach leverages the velocity gradient partitioning, which provides an expressive, broadly applicable, and Galilean invariant description of local flow features. The error metric we define represents the deviation of the partitioning far from the wall in WMLES from the partitioning associated with isotropic turbulence. Using simulations of turbulent channel flow, we compare the convergence trends of this metric with those of conventional metrics representing turbulence kinetic energy errors. The partitioning metric effectively captures the response of resolved small-scale flow features to grid resolution and various other mesh and simulation parameters. It is particularly sensitive to the mesh-cell aspect ratio and subgrid-scale modeling; however, it is less sensitive to the friction Reynolds number and wall boundary conditions, which primarily impact large-scale flow features. Hence, the metric provides valuable and inexpensive insight into the sensitivity of small-scale flow features. While it alone cannot be used to holistically validate the accuracy of simulation results, coupling it with a metric that diagnoses large-scale flow features would provide a more complete picture. Beyond turbulent channel flow, preliminary results for flow over a Gaussian bump highlight the potential for the present approach to be applied to WMLES in more complex geometries.

In the context of transition analysis, traditional linear input-output analysis determines worst-case disturbances to a laminar base flow based on a generic right-hand-side volumetric/boundary forcing term. The worst-case forcing is not constrained to be physically realizable. To connect input-output analysis with physical receptivity, we use a scattered-wave ansatz developed in Kamal et al. (2023) to restrict input-output analysis to forcings that are associated with a specified set of free-stream disturbances. In this paper, we apply this technique to parametrically study the optimal free-stream waves, i.e. ones that maximally excite boundary-layer instabilities, on a Mach-6 flat-plate boundary layer across a range of frequencies, wall-temperature ratios, and spanwise wavenumbers. This parameter space gives rise to the first mode as the dominant instability at low frequencies and transitions into the second mode at higher frequencies, while also giving rise to the supersonic mode under high wall cooling. We also closely study the receptivity of the oblique first mode, second mode, and the supersonic mode and characterize the spanwise wavenumbers that yield the largest disturbance amplification.

The velocity gradient tensor can be decomposed into normal straining, pure shearing and rigid rotation tensors, each with distinct symmetry and normality properties. We partition the strength of turbulent velocity gradients based on the relative contributions of these constituents in several canonical flows. These flows include forced isotropic turbulence, turbulent channels and turbulent boundary layers. For forced isotropic turbulence, the partitioning is in excellent agreement with previous results. For wall-bounded turbulence, the partitioning collapses onto the isotropic partitioning far from the wall, where the mean shearing is relatively weak. By contrast, the near-wall partitioning is dominated by shearing. Between these two regimes, the partitioning collapses well at sufficiently high friction Reynolds numbers and its variations in the buffer layer and the log-law region can be reasonably modelled as a function of the mean shearing strength. Altogether, our results highlight the expressivity and broad applicability of the velocity gradient partitioning as advantages for turbulence modelling.

Time-resolved (TR), three-dimensional (3D) velocity fields of screeching, underexpanded jets are estimated using non-time-resolved particle image velocimetry and simultaneous TR microphone measurements. Specifically, we aim to reconstruct TR 3D velocity fluctuation fields associated with the A2, B, and C modes of a screeching jet using a linear regression model and to analyze screech dynamics of these modes. The linear regression model is constructed on the basis of a linear relationship between the velocity and acoustic fields. Three nozzle pressure ratios (NPRs) of 2.30, 2.97, and 3.40 are employed. The dominant azimuthal modes for three cases are investigated using azimuthal Fourier coefficients of the acoustic data obtained by the azimuthal array of eight microphones placed near the nozzle exit. The dominant azimuthal modes at NPRs of 2.30, 2.97, and 3.40 are <a:math xmlns:a="http://www.w3.org/1998/Math/MathML"> <a:mrow> <a:mi>m</a:mi> <a:mo>=</a:mo> <a:mn>0</a:mn> </a:mrow> </a:math> , 1, and 1, respectively. The first two proper orthogonal decomposition (POD) modes in these azimuthal modes are dominant at all NPRs and are associated with screech. 3D velocity fluctuation fields associated with screech are reconstructed from these leading POD modes of the acoustic data. The reconstructed 3D velocity fluctuation fields at NPRs of 2.97 and 3.40 exhibit two helical structures with opposite rotation directions. The present results demonstrate that, in the B mode, the flapping structure exhibits random clockwise and counterclockwise rotations over an extended time domain, while maintaining a consistent direction within short time domains. In addition, in the C mode, two helical structures with opposite rotation directions, as well as the flapping structure, are observed. Published by the American Physical Society 2024

We extend the One-Way Navier Stokes (OWNS) approach to support nonlinear interactions between waves of different frequencies, which will enable nonlinear analysis of instability and transition. In OWNS, the linearized Navier-Stokes equations are parabolized and solved in the frequency domain as a spatial initial-value (marching) problem. OWNS yields a reduced computational cost compared to global solvers, while also conferring numerous advantages over the parabolized stability equations (PSE), despite its higher computational cost relative to PSE, that we seek to extend to nonlinear analysis. We validate the nonlinear OWNS (NOWNS) method by examining nonlinear evolution of two- and three-dimensional disturbances in a low-speed Blasius boundary layer compared to nonlinear PSE (NPSE) and direct numerical simulation (DNS) results from the literature. We demonstrate that NOWNS can be used to simulate flows involve blowing/suction strips, is more robust to numerical noise, and converges for stronger nonlinearities, as compared to NPSE.

We solve the incompressible Navier-Stokes equations using a lattice Green's function (LGF) approach, including immersed boundaries (IB) and adaptive mesh refinement (AMR), for external flows with one homogeneous direction (e.g. infinite cylinders of arbitrary cross-section). We hybridize a Fourier collocation (pseudo-spectral) method for the homogeneous direction with a specially designed, staggered-grid finite-volume scheme on an AMR grid. The Fourier series is also truncated variably according to the refinement level in the other directions. We derive new algorithms to tabulate the LGF of the screened Poisson operator and viscous integrating factor. After adapting other algorithmic details from the fully inhomogeneous case, we validate and demonstrate the new method with transitional and turbulent flows over a circular cylinder at $Re=300$ and $Re=12,000$, respectively.

A two-dimensional body that moves suddenly in a viscous fluid can instantly generate vortices at its sharp edges. Recent work using inviscid flow theory, based on the Birkhoff–Rott equation and the Kutta condition, predicts that the ‘starting vortices’ generated by the sharp and straight edges of a body – i.e. the vortices formed immediately after motion begins – can be one of three distinct self-similar types. We explore the existence of these starting vortices for a flat plate and two symmetric Joukowski aerofoils immersed in a viscous fluid, using high-fidelity direct numerical simulations (DNS) of the Navier–Stokes equations. A lattice Green's function method is employed and simulations are performed for chord Reynolds numbers ranging from 5040 to 45 255. Vortices generated at the leading and trailing edges of the flat plate show agreement with the derived inviscid theory, for which a detailed assessment is reported. Agreement is also observed for the two symmetric Joukowski aerofoils, demonstrating the utility of the inviscid theory for arbitrary bodies. While this inviscid theory predicts an abrupt transition between the starting-vortex types, DNS shows a smooth transition. This behaviour occurs for all Reynolds numbers and is related to finite-time effects – there is a maximal time for which the (self-similar) starting vortices exist. We confirm the inviscid prediction that the leading-edge starting vortex of a flat plate can be suppressed dynamically. This has implications for the performance of low-speed aircraft such as model aeroplanes, micro air vehicles and unmanned air vehicles.

This work seeks rapid prediction of laminar-to-turbulent transition based on first principles, rather than commonly employed empirical correlations. The nonlinear one-way Navier-Stokes (NOWNS) equations were recently applied to the early stages of boundary layer transition where it was demonstrated that NOWNS can accurately replicate direct numerical simulation (DNS) results with similar accuracy to the nonlinear parabolized stability equations (NPSE). While having greater computational cost than NPSE, NOWNS is a more robust, convergent parabolization of the governing equations and is expected to succeed for stronger nonlinearity where NPSE fails. We demonstrate that NOWNS succeeds for stronger nonlinearity by applying it to an oblique-wave breakdown case where NPSE fails. In addition, we demonstrate that NOWNS supports non-modal disturbances in the form of random noise applied to the inlet boundary condition of the oblique-wave breakdown case, and a blowing/suction strip for a K-type (fundamental) transition case.

The velocity gradient tensor can be decomposed into normal straining, pure shearing and rigid rotation tensors, each with distinct symmetry and normality properties. We partition the strength of turbulent velocity gradients based on the relative contributions of these constituents in several canonical flows. These flows include forced isotropic turbulence, turbulent channels and turbulent boundary layers. For forced isotropic turbulence, the partitioning is in excellent agreement with previous results. For wall-bounded turbulence, the partitioning collapses onto the isotropic partitioning far from the wall, where the mean shearing is relatively weak. By contrast, the near-wall partitioning is dominated by shearing. Between these two regimes, the partitioning collapses well at sufficiently high friction Reynolds numbers and its variations in the buffer layer and the log-law region can be reasonably modelled as a function of the mean shearing strength. Altogether, our results highlight the expressivity and broad applicability of the velocity gradient partitioning as advantages for turbulence modelling.

When the jet noise is cast in the form of a linear input-output problem, the cross-spectral-density (CSD) matrix of the sound field is forced by the CSD of non-linear interactions via the resolvent operator. The forcing CSD is difficult to compute or measure, but its projection onto the resolvent input space can be identified from the CSD of the sound field. In a recent study using data from numerical simulation, such an identification showed that a low-rank truncation of the forcing-projection matrix can be used to reconstruct the acoustic field of transonic and supersonic turbulent jets with Mach numbers 0.9 and 1.5 respectively. Encouraged by this result, we performed two-point acoustic measurements in isothermal turbulent jets over a broad range of subsonic Mach numbers. In the present work, these CSD matrices are used to obtain an empirical model of the forcing-projection matrix, following the procedure reported in Pickering et al., 2021. We investigate the parametric dependence of the low-rank, forcing-projection matrix on Mach number, frequency, and azimuthal mode, and we propose a model that captures this dependence and, allowing computation of downstream radiation with 1.5dB precision for subsonic jets in the Mach-number range 0.7 to 0.9.

Modes from global resolvent analyses have been shown to accurately model the frequencies and spatial structure of the dominant coherent structures in several turbulent flows. However, resolvent-mode forcing models must be developed to predict the amplitude of the structures or other flow statistics, including the radiated noise. The present research aims to apply data-driven approaches to learn forcing coefficients from lower-order statistics available from Reynolds-averaged Navier-Stokes (RANS) predictions. As a first step towards this goal, we present a novel localized resolvent framework that reconstructs global quantities at low rank through spatially restricting the resolvent forcing and response domains. To illustrate the flexibility and robustness of the proposed framework, we initially utilize localized resolvent modes to reconstruct the spectral proper orthogonal decomposition (SPOD) modes of an isothermal Mach 0.4 jet at ���� = 450, 000. The results showcase the flexibility localized resolvent modes provide in the construction of global SPOD, while using 10 or fewer total localized modes total across ���� = [0.05, 1.00]. Furthermore, we employ localized resolvent modes to reconstruct second-order statistics, comparing their performance with that of global modes. At low reconstruction error, it is shown that about twice as many global modes are needed to achieve comparable errors.

Bispectral mode decomposition is used to investigate triadic interactions within a Mach 0.4 turbulent jet. We explore its potential to identify dominant triadic interactions and their associated spatial structures in an unforced turbulent jet. The bispectral measure is broadband in frequency for each azimuthal wavenumber triad. The [1,1,2] and [0,0,0] azimuthal wavenumber triads are dominant, emphasizing the importance of the self-interactions of the helical and axisymmetric components. Bispectral modes reveal that streaky structures are fed by the interaction of a Kelvin-Helmholtz wavepacket with its conjugate. Streaks are also observed in other frequency interactions, occurring in regions where the structures of these frequencies are spatially active. Furthermore, integral interaction maps and nonlinear transfer terms are computed to determine the direction of energy transfer and to pinpoint the spatial regions where nonlinearity is most active. As the shear layer develops, small scales interact nonlinearly, transferring energy to larger scales. Moving downstream, near the potential core closure, larger scales become more active, resulting in a forward energy cascade.

Coherence decay has been understood to be a key quantity to predict acoustic noise emitted by wavepackets in subsonic turbulent jets. Frequency-domain frameworks such as input-output and resolvent analyses are able to predict accurately the spatial structure of wavepackets turbulent flows compared to coherent structures educed from simulation data (as for example identified using spectral proper orthogonal, SPOD). However, at least at reduced-order, they are unable to capture two-point statistics such as coherence. A missing piece is the modelling of variability induced by the turbulence, which jitters (disorganises) the coherent structures and leads to stronger noise radiation. The aim of the present study is to consider the impact of turbulence on jet wavepackets through stochastic modelling under location uncertainty. This framework considers the conservation of mass and momentum of fluid parcels submitted to a stochastic transport, representing here the effect of turbulence. By linearising the resulting generalised stochastic Navier-Stokes equations and expressing it in the Fourier domain, a stochastic linear model (SLM) is obtained. We explore in this paper that ability of SLM to predict the two point coherence of the wavepackets in turbulent jets, and show its impact on acoustic emissions.

Multiphase compressible flows are often characterized by a broad range of space and time scales. Thus entailing large grids and small time steps, simulations of these flows on CPU-based clusters can thus take several wall-clock days. Offloading the compute kernels to GPUs appears attractive but is memory-bound for standard finite-volume and -difference methods, damping speed-ups. Even when realized, faster GPU-based kernels lead to more intrusive communication and I/O times. We present a portable strategy for GPU acceleration of multiphase compressible flow solvers that addresses these challenges and obtains large speedups at scale. We use OpenACC for portable offloading of all compute kernels while maintaining low-level control when needed. An established Fortran preprocessor and metaprogramming tool, Fypp, enables otherwise hidden compile-time optimizations. This strategy exposes compile-time optimizations and high memory reuse while retaining readable, maintainable, and compact code. Remote direct memory access, realized via CUDA-aware MPI, reduces communication times. We implement this approach in the open-source solver MFC. Metaprogramming-based preprocessing results in an 8-times speedup of the most expensive kernels, 46% of peak FLOPs on NVIDIA GPUs, and high arithmetic intensity (about 10 FLOPs/byte). In representative simulations, a single A100 GPU is 300-times faster than an Intel Xeon CPU core, corresponding to a 9-times speedup for a single A100 compared to the entire CPU die. At the same time, near-ideal (97%) weak scaling is observed for at least 13824 GPUs on Summit. A strong scaling efficiency of 84% is retained for an 8-times increase in GPU count. Collective I/O, implemented via MPI3, helps ensure negligible contribution of data transfers. Large many-GPU simulations of compressible (solid-)liquid-gas flows demonstrate the practical utility of this strategy.

Many turbulent flows exhibit time-periodic statistics. These include turbomachinery flows, flows with external harmonic forcing and the wakes of bluff bodies. Many existing techniques for identifying turbulent coherent structures, however, assume the statistics are statistically stationary. In this paper, we leverage cyclostationary analysis, an extension of the statistically stationary framework to processes with periodically varying statistics, to generalize the spectral proper orthogonal decomposition (SPOD) to the cyclostationary case. The resulting properties of the cyclostationary SPOD (CS-SPOD for short) are explored, a theoretical connection between CS-SPOD and the harmonic resolvent analysis is provided, simplifications for the low and high forcing frequency limits are discussed, and an efficient algorithm to compute CS-SPOD with SPOD-like cost is presented. We illustrate the utility of CS-SPOD using two example problems: a modified complex linearized Ginzburg–Landau model and a high-Reynolds-number turbulent jet.

We study coherent structures in subsonic turbulent jets subject to a flight stream. A thorough characterisation of the effects of a flight stream on the turbulent field was recently performed by Maia et al. ( Phys. Rev. Fluids , vol. 8, 2023, 063902) and fluctuation energy attenuations were observed over a broad range of frequencies and azimuthal wavenumbers. The Kelvin–Helmholtz, Orr and lift-up mechanisms were all shown to be weakened by the flight stream. Here we expand upon that study and model the changes in the dynamics of jets in flight using global resolvent analysis. The resolvent model is found to correctly capture the main effects of the flight stream on the dynamics of coherent structures, which are educed from a large-eddy simulation database using spectral proper orthogonal decomposition. Three modifications of note are: the damping of low-frequency streaky/Orr structures that carry most of the fluctuation energy; a degradation of the low-rank behaviour of the jet in frequencies where modal instability mechanisms are dominant; and a rank decrease at very low Strouhal numbers. The latter effect is underpinned by larger gain separations predicted by the resolvent analysis, due to a reduction in the wavelength of associated flow structures. This leads to a clearer relative dominance of streaky structures generated by the lift-up mechanism, despite the fact that the lift-up mechanism has been weakened with respect to the static jet.

We simulate the head-on collision between vortex rings with circulation Reynolds numbers of 4000 using an adaptive, multiresolution solver based on the lattice Green's function. The simulation fidelity is established with integral metrics representing symmetries and discretization errors. Using the velocity gradient tensor and structural features of local streamlines, we characterize the evolution of the flow with a particular focus on its transition and turbulent decay. Transition is excited by the development of the elliptic instability, which grows during the mutual interaction of the rings as they expand radially at the collision plane. The development of antiparallel secondary vortex filaments along the circumference mediates the proliferation of small-scale turbulence. During turbulent decay, the partitioning of the velocity gradients approaches an equilibrium that is dominated by shearing and agrees well with previous results for forced isotropic turbulence. We also introduce new phase spaces for the velocity gradients that reflect the interplay between shearing and rigid rotation and highlight geometric features of local streamlines. In conjunction with our other analyses, these phase spaces suggest that, while the elliptic instability is the predominant mechanism driving the initial transition, its interplay with other mechanisms, e.g. the Crow instability, becomes more important during turbulent decay. Our analysis also suggests that the geometry-based phase space may be promising for identifying the effects of the elliptic instability and other mechanisms using the structure of local streamlines. Moving forward, characterizing the organization of these mechanisms within vortices and universal features of velocity gradients may aid in modelling turbulent flows.