This paper introduces a new algorithm to improve the accuracy of numerical phase-averaging in the oscillatory, multiscale, differential equations that govern geophysical fluid flow. Phase-averaging is a timestepping method that averages a mapped modulation variable in which the equations no longer contain a highly oscillatory linear term. This allows phase-averaging to retain the main contribution of fast waves on the low frequencies without explicitly resolving the rapid oscillations, thus allowing larger timesteps for explicit timesteppers. Whilst the modulation variable mapping is also used in operator splitting methods, and is accurate for asymptotic regimes, the advantage of phase-averaging is that low-frequency oscillations are retained, allowing for accuracy in systems further from the asymptotic limit. However, this comes at the cost of introducing an averaging error when computing phase-averages using a smooth kernel. This paper proposes a method to offset some averaging error through a modified mapping that includes a mean correction term to capture an average measure of nonlinear oscillation. We show that the mean corrected algorithm reduces phase-averaging errors in the swinging spring system and the one-dimensional rotating shallow water equations, which are an important geophysical fluid system for weather and climate applications. We discuss potential new directions for the method based on the outcome of the numerical experiments.
As a starting point to gain an understanding of settling porous objects with complex geometries, such as marine snow and biofouled microplastics, we study the settling of a thin porous mesh. Motivating the theoretical-numerical work presented here are laboratory experiments that show the meshes settle without tumbling and the settling speed is effectively independent of the mesh lateral extent. The numerical simulations model flow through the mesh using a time-relaxation method. This captures the qualitative features of the laboratory experiments and gives insight into the dynamics governing the steady state settling velocity. Specifically it shows that a steep pressure gradient develops across the thickness of the mesh equal in magnitude to the reduced gravity, which itself equals the difference of the fluid and mesh velocity divided by the relaxation time scale. In effect, fluid is sucked through the mesh which provides the balance between buoyancy forces and drag controlled primarily by the mesh thickness, not its lateral extent. This insight justifies a dimensional analysis approach which successfully collapses measurements of the steady state flow across the mesh and mesh settling speeds thus giving semi-empirical predictions for these speeds as they depend upon the reduced gravity, mesh thickness and relaxation time scale.
The efficient and accurate rational approximation of the Schwarz function in the complex plane by means of the recently developed AAA (adaptive Antoulas-Anderson) algorithm is applied here to the dynamics of uniform vorticity patches. The key observation is realising that the AAA algorithm readily computes the poles of the rational approximation, and that this enables a natural partitioning of the Schwarz function in such a way that the velocity field exterior to a patch is given in terms of the poles interior to the patch. In particular, knowing the velocity on the patch boundary and exterior to the patch allows both vortex patch equilibria, including configurations with and without satellite point vortices surrounding the vortex patch, and time-dependent motions to be computed. Examples are presented and the performance of the new numerical approach compared to previously established numerical methods, including the method of contour dynamics and the recently developed AAA-least squares algorithm.
A new vortex-dominated flow problem was defined to reveal some flow physics and to aid the design and assessment of turbulence models. The case is rich in turbulence and vorticity and physically complex although two-dimensional in the mean and free of solid boundaries. The z direction is homogeneous, and used for averages in the DNS so that the turbulence quantities are functions of (x, y, t). The flow is periodic in the x direction. The initial condition is a thin shear layer straddling the (x, z) plane and given a slight wave in the (x, y) plane. It is seeded with random perturbations, and turbulence develops before the Kelvin-Helmholtz instability creates two primary vortices. Each vortex entrains the shear layer into spiral sheets around itself, like vortices do over wings. Most RANS models are inaccurate in these situations, which is confirmed in the present case; they miss how the turbulence rapidly decays in and near the young vortex. The two vortices later collide and merge. The flow is driven by concurrent inviscid and turbulent phenomena, and RANS is more accurate for the former than the latter. A variety of measures of the numerical quality of the DNS are presented. We seek evidence that the vorticity exceeds the bounds of its initial condition. This cannot happen in (2D) laminar flow, but is common in RANS: vorticity of the sign opposite to the initial condition is generated by many models. We see strong evidence that this is not meaningfully present in the DNS, and attribute the small opposite excursions to the finite spatial-averaging interval and residual numerical errors. This matches our intuition, although we know of no theorem stating that opposite vorticity is unphysical. Part II contains comparisons with models, including eddy viscosity with and without rotation corrections and Reynolds-stress transport.
To predict the amplification factor of the Mack second mode instability in a hypersonic boundary layer, an amplification factor transport equation was developed within a RANS framework. Its production term is based on the approximate N-envelope method and modeled using compressible boundary layer profiles and linear stability theory (LST) analysis data. The transport equation is designed to use only local variables and does not require wall temperature as an explicit parameter. Furthermore, it accounts for the influence of wall temperature on the growth of the Mack second mode instability, which is a key contribution of this study. The amplification factor and the modified turbulence intermittency transport equations are coupled with a k- ω SST turbulence model as the transition model. The model was implemented in SU2, and grid experiments were performed. The predicted amplification factor was compared with the LST results with differences typically ranging from 0.4
Acoustic-gravity waves (AGWs) are sound waves in the ocean slightly modified by a gravity wave boundary condition at the surface. Typical frequencies are in the range 0.1-1 Hz . In the study of such waves, the acoustic properties of the seabed are usually not discussed, and the sea floor is taken to be insulating to sound waves. In this paper, we allow the seabed to be acoustically conducting such that a small part of the acoustic wave energy may propagate though the sea floor. To model this process in a principal way, we apply a Robin boundary condition for the velocity potential. Because of the energy loss at the bottom, the AGW motion becomes slightly damped in the fluid. Nonlinearly, this leads to a nonzero forcing term for the mean Eulerian velocity (the acoustic streaming), when averaged over the wave period. We utilize the fact that the planetary rotation is very slow compared with the AGW period and derive the Coriolis-Stokes force for the mean Eulerian motion in the cross-wave direction. We use that the mean Lagrangian (particle) motion is the sum of the Eulerian mean motion and the Stokes drift. We then find, by averaging over the inertial period, that the Stokes drift in the wave propagation direction is balanced by a negative Eulerian mean current (an anti-Stokes flow) at all depths. The mean Lagrangian drift in the cross-wave direction is a result of a balance between the Coriolis force and the gradient of the mean Reynolds stresses, constituting the Eulerian acoustic forcing in this case. It appears that the associated Lagrangian mass flux could be of importance for the circulation in the deep ocean.
This investigation examines instability mechanisms and turbulent transition in the stationary disk boundary layer of a rotor-stator cavity through an integrated approach combining two-dimensional direct numerical simulations (2D DNS), global linear stability analysis, and three-dimensional direct numerical simulations (3D DNS). The base flow characteristics, the global linear stability properties, and the nonlinear evolution processes leading to turbulence are explored. The research examines an annular, enclosed rotor-stator cavity with curvature parameter Rm = (b^*+a^*)/(b^*-a^*)=1.8 , and an aspect ratio L=(b^*-a^*)/(2h^*) = 5 , where a^* and b^* are the minimum and maximum radius, respectively, and 2h^* denotes the inter-disk spacing. The current study examines Reynolds numbers Re_φ = Ω ^*_db^*2/ν ^* ranging from 10^4 to 10^5 , where Ω ^*_d is the angular velocity of the rotor, ν ^* is the fluid’s kinematic viscosity. The results reveal that for Re_φ≤ 3.6 × 10^4 , the base flow converges to a stable equilibrium state. However, when Re_φ exceeds 3.7 × 10^4 , the base flow destabilizes, giving rise to persistent velocity oscillations within the stationary disk boundary layer. Global linear stability analysis confirms these observations. For axisymmetric disturbances (azimuthal wavenumber β = 0 ), the critical Reynolds number for global linear instability is Re_φ = 3.67 × 10^4 , precisely demarcating the threshold for self-sustained oscillations in the base flow. For non-axisymmetric disturbances with β = 15 , a more complex dynamics is identified: the maximum linear temporal growth rate occurs at Re_φ = 2.5 × 10^4 , while the critical threshold for global linear instability emerges at a lower Re_φ = 1.58 × 10^4 , marking the onset conditions for a self-sustained spiral wave mode. The 3D DNS results reveal sophisticated nonlinear dynamics where spiral wave energy is transferred from higher to lower radial positions through complex mode interactions. When circular waves dominate in the flow field, their interaction with spiral waves across various radial positions generates a rich spectral composition with multiple distinct dominant frequencies within the boundary layer. At Re_φ≥ 7.0 × 10^4 , the flow transitions to a chaotic state with turbulent characteristics. Remarkably, even in this strongly chaotic regime, coherent circular wave structures remain observable within the chaotic flow field, although distinct dominant frequencies are no longer present in the spectral analysis.
Compression corner flows exhibit a shock-wave/boundary-layer interaction that induce flow separation and potential laminar-to-turbulent transition that may produce high surface heating in the reattachment region. Laminar separated flows can sustain the amplification of global instabilities that can be the origin for boundary-layer transition. A computational study of double wedge and cone-flare geometries at zero degrees angle of attack is performed to investigate the effects of geometrical parameters, wall temperature, and freestream conditions on the laminar separation and its global instability characteristics, with the objective of deriving a simple correlation between the separated flow properties and the global instability onset. Regardless of the nosetip/leading edge bluntness, the global instability is first dominated by a stationary mode with low azimuthal wavenumber that extends across the entire separated region. To examine the correlation between separation bubble strength and the onset of global instability, we created and analyzed a database spanning a wide range of parameters, including freestream Mach number, Reynolds number, surface temperature, nosetip/leading edge radius, an intermediate cylindrical section ahead of the flare, radius of curvature at the cone-flare corner, and flow turning angles. Multiple normalized quantities are investigated as potential metrics toward an empirical criterion for global instability onset. Generally, the separation bubble becomes globally unstable when the maximum reverse streamwise velocity exceeds 10
A species-specific multi-vibrational model for hypersonic non-equilibrium flows is presented in this work, along with a systematic analysis of its impact across multiple cases, ranging from zero-dimensional adiabatic relaxation to rapid expansion flows and axisymmetric blunt-body configurations. This model assumes a single temperature to characterise both translational and rotational energies, while a separate vibrational temperature is considered for each species. It explicitly accounts for vibrational-translational (V-T) relaxation of individual species and vibrational-vibrational (V-V) relaxation between unlike molecules. The V-V relaxation probabilities for N_2 - O_2 and N_2 -NO collisions are determined by fitting experimental and Schwartz-Slawsky-Herzfeld (SSH) data from the literature, whereas, the probability of O_2 -NO is assumed to follow the N_2 -NO system. The proposed model is implemented in an in-house hypersonic computational fluid dynamics (CFD) solver and validated against the numerical simulations for a non-reacting one-dimensional system, and compared with the high-fidelity state-to-state (STS) and direct simulation Monte Carlo (DSMC) approaches for a reacting flow. Further validation is performed against experimental measurements and DSMC for rapid expansion flow and flow over a cylinder. Comparisons with the single-vibrational temperature model reveal noticeable differences in vibrational temperature profiles, shock stand-off distance and vibrational freezing values. The additional computational cost associated with this formulation is also detailed, supporting its applicability as an intermediate-fidelity modelling approach for hypersonic CFD simulations.
Analytical models are among the most direct and efficient methods for predicting jet behavior in engineering and practical applications. Based on Reichardt’s hypothesis, the momentum transformation equation for a single free jet was derived. Using a superposition approach, the momentum distribution for multiple jets was determined, enabling the development of a three-dimensional analytical model for multiple parallel jets. The effectiveness of this superposition technique in predicting the mean streamwise velocity components of multiple jets was validated through a combination of experimental and numerical simulations. For non-parallel jets, interactions between jets result in deflections as they enter the flow field, introducing the concept of a“velocity deflection interface”. Due to experimental limitations, accurately determining the position of the velocity deflection interface is challenging, making numerical simulations the preferred approach. The Shear-Stress transport (SST) k- turbulence model was employed, and jet angles between 1 ^∘ and 20 ^∘ , commonly used in industrial applications, were analyzed. By fitting data along the flow direction, the position of the velocity deflection interface was identified, and its relationship with the jet angle was established. In the upstream region of the velocity deflection interface, the relationship between the velocity deflection interface and the initial jet exit characteristics was quantified. In the downstream region, the parallel jet theory based on Reichardt’s hypothesis was extended to derive analytical equations for three-dimensional non-parallel jets. Finally, a comparison of jet analysis results with experimental data and computational fluid dynamics (CFD) simulations confirmed the validity of the analytical model. It provides theoretical support for predicting the behavior of multiple jets in engineering applications.
In this study, we investigate how the asymmetric wetting conditions affect the droplet dynamics in microfluidic T-junctions using a three-dimensional multicomponent lattice Boltzmann method based on the two-component Shan-Chen model. The model has been validated against theoretical and benchmark cases, including the Young-Laplace law and classical symmetric droplet breakup. Through a series of simulations, we have explored the effects of contact angle contrast, droplet length, capillary number, channel geometry, and viscosity ratio on the margins that split breakup and non-breakup regimes. The results demonstrate that spatial wettability asymmetry significantly alters droplet evolution, which promotes controlled redirection or suppression of breakup. Increasing the upper-lower contact angle difference promotes droplet steering toward the more wettable side, whereas higher capillary numbers and longer droplets tend to experience more breakup. Additionally, geometric parameters such as aspect ratio and side-to-main channel width ratio regulate the competition between surface and hydrodynamic forces, critically influencing droplet final state. The viscosity ratio affects the droplet dynamics in two ways: impacting internal recirculation and interfacial stress resistance. This work aims to provide a physics-based framework for passive droplet control by utilizing wettability engineering and geometrical design to offer practical insights for lab-on-a-chip devices and multiphase flow systems. Future work will incorporate non-Newtonian effects and dynamic wetting to extend applicability to more complex flow regimes.
In this work, a new model-variant in the Bautista-Manero-Puig (BMP) model is proposed, with the peculiarity of providing true yield-stress features under vanishing deformation rates. This proposed model is devised to predict an Elasto-Visco-Plastic (EVP) rheological response through the assumption of a null fluidity (infinite viscosity) limit, and inherits the BMP ability to track the material-structure temporal evolution through a thixotropic kinetic equation. This new model, termed BMP-EVP, is systematically compared against the parent original BMP model, which provides solid-like features in an apparent yield-stress modality, through diminishing solvent fractions. The rheological fingerprints of both BMP and BMP-EVP models are exposed and contrasted under typical steady and transient rheometrical tests, under shearing and extensional deformations, alongside oscillatory protocols under Small and Large Angle Oscillatory Shear (SAOS and LAOS, respectively). Here, plasticity is measured through the definition of critical deformations for material fluidisation that appeared explicitly and directly correlated with the BMP thixotropic parameters. Conventional steady-state simple shear and uniaxial extensional flows render the BMP and BMP-EVP models comparable in their material-function predictive capabilities. In contrast, non-linear transient tests show attractive differences across these models, for which secondary loops in LAOS and non-monotonic stress-growth coefficient in start-up flows, evidence the interplay between the thixo-viscoelastic features of these models and their distinct plastic responses.
In this study, the spectral difference method is implemented for the numerical solution of the discrete Boltzmann equation to provide a robust and high-order accurate compressible gas kinetic scheme for effectively computing compressible rarefied gas flows. To this aim, the discrete Boltzmann equation with the Shakhov model is considered and the spatial discretization in the resulting equation is performed by the spectral difference method and the fourth-order Runge-Kutta method is used for the temporal discretization. Different one- and two-dimensional test cases are simulated to examine the accuracy and performance of the present methodology based on the spectral difference solution of the Boltzmann equation (SDBE). At first, the two-dimensional incompressible flow problems, namely, the Taylor-Green vortex flow and the cavity flow are simulated by applying the third-order SDBE and the results obtained are compared with the analytical solution and the available gas-kinetic and direct simulation Monte Carlo (DSMC) results which exhibit good agreement. Then, some compressible flow problems including the one-dimensional Riemann shock tube, the one-dimensional normal shock structure and the supersonic flow over a cylinder problem are computed to better examine the accuracy and robustness of the present method by applying the SDBE in different conditions. To further assess the accuracy and performance of the third-order SDBE, the simulations are also performed by the third-order upwind finite-difference solution of the Boltzmann equation (UFDBE) and the results of these two numerical schemes are thoroughly compared with each other. It is indicated that the high-order gas kinetic scheme implemented based on the spectral difference solution of the Boltzmann equation (SDBE) can be applied for accurately and effectively computing rarefied gas flows in a wide range of Knudsen numbers.
We propose an immersed boundary method (IBM) with a curvature-dependent, area-preserving correction algorithm for binary immiscible incompressible fluid flows. The IBM was first proposed to solve biofluid dynamics in complex geometries, and researchers later used it for multiphase fluid flows due to its direct and simple representation of complex interfaces. In this method, two types of grids are needed: an Eulerian formulation for the computation of fluid flow and a Lagrangian representation for the movement of the immersed boundary. If the conventional method is used, area loss or area increases occur due to numerical discretization errors. In the proposed correction method, the positions of the Lagrangian interface points are adjusted along the normal direction in proportion to the local curvature. Numerical examples of droplet deformation under various flow conditions show that the proposed algorithm can discretely preserve the initial area while the interface undergoes large deformation.
We use three-dimensional numerical simulations based on the lattice Boltzmann method to study how an ellipsoidal microscopic swimmer moves through a square microchannel. Three key factors are varied: the strength of inertial effects in the flow, the swimmer’s shape (from spherical to three times longer than it is wide), and the degree of confinement by the channel walls (from narrow channels about three swimmer diameters across to wider channels about eight diameters). We consider both pushers (which drive the fluid backward with their rear end like sperm cells) and pullers (which pull the fluid forward with their front end like algae). As inertia increases, pushers swim faster whereas pullers slow down, and the change in speed is much stronger for pushers. The swimming speed can be captured by a simple quadratic trend when expressed in terms of a single combined measure of inertia and swimming stroke. More elongated swimmers move faster overall and are less influenced by inertia, while nearly spherical swimmers are the most sensitive to changes in inertia. As the channel becomes wider, the walls constrain the swimmer less, and variations in inertia have a more pronounced impact on the swimming speed.
Rayleigh-Bénard convection is a canonical problem in fluid mechanics, where an adverse temperature gradient between the opposing boundaries induces instabilities that drive natural convection. For viscoplastic fluids, a minimal perturbation is required to initiate the flow. This study presents a numerical investigation of the three-dimensional Rayleigh-Bénard convection within a cubic cavity heated from below, with lateral walls subject to a linear temperature profile. The fluid behavior is modeled using the Bingham constitutive model. The moment-based Lattice Boltzmann Method was employed as the numerical method to solve the mass and momentum transport equations, with an extended formulation to incorporate the energy transport equation using a local diffusion coefficient approach. Simulations are performed for Rayleigh numbers between 104 and 107. Within this range, we observed a region of Yield numbers, between 0.004 and 0.007, that fluid plastifies. Increasing the Rayleigh number led to a transition from a stationary to a chaotic state, while larger Prandtl numbers damped the fluctuations in the velocity field. Notably, the imposition of a linear temperature profile on the lateral boundaries enhances flow instability, thereby amplifying plastic instabilities as the critical Yield number is approached.
We present M-OWNS, a spatial marching method that combines the carrier-wave factoring of the parabolised stability equations (PSE) with a recursive one-way Navier–Stokes (OWNS-R) projection framework. A distinct numerical resolution and efficiency advantage is offered by the approach, in modelling disturbance and instability state evolution. A spectral resolution comparison analysis shows that to leading order, for any excited eigenfunction whose eigenvalue lies closer to the carrier wavenumber than to the origin, the wave-factored system resolves the mode at a coarser streamwise numerical step size relative to the unfactored system. A non-iterating variant, with the carrier wavenumber determined from the base flow, temporal frequency and spanwise wavenumber alone, achieves equivalent resolution accuracy at identical per-step cost to unfactored OWNS. For the fixed-carrier variant, M-OWNS reduces the total solve count by factors of two to eight relative to unfactored OWNS across the test cases considered, with larger reductions possible when the iterated closure condition of PSE is suitable. The method is validated across incompressible and subsonic flat-plate boundary-layers, three-dimensional crossflow disturbances, and a Mach 4.5 hypersonic boundary-layer with four forcing configurations: eigenfunction inlet forcing, wall suction/blowing, multi-mode freestream forcing and randomised inlet forcing. The wall suction/blowing case is validated against a fully elliptic linear harmonic Navier–Stokes solver. For deterministic forcing scenarios, M-OWNS captures disturbance amplitudes, acoustic radiation fields, and modal synchronisation sequences at coarser streamwise resolution than unfactored OWNS. Under broadband randomised forcing, M-OWNS resolves mixed-mode disturbance development at half the numerical cost relative to standard OWNS.
We develop a model for steady, laminar boundary layers over small-scale textured surfaces. Although the texture is small relative to the boundary-layer thickness, it modifies the flow via a slip length. We use matched asymptotic expansions to simplify the problem, dividing the flow into outer, boundary-layer and inner regions. The far-field behaviour of the inner problem yields a slip boundary condition for the boundary layer. We derive an asymptotic solution valid when the slip length is small, and for arbitrary slip lengths, we develop a numerical method combining Chebyshev collocation and finite differences. We apply this framework to canonical small-scale textured surfaces, including superhydrophobic surfaces and riblets, and utilise existing analytical slip formulae. However, the framework is expected to extend to liquid-infused, porous, compliant or deformable surfaces with a variety of regular or random textures. We demonstrate how slip modifies the boundary layer’s velocity field, wall shear stress and displacement thickness across a range of surface configurations, and examine the linear stability of the resulting slip-modified boundary layers. Our approach enables computationally inexpensive modelling of a wide range of small-scale textured surfaces within laminar boundary-layer flows, providing predictive capability for drag, boundary-layer growth and transition across applications ranging from microfluidics to turbo-machinery and marine transport.
This study investigates the impact of elasticity and plasticity on two-dimensional flow past a circular cylinder at Reynolds number Re = 100 . Ten direct numerical simulations were performed using the Saramito-Herschel–Bulkley model to represent viscoelastic and elastoviscoplastic (EVP) fluids. The flow evolves from a periodic von Kármán vortex street to chaotic-like regimes. Proper Orthogonal Decomposition (POD) and Higher Order Dynamic Mode Decomposition (HODMD) are applied to extract dominant flow structures and their temporal dynamics. For viscoelastic fluids, increasing the Weissenberg number Wi elongates the recirculation bubble and shifts it downstream, resulting in more intricate but still periodic behavior. In EVP fluids, seven cases explore variations in Bingham number Bn, solvent viscosity ratio β _s , and power law index n, aiming to qualitatively assess their influence rather than determine critical thresholds. Results indicate that stronger plastic effects, especially with n ≥ 1 , lead to increased flow complexity. Three dynamic regimes are identified: (i) periodic; (ii) transitional, with elongated recirculation and disrupted periodicity; and (iii) fully complex, with breakdown of recirculation. Overall, the study highlights the interplay between inertia, elasticity, and yield stress in non-Newtonian flows past obstacles and identifies key parameters driving the transition from periodic to complex regimes.
This study explores the effects of a non-conductive gas layer flowing concurrently with a conductive liquid on the two-phase flow characteristics in wide horizontal ducts under a constant vertical magnetic field. To this end, analytical solutions for the velocity profile and induced magnetic field are presented for laminar gas-liquid stratified magnetohydrodynamic (MHD) flow between two infinite plates of various conductivities. The contributions of the Lorentz force and wall shear stresses to the pressure gradient are examined. To the best of our knowledge, it is shown for the first time that, unlike the single-phase Hartmann flow, the velocity profiles in two-phase flow differ significantly depending on whether the bottom wall is conducting or insulating. In the case of an insulating bottom wall, the gas lubrication effect and potential pumping power savings are significantly greater, regardless of the magnetic Reynolds number. This conclusion also holds for gas-liquid MHD flows in rectangular ducts with finite width-to-height aspect ratios. To assess the applicability of the Two-Plate (TP) model to wide ducts, numerical solutions of the two-dimensional problem are used to investigate the influence of side walls on the two-phase flow characteristics, considering various combinations of bottom and side wall conductivities. In all cases, the results for high aspect ratios converge to the analytical solution obtained from the TP model with the same bottom wall conductivity. However, the influence of insulating side walls remains significant even at large aspect ratios when the bottom wall is conducting. Unexpectedly, in such cases, the change in the induced magnetic field due to the presence of side walls has a dramatic effect on the velocity profile, leading to a reduced pressure gradient compared to that predicted by the TP model.
