ABSTRACT. Low, ideally zero, numerical dissipation is essential for high-fidelity, scale-resolving flow simulations, as it disrupts the accurate representation of inviscid kinetic energy and entropy conservation, both fundamental to compressible turbulence. In practice, however, numerical schemes for compressible flows rely heavily on dissipation to ensure stability, creating an inherent tension between accuracy and robustness. To address this challenge, we develop a physics-compatible numerical framework that preserves both kinetic energy and entropy at the discrete level. In this talk, we present the construction of robust, non-dissipative kinetic energy and entropy preserving (KEEP) schemes, along with their application to high-Reynolds-number aircraft flows and their extension to compressible multiphase flows.

ABSTRACT. Machine learning surrogate models can produce predictions orders of magnitude faster than classical Computational Fluid Dynamics (CFD) solvers, which makes them attractive for industrial applications where many design iterations are needed. Graph Neural Networks (GNNs) are particularly well suited for CFD since they naturally handle unstructured meshes. However, their local message passing operation leads to under-reaching, limiting their ability to model the long-range interactions common in fluid dynamics. To compensate, GNNs typically require good initial conditions, incurring extra computational cost to generate them. Transformer architectures offer an alternative: their attention mechanisms enable global information exchange, but at the cost of quadratic complexity, and they lack the geometric inductive bias that GNNs possess.

We investigate hybrid architectures incorporating attention layers into GNNs to leverage the strengths of both approaches: the geometric structure awareness of GNNs and the global modeling capability of attention for long-range dependencies. We explore different configurations and demonstrate improved accuracy over GNNs alone through numerical experiments on practical datasets.

ABSTRACT. This paper presents an optimized low-storage variational reconstruction for high-order finite volume methods on unstructured grids for compressible flow simulations. The reconstruction polynomial in each cell is obtained by minimizing a cost function defined through interfacial jumps. Unlike conventional formulations based on surface integrals of jump terms, the proposed method evaluates the weighted sum of the jumps of the reconstruction polynomial and its spatial derivatives only at interface centers. This enables efficient matrix-free evaluation of the right-hand side of the local linear reconstruction equation system, avoiding the storage of large matrices and substantially reducing memory consumption. The derivative weights in the reconstruction are optimized using an ensemble Kalman inversion framework. The optimization maximizes an efficiency metric defined as the reciprocal of the product of solution error and computational cost, evaluated using a linear advection problem with a broadband wave-packet initial condition. This design promotes both high resolution and computational efficiency. Numerical results for inviscid flow problems show that the optimized low-storage variational finite volume scheme is remarkably more efficient than both reference and non-optimized schemes, demonstrating the effectiveness of the proposed optimization strategy.

ABSTRACT. Accelerating computational simulations remains an active area of research, pursued from multiple directions across the CFD and numerical analysis communities. In all cases, the ability to generate results more rapidly directly strengthens the role of simulations in time-critical engineering workflows. One promising approach focuses on reducing the cost of core components of the computation, particularly the repeated solutions of nonlinear algebraic systems whose linearizations yield large sparse linear systems. Computational efficiency can be improved by initializing both the nonlinear and linear solver iterations with estimates that are closer to the eventual solution, thereby reducing the number of iterations required for convergence and shortening the overall time to solution. In the present work, this acceleration is achieved by leveraging data already available in memory during a simulation to train a neural network capable of predicting updates for the nonlinear and linear solvers. These predictions are then used to initialize the solvers. A critical requirement, however, is that the machine learning components themselves must execute in less time than the savings they provide. This balance can be maintained by training a lightweight neural network to achieve sufficient, but not necessarily high, accuracy. Such models are smaller, faster to train, and still effective at improving solver convergence. The approach demonstrated here employs in-situ graph neural network training and inference to initialize both nonlinear and linear solvers. Results from one-dimensional and two-dimensional test cases show improved convergence behavior and, importantly, an overall reduction in total simulation time. The effort will expanded to investigate production-oriented solvers for the final paper.

ABSTRACT. We improve the performance of Physics-Informed Neural Networks (PINNs) for modeling discontinuous solutions in hydrodynamics and relativistic hydrodynamics. Standard PINNs often struggle with convergence and accuracy near shocks. Building on Locally-Linearized PINNs (LLPINNs), which modify the Jacobian in regions of strong compression, we remove the need for a priori shock speeds by dynamically computing them from neighboring states and enforcing entropy conditions. We also introduce Locally-Roe PINNs (LRPINNs), which employ a Roe-averaged Jacobian to enhance conservation and shock resolution. The approach is extended to two-dimensional Riemann problems using divergence-based shock detection and dimensional splitting. Compared with a high-order WENO solver, our method captures sharper shock transitions, though small-scale vortical structures remain less resolved, motivating future improvements.

ABSTRACT. Introduction

Over the last decades, due to the growing availability of computational resources, high-fidelity simulations have gained increasing interest in the prediction of compressible turbulent flows over complex geometries, such as airfoils, compressor intakes, and turbine cascades. The simultaneous presence of pressure gradients and surface curvature has a strong influence on the resulting Shock Wave-Boundary Layer Interaction (SWBLI) which can be investigated on the benchmark case of the bump. Despite its geometric simplicity, the flow over a bump has been widely studied both numerically and experimentally since it reproduces many of the key physical mechanisms observed in the SWBLI in transonic flows over supercritical airfoils. Accurately capturing these complex flow features represents one of the major challenges in the numerical simulation.

Reynolds Averaged Navier-Stokes-based methods have shown difficulties in predicting the turbulent statistics and the onset of separation, while Direct Numerical Simulation (DNS) cannot be employed at high Reynolds conditions due to the prohibitively large number of grid points required. Although Wall-Resolved Large Eddy Simulations (WRLES) are more affordable than DNS, the resolution of the inner part of the boundary layer is still computationally demanding. For this reason, Wall-Modeled Large Eddy Simulation has emerged as a suitable compromise between accuracy and computational cost, solving the outer layer as a classical LES and modeling the near-wall region.

While the large majority of CFD solvers use body-fitted grids to perform accurate simulations on complex geometries, Immersed Boundary Method (IBM) based solvers, which rely on Cartesian meshes, have shown significantly improved efficiency on modern High Performance Computing (HPC) systems. However, IBM has usually been applied in low Reynolds conditions as a result of the difficulty in clustering points in the near-wall region.

In this work, IBM-based WMLES are applied to investigate the transonic shock boundary layer interaction on a circular bump. Different simulations are performed to evaluate the capability of the current methodology in capturing the main flow physics of the phenomenon. By comparing results with experimental data, different grid resolutions are employed to effectively assess the accuracy of the approach. Finally, an analysis of the SWBLI, pressure distributions, separation, and turbulent statistics is carried out in order to identify the key physical mechanisms governing the flow at a high Reynolds number.
 

Methodology

The solver used in the present study is URANOS (Unsteady Robust All-around Navier-StOkes Solver), which integrates a dynamic wall model approach with IBM, providing an efficient and scalable framework for GPU-based architecture. The solved governing equations are the compressible Navier-Stokes equations written in conservative form, closed through the ideal gas equation of state and the constitutive expression for the internal energy.

The subgrid scale viscosity model adopted is the wall-adaptive large-eddy (WALE), while spatial discretization is performed on a structured Cartesian grid using high-order finite-difference schemes. Convective fluxes are evaluated with shock-capturing formulation, while viscous terms are discretized using central schemes. Time advancement is carried out via an explicit multi-stage Runge–Kutta method where the time-step is dynamically varied to satisfy the Courant–Friedrichs–Lewy (CFL) stability condition. In order to avoid prohibitively high computational cost related to the resolution of the inner layer, an equilibrium-based wall-stress model is used. Convective and pressure terms are assumed to be zero near the wall, so the mean wall-parallel velocity and mean wall temperature distribution are determined from the resulting ODE system. This system of equations is solved on an independent grid starting from the wall, where velocity and temperature are set by boundary conditions, to an outer interface where the value of these variables is prescribed by LES. Once the ODE system is solved, the mean wall-parallel velocity and mean wall temperature distribution are used to evaluate the wall shear stress and wall heat flux.

Regarding the surface treatment in IBM, the Ghost Point Forcing Method is employed to represent the effects of the embedded geometry. First, a properly defined mapping function classifies each grid point as fluid, ghost, or solid. For each ghost point, a corresponding image point is built in the direction normal to the body boundaries, where the flow field variables are interpolated from surrounding points. Subsequently, this value is used to enforce the wall condition by updating the value at the associated ghost point. Then, the image point can be used as the initial location of the interface between LES and wall-model allowing the simulation at high Re numbers. Further details on the algorithm adopted in the implementation of the dynamic wall-model with the IBM are shown in previous works.

Results and Discussion

Simulations were performed on a circular arc bump and then compared the results with a previous experimental study at Mach number M=0.72 and Reynolds number Re=1.6∙10⁶ based on the bump length.

The thickness of the boundary layer at the inflow is used as reference length, so computational domain is a tridimensional box with size L_xxL_yxL_z= (240x65.5x10). On the y direction the domain is extended to include the IBM block which occupies the lower part of the domain and places a bump in the center of the lower surface. The height of the bump is h = 5, and the radius is equal to 163 , with a total length of 80. The length of the up- and downstream flat plate is equal to 80. Preliminary studies highlighted the importance of an appropriate resolution in the shock and separation region. For this reason, the grid is locally refined in the near wall region spanning over the entire bump height as well as near the trailing edge in the streamwise direction. In the vertical direction, a uniform fine grid is used near the wall, transitioning to a coarser distribution in the outer region. Grid points are distributed uniformly in the spanwise direction. This choice is motivated by the necessity to ensure the correct resolution across the bump, since it is not possible to refine in normal-to-wall direction when the IBM is adopted.

The proposed IBM-WMLES approach captures the main physical phenomena of the transonic flow over a circular bump. Flow accelerates from the leading edge of the bump due to the favorable pressure gradient originating a supersonic region which terminates in the rear part of the bump with a shock wave. Close to the wall, the shock presents -shaped structure with the left-leg at x ≈ 20, whereas further from the wall it becomes nearly normal to the incoming flow direction. The pressure rise associated with the surface curvature and the interaction with the shock, determine the separation and downstream reattachment of the boundary layer. Despite the simplicity of the wall-model adopted, the analysis of mean wall pressure reveals that the shock position is predicted reasonably well in comparison with experimental data. However, the major discrepancy with respect to the reference case is in the post-shock region where equilibrium assumption of the wall-model is not satisfied.

After an initial validation with a coarse grid of approximately 37 million grid points, the accuracy and the robustness of the current approach is investigated performing simulations with two finer resolutions. Grid resolution is increased in the three dimensions doubling and triplicating the total number of points. An analysis of the turbulent statistics, pressure distribution and skin friction coefficient indicates that wall model guarantees a proper resolution of the flow field with affordable computational cost.

ABSTRACT. Introduction Aerodynamic decelerators are essential components for ensuring feasible and controlled atmospheric entry trajectories of space vehicles. Parachutes, in particular, exhibit complex dynamic behaviors during deployment and inflation, which have a strong impact on the stability of the vehicle. Despite past and ongoing efforts, correlating the performance of the parachute with the operating conditions remains a significant challenge [1]. Therefore, research efforts are still necessary to improve the understanding of parachute dynamics in light of the evolving characteristics of modern spacecraft and the wide range of atmospheric conditions associated with mission destinations. The objective of this study is to validate the recently developed numerical methods for simulating three dimensional bodies in compressible flow through benchmark cases.

Methods The computational framework is based on STREAmS 2, a Cartesian finite-difference solver for direct numerical simulation of high-speed turbulent flows [2]. The spatial discretization of the convective terms in the Navier–Stokes equations relies on a combination of a robust energy-preserving scheme and a shock-capturing scheme. The flow solver is advanced in time using a low-storage third-order Runge-Kutta scheme.

The immersed body is introduced into the framework as a triangulated (and possibly open) surface. Boundary conditions on the body are imposed at each time step using a sharp-interface Immersed Boundary (IB) method, specifically designed for thin-walled open structures. The discretization scheme is properly modified in the vicinity of the body to prevent any finite-difference stencil from containing nodes from both sides of the immersed interface. The solver relies on a multi-task implementation using Open MPI, and it is accelerated on GPUs through CUDA.

Results Figure 1 shows the pressure coefficient for a supersonic flow past circular cylinder (Ma=2, Re=300). The results agree well with those obtained in [3]. Further comparisons on both 2D and 3D test cases will be presented.

References [1] Luca Placco, Giulio Soldati, Matteo Bernardini, and Francesco Picano. On flight instabilities of capsule-rigid parachute system during supersonic planetary descent. Aerospace Science and Tech-nology, 160:110026, 2025. [2] Matteo Bernardini, Davide Modesti, Francesco Salvadore, and Sergio Pirozzoli. Streams: A high-fidelity accelerated solver for direct numerical simulation of compressible turbulent flows. Computer Physics Communications, 263:107906, 2021. [3] Daniel Canuto and Kunihiko Taira. Two-dimensional compressible viscous flow around a circular cylinder. Journal of Fluid Mechanics, 785:349–371, 2015.

ABSTRACT. To implement the k-omega shear stress transport (SST) model with the immersed boundary method, the non-plausible behaviors of the transport variables in the near-wall region are challenging. To avoid numerical problems, we adopted the inverse of omega as the transport variable and derived the model equation from the original omega equation. Additionally, we introduced a wall-damping function so that the near-wall solutions of the transport variables are linear or constant with respect to the wall distance. The derived model was implemented in the Cartesian-grid-based flow solver UTCart and validated through flow simulations around a NACA 0012 airfoil. The results of the validation problem showed good agreement with the reference results using the original k-omega SST model with body-conforming grids. The presentation will include other validation problems, e.g., the flow simulations around the full aircraft configuration, and comparisons to the results using SA turbulence model.

ABSTRACT. One of the most important and time consuming parts of a Computational Fluid Dynamics workflow is the generation of the computational mesh. Cartesian cut-cell methods circumvent this problem by automatically generating meshes without user intervention. Traditionally, cut-cell methods have struggled to provide solutions to the Reynolds-Averaged Navier-Stokes equations for high Reynolds number flows. Firstly, the irregular cell sizes near no-slip walls due to cut cells introduce noise into gradient reconstruction in the boundary layer. Secondly, satisfying the off-wall meshing requirements for wall resolved flows places too many isotropic cells in the wall-transverse direction; incurring significant cost compared to boundary-conforming methods. In this work, we compare two interior forcing point methods: the Spalart-Allmaras wall function and an ordinary differential equation (ODE) wall model. We use three NASA turbulence modeling resource verification cases: the 2D turbulent bump, the NACA 0012 airfoil, and the high-lift multielement airfoil. The ODE wall model demonstrates significant improvement over the SA wall function and comparable error convergence rates to state of the art boundary conforming methods.

ABSTRACT. We develop Global Stability Analysis (GSA) of multi-step integration methods to analyze the occurrence of physical and spurious modes arising in spectrally discretized convection-diffusion dynamics. Stability conditions are obtained and further quantified concerning the presence of both physical and spurious modes, as well as physical modes only. The Adams-Bashforth and Adams-Moulton methods provide specific examples for which GSA yields explicit conditions on the Peclet and CFL numbers that need to be respected to achieve simulations free of spurious modes.

ABSTRACT. Convection-diffusion equations constitute the mathematical foundation of many governing laws in Computational Fluid Dynamics (CFD), including the Navier-Stokes and Reynolds-Averaged Navier-Stokes (RANS) equations. To numerically resolve these equations, High-Order (HO) schemes have gained significant popularity in recent years due to their computational efficiency in resolving smooth regions of the solution.

Despite these advantages, achieving robust HO discretizations remains challenging. While stabilization for hyperbolic operators is well established, robust treatment of diffusion—especially for strongly anisotropic diffusion on distorted meshes—remains an active research area. Moreover, the parabolic nature of the problem imposes severe time-step restrictions on explicit schemes due to the von Neumann stability condition, making fully implicit formulations essential.

In this work, we present a theoretical analysis of a fully discrete, time-implicit DGSEM based on the Bassi and Rebay 2 (BR2) formulation, focusing on domain invariance and entropy stability. We derive robustness conditions and introduce an algebraic graph viscosity for stabilization. A Flux-Corrected Transport (FCT) limiter is used to couple a robust low-order scheme with the high-order operator, preserving accuracy while enhancing stability. These properties are extended to generic distorted meshes, and numerical benchmark tests are provided to validate the theoretical results.

ABSTRACT. This paper constructs a multistep, multistage continuous-explicit Runge–Kutta (MSCERK) temporal predictor and applies it within a time-accurate local time-stepping discontinuous Galerkin (LTS-DG) framework. The motivation comes from temporally multiscale unsteady flows, where a global explicit time step is often dictated by the most restrictive local time scale and therefore leads to redundant updates in the rest of the domain. In the LTS-DG predictor–corrector framework, each element requires a continuous-in-time local predictor so that interface fluxes can be integrated consistently from neighboring states at asynchronous time levels. The constructed MSCERK predictor uses both current-step and previous-step information to build this continuous representation. In particular, the two-step MSCERK construction reduces the number of currentstep predictor stages required for a prescribed continuous order, making it well suited for LTS-DG calculations in which stage residual evaluations dominate the cost. The order conditions of the two-step continuous extension are summarized, and representative numerical tests are reported to demonstrate accuracy and efficiency.

ABSTRACT. Developments toward entropy-stable high-order methods have brought high-fidelity computational fluid dynamics simulations closer to industrial application. However, the best choice of temporal integration method remains unclear. This presentation will compare three temporal integration schemes. The first is an explicit Runge-Kutta solver using error-adaptive time step size selection. Second, we use diagonally implicit Runge-Kutta, tackling the implicit solve using a Newton-Krylov solver with the exact Jacobian assembled via automatic differentiation. Finally, we use an entropy-stable space-time method, where space and time are solved in a fully coupled fashion. The three methods are compared on an Euler test case that mimics the challenges of wall-bounded flows using an anisotropic grid. On the selected test case, implicit Runge-Kutta is about six times faster than explicit Runge-Kutta, suggesting that implicit methods equipped with modern solvers can be highly competitive for computational aerodynamics.

ABSTRACT. 1. Introduction The predictive simulation of reacting flows is essential for designing next-generation aerospace propulsion systems, including scramjets, rotating detonation engines, and hypersonic propulsion concepts. These applications involve strong coupling between compressible fluid dynamics, heat release, molecular transport, and finite-rate chemical kinetics, leading to thin reaction zones, sharp thermochemical gradients, shock waves and multiscale dynamics. Accurate numerical resolution is critical for predicting flame dynamics under extreme operating conditions. Yet, conventional low-order discretization methods often suffer from excessive numerical dissipation and dispersion that limiting predictive accuracy for predicting the relevant combustion regimes.

High-order discretization methods offer a promising pathway to overcome these limitations by providing superior accuracy per degree of freedom and reduced numerical dissipation, enabling improved resolution of steep gradients and fine-scale thermochemical structures. In particular, p-refinement achieves exponential convergence and lower cost compared to h-refinement, driving increasing adoption of high-order approaches. Among these, the Spectral Element Method (SEM) combines high accuracy with computational efficiency. Its element-local formulation yields favorable properties, including diagonal mass matrices and high arithmetic intensity, making it well suited for accelerated architectures based on graphics processing units (GPUs).

GPUs have emerged as a key enabling technology, providing massive parallelism and high memory bandwidth well suited for compute-intensive, element-local operations and stiff chemistry integration, which typically dominates runtime. High-order spectral element methods are particularly well aligned with GPU architectures due to their localized data access and high arithmetic intensity, enabling efficient acceleration. Efficient GPU-based integration of finite-rate chemical kinetics is also essential to fully exploit modern heterogeneous architectures. Leveraging GPUs is therefore essential to achieve the performance and scalability required for predictive high-order combustion simulations. These challenges are particularly evident across combustion regimes. Deflagrations involve thin reaction zones governed by diffusive and reactive transport, demanding high-order accuracy to resolve flame structure and propagation. In contrast, detonations involve tightly coupled shock–reaction interactions, producing strong discontinuities and complex wave dynamics that impose stringent requirements on stability, shock resolution, and computational efficiency.

While several GPU-accelerated solvers exist for discontinuous Galerkin methods and incompressible flows, GPU-native high-order continuous Galerkin solvers capable of handling compressible reacting flows with finite-rate chemistry remain largely unexplored. The present work addresses this gap by extending the GPU-accelerated spectral element solver SOD2D to simulate compressible multi-component reactive flows with finite-rate chemistry. The approach combines high-order spatial discretization with entropy-viscosity shock capturing and projection-based stabilization to ensure numerical robustness while preserving accuracy. Chemical kinetics integration is fully offloaded to GPUs through a dedicated acceleration library, enabling efficient treatment of stiff reaction mechanisms. The resulting framework enables accurate and scalable simulations of reactive flows across combustion regimes, from laminar premixed flames to detonation waves, providing a foundation for predictive simulations of high-speed propulsion systems on modern GPU-accelerated supercomputers.

2. Methodology 2.1 Numerical Methods The present work employs SEM as the spatial discretization framework. A key advantage lies in its computational efficiency: unlike the sparse banded mass matrix in traditional FEM, SEM's mass matrix is purely diagonal, significantly reducing memory requirements and computational cost. This property arises from the use of Lagrange basis functions, with interpolation and integration collocated at Gauss-Lobatto-Legendre (GLL) nodes. Since the basis functions satisfy the Kronecker delta property, for the weak form of the mass matrix, only the diagonal terms with i=j yields non-zero contributions, thus leading to a diagonal mass matrix.

The implementation employs a hybrid approach combining CUDA and OpenACC directives, balancing performance with portability and maintainability. The element-local nature of SEM naturally lends itself to parallel execution: nodal operations require only data from the corresponding element in most situation, thus minimizing global memory traffic.

The chemistry integration, which dominates cost due to stiff ODE solving, is fully offloaded to GPUs through ChemInt—a newly developed in-house library for GPU-accelerated stiff chemical kinetics. All the development and performance benchmarking above are conducted on the MareNostrum 5 at Barcelona Supercomputing Center, utilizing NVIDIA H100 GPUs for multi-GPU strong scaling.

2.2 Governing Equations and Stabilization Techniques The governing equations for compressible multi-component reactive flows comprise the Navier-Stokes equations along with the species transport equations for every species s, where t denotes the physical time, \rho, \boldsymbol{u}, p, T, \boldsymbol{\tau}, \lambda, \boldsymbol{V_k}, Y_s, \boldsymbol{V_s} and \dot{\omega}_s represent the density, velocity vector, pressure, temperature, viscous stress tensor, thermal conductivity, correction velocity, and the mass fractions, diffusion velocities, and reaction rates of species s, respectively. Additional entropy viscosity (EV) terms are added on the right-hand sides of the governing equations. These artificial diffusion terms rely on an additional viscosity \mu_\mathrm{EV}, which is constructed following the entropy residual approach described in the work of -L. Guermond, et al. Therefore, these EV terms capture shocks and only be activated nearby, thus stabilizing the flow precisely. This local low-order stabilization is completed with global high-order "Local Projection Stabilization" (LPS) terms, included following the procedure described in, when computing the equations weak form integrals within the continuous Galerkin spectral element method framework. This projection-based stabilization scales with O (h^{p+1}), becoming increasingly targeted as polynomial order increases and thus minimizing interference with resolved flow features.

Finally, the system is closed by thermodynamic relations and temperature-dependent transport properties for an ideal gas mixture. To further enhance numerical stability, a flux splitting strategy is applied to the convective terms following the work of Kennedy and Gruber. The time integration of the convective and diffusive fluxes is performed using a third-order strong stability-preserving explicit Runge-Kutta schemes, while operator splitting is applied to integrate the chemical reactions rates \dot{\omega}_s.

3. Results 3.1 One-dimensional Hydrogen Laminar Premixed Flame First, a one-dimensional laminar premixed hydrogen flame is simulated for a stoichiometric H2/air mixture, at initial conditions of 300K and 1atm, and with the reduced mechanism of P. Boivin et al. involving 9 chemical species and 12 reactions. The computational domain has been discretized with third-order spectral elements of length h_e=50~\mathrm{\mu} m, while the time step has been computed from Courant-Friedrich (CFL) and Fourier (Fo) numbers set to 0.9 and 0.3, respectively. The computed flame structure and propagation speed show excellent agreement with reference solutions obtained from Cantera (seen in Figure 1), validating the accurate implementation of transport properties and chemical kinetics.

3.2 Two-dimensional Hydrogen/Oxygen/Argon Detonation Second, to evaluate the solver's performance in capturing supersonic combustion regimes, a two-dimensional H2/O2/Ar detonation at low pressure is simulated with the same reduced mechanism from P. Boivin et al. The space-time discretization is here set-up with third-order quadratic spectral elements of size h_e=75~\mathrm{\mu} m, as well as (CFL,Fo)=(0.35,0.1) conditions. The numerical simulation is initialized following the work of R. F. Johnson et al., allowing the development of a self-sustained detonation wave. The resulting cellular detonation structure (seen in the Figure 2), characterized by triple points and transverse waves, is successfully captured, demonstrating the solver's ability to simulate the detonation initiation and quenching processes, and the ability to resolve complex shock-chemistry interactions inherent in hydrogen detonations.

4. Conclusions A GPU-accelerated reactive flow solver has been developed by extending the spectral element framework SOD2D with finite-rate chemistry capabilities. Advanced stabilization techniques ensure numerical robustness for supersonic flows with strong discontinuities while preserving high-order accuracy. Validation cases demonstrate accurate resolution of both laminar flame propagation and detonation structures, establishing a promising computational framework for predictive simulations of high-speed propulsion systems on modern GPU-accelerated supercomputers.

ABSTRACT. The growing recognition of the importance of detailed reaction mechanisms in high-speed combustion flows has made their application in numerical simulations increasingly essential. To address the challenges associated with implementing these detailed mechanisms, this paper proposes a three-tier collaborative computational framework. This framework dynamically reduces the mechanism scale during computation while maintaining high-performance parallel processing. In a homogeneous reactor model, the parallel implementation of this framework achieved a 30% reduction in computational time. Furthermore, when applied to an oblique detonation flow field, the framework reduced the reaction computation time by 93%.

ABSTRACT. Initiated chemical vapour deposition (iCVD) is an increasingly important technique for the fabrication of conformal polymer thin films, offering a solvent-free process with excellent compositional control and compatibility with temperature-sensitive substrates. Despite these advantages, the iCVD reactor performance is highly sensitive to operating conditions, and process optimisation is often carried out through extensive experimental trial-and-error. In this work, a multiphysics 3-D model of a laboratory-scale iCVD reactor is developed to provide mechanistic insight into the coupled transport and reaction phenomena governing polymer film growth. The model combines laminar non-isothermal flow, multicomponent mass transport, electric potential, and detailed surface polymerisation reactions. The gas-phase initiator decomposition, precursor transport, and surface adsorption, initiation, propagation, and termination reactions are explicitly accounted for. Simulations allow to determine the spatial distributions of velocity, temperature, and chemical concentrations under low-pressure operating conditions, typical of these reactors. The radiative heat transfer from the filament is shown to dominate the thermal field, strongly influencing initiator decomposition and surface reaction rates. Predicted deposition rate profiles exhibit near-radial symmetry and high uniformity across the substrate, with minor variations of approximately 2 nm/min. A systematic comparison between simulated and experimental deposition rates over a range of pressures demonstrates good agreement, validating the modelling framework. A sensitivity analysis indicates that the surface propagation rate constant is the dominant kinetic parameter affecting deposition rates, while reactor pressure, filament voltage, cooling stage temperature, and inlet flow rates significantly influence substrate temperature. All in all, the model provides valuable insight into iCVD reactor behaviour and offers a robust tool for guiding reactor design and preliminary process optimisation, reducing reliance on the costly experimental iterations.

ABSTRACT. Flame control strategies can be broadly classified into passive and active approaches. Passive methods are simple and require no external energy input, however, they lack adaptability to changing flow conditions. Typically, they involve geometric modifications of the injection system to stabilize the flame by altering the flow field.

Active flow control methods offer greater flexibility, as they do not rely solely on optimizing the flow domain for a specific operating regime. Their primary advantage lies in their ability to dynamically adjust control parameters, such as excitation amplitude and frequency, enabling real-time responses to variations in flow conditions, for example, changes in inlet velocity or temperature.

In the present work, we employ both passive and active flow control strategies to improve the combustion characteristics of nitrogen-diluted hydrogen and ammonia–hydrogen flames, which are considered eco-friendly, carbon-free energy carriers. Injection systems incorporating cylindrical and star-shaped bluff bodies are investigated, with harmonic forcing applied to exploit their combined effects and enhance fuel–oxidizer mixing. The study is conducted using Large Eddy Simulation (LES) and an advanced numerical framework based on a two-stage procedure that involves two computational codes, each with second-order finite-volume and sixth-order compact difference schemes. The primary objective is to assess how the bluff-body shape and excitation parameters (frequency and amplitude) influence the formation of large- and small-scale vortical structures, and how the resulting flow modifications affect flame characteristics, including flame shape, fuel consumption, temperature distribution, and NOx formation.

ABSTRACT. The HYPERSONICS Mini-Symposium aims to bring together specialists in the numerical simulation of hypersonic flows in order to review the state of the art in numerical methods suitable for resolving the complex physical phenomena involved in such flows. The flight domain considered in this MS corresponds to altitudes at which the flow regime remains continuous. Under this assumption, the flow can be described within the framework of continuum mechanics using the Navier–Stokes equations, which express the conservation of mass, momentum, and energy in the fluid. Hypersonic flows are characterized by extreme physical phenomena, such as detached shock waves across which the physical variables describing the flow undergo abrupt discontinuities, including a sharp conversion of kinetic energy into internal energy that manifests as a steep rise in temperature. Moreover, in a very thin region adjacent to the wall, known as the boundary layer, strong gradients of velocity and temperature in the direction normal to the wall induce the transfer of momentum and energy to the surface. Accurate prediction of viscous drag and wall heat flux is therefore a critical issue for trajectory control and for an accurate design of the thermal protection systems of hypersonic vehicles. Given the severity of the physical phenomena and their multiscale nature, numerical methods face a significant challenge: they must simultaneously provide robustness and accuracy to capture both intense shock waves and the detailed flow structure within the boundary layer, while remaining as insensitive as possible to mesh quality. Achieving this balance between accuracy and robustness requires precise control of the numerical dissipation inherent in the discretization methods used to solve the partial differential equations governing the flow.

ABSTRACT. In recent years, there has been a renewed interest in the numerical simulation of hypersonic flows, after a few decades of sleep. Since the 90’s, significant progress has been made in the development of numerical schemes for compressible flows. The use of unstructured meshes has become routine, and numerous high-order schemes (i.e., higher than second order) are now available. In this paper, we want to present, and critically assess, a particular methodology: the so-called residual distribution schemes. We will consider the steady version of these methods. At the conference, we intend to present a fix to the scheme that avoid negative pressure and densities. We would also understand the iterative convergence issues, show results wih formal accuracy of 3, and hopefully some preliminary results for viscous problems.

ABSTRACT. This work focuses on the development of robust and accurate numerical methods for the simulation of three-dimensional hypersonic flows around complex vehicles. The flow regime under consideration corresponds to flight altitudes at which the continuum hypothesis is valid and the fluid behavior can be described by the compressible Navier–Stokes equations, expressing the conservation of mass, momentum, and energy. Hypersonic flows are hypervelocity flows,typically several times the speed of sound, and are characterized by extreme physical phenomena. These include strong detached shock waves, across which the flow variables experience abrupt discontinuities, leading to an intense conversion of kinetic energy into internal energy and, consequently, to very high temperatures. In addition, within a very thin region adjacent to the vehicle surface, known as the boundary layer, strong velocity and temperature gradients develop in the direction normal to the wall. These gradients induce significant transfers of momentum and energy toward the wall, resulting in severe aerodynamic and aerothermal loads. The accurate numerical prediction of these loads, and in particular of wall heat fluxes, remains a critical challenge for the design of thermal protection systems of hypersonic and space vehicles.

We aim at developping a robust and accurate subface-based, cell-centered finite volume method for solving the three-dimensional compressible Navier–Stokes equations in hypersonic conditions on hybrid unstructured grids. The use of unstructured meshes is essential for handling complex three-dimensional geometries, but it also imposes stringent requirements on the robustness, stability, and accuracy of the numerical discretization. A major difficulty lies in the simultaneous resolution of strong curved bow shocks ahead of the vehicle and of the thin boundary layers along its surface. While sufficient numerical dissipation is required in the inviscid fluxes to stabilize shock waves, excessive dissipation may severely degrade the prediction of viscous stresses and heat transfer within the boundary layer. Balancing these conflicting requirements is one of the central objectives of this work.

ABSTRACT. Hypersonic turbulent flows with shock–boundary-layer interactions remain challenging for Reynolds-averaged Navier–Stokes (RANS) modeling, particularly in terms of accuracy, but also of robustness, especially when high-order numerical schemes are employed to discretize the governing equations. This work presents a DNS-driven strategy to improve turbulence modeling for hypersonic applications. First, a reformulation of the k–ω SST model based on a logarithmic transformation of the conservative dissipation variable is introduced, significantly enhancing numerical stability in highly compressible flows while preserving predictive accuracy of the baseline model. Second, the model is augmented with established compressibility corrections and data-driven corrections obtained via symbolic regression from direct numerical simulation databases. The proposed approaches are assessed on hypersonic boundary layers and compression-corner configurations, demonstrating improved robustness and more accurate predictions of skin friction and heat transfer.

ABSTRACT. Deterministic Reynolds-averaged Navier–Stokes (RANS) simulations remain a standard tool for hypersonic flow analysis, yet their predictive reliability is often limited by the neglect of uncertainties in inflow and boundary conditions and in modeling parameters. Even for canonical configurations such as large cone–flare geometries, small input variations can produce significant differences in predicted aerothermal loads. This work presents a systematic and reproducible uncertainty quantification (UQ) study of Run 33 of the large cone–flare benchmark of Hoste et al., corresponding to a code validation campaign originally introduced by Holden et al., to assess the impact of aleatory uncertainties on key quantities of interest for hypersonic applications. Seven uncertain inputs are considered, including freestream Mach number, pressure and temperature, wall temperature, laminar and turbulent Prandtl numbers, and nose bluntness. Probability distributions are assigned based on experimental variability, measurement uncertainty, and manufacturing tolerances. Latin Hypercube Sampling is used to generate an ensemble of high-fidelity RANS simulations performed with the SU2 solver and the SST turbulence model. Kriging surrogate models are constructed and adaptively refined to enable efficient uncertainty propagation, statistical analysis, and global sensitivity assessment. The results reveal substantial variability in predicted aerothermal responses, highlighting the limitations of deterministic CFD. In particular, uncertainty in the turbulent Prandtl number and nose bluntness is found to dominate the variance of most outputs, underscoring the need for uncertainty-aware hypersonic CFD analyses.

ABSTRACT. The motivation of this work is to design computationally efficient methods for hypersonic stationary flows simulation for reentry vehicles. This is of interest for fast tools using low fidelity description of 3D objects, or for high fidelity large scale simulations where parallelisation efficiency and robustness are crucial. Unstructured meshes in 3D are considered, and the scheme should be robust with no condition on them. The point is to design explicit schemes, stable when using large time steps to reach quickly the stationary solution. The leading idea is to use robust non conservative schemes, but achieving a conservative form when convergence is reached. This is the case of the point implicit relaxation scheme used by Gnoffo et al for instance. In this work, an Advection Upstream Splitting Method like scheme is proposed for the Euler hyperbolic part of Navier Stokes equations. These schemes are successfully used in the hypersonic aerodynamic community and the convection and pressure terms of the flux will be treated such that they will be robust and convergent for large CFL numbers. Moreover, an isentropic approximation will be used to implicit pressure in the acoustic part of the scheme, recovering the physical entropic solution at convergence. Momentum dissipative terms and thermal diffusion in Navier Stokes equations will be treated using an operator splitting and considered as two separate hyperbolic systems, as in the Hyperbolic Navier Stokes HNS20 model developed by Nishikawa et al. This strategy is used for its capability to smooth solutions on arbitrary unstructured meshes and to accurately compute physical variables as density, temperature, velocities, and their gradients computed potentially at same order of accuracy. In particular viscous forces and thermal fluxes may have a proper behavior even with a simple scheme. An explicit unconditionally stable scheme is designed for these Navier stokes terms, as well with a point implicit relaxation strategy. First results for the Euler system and full Navier Stokes equations will be presented, but validation and further computational optimization of the code are still ongoing work.

ABSTRACT. We applied the matrix stability analysis of the carbuncle phenomenon to ideal MHD equations for the first time. Perpendicular and parallel shocks are analyzed for Roe, HLL and HLLD schemes. As a result, HLL scheme is always stable in both perpendicular and parallel shocks, while Roe and HLLD schemes are always unstable in parallel shocks, and the stability of them depends on the relationship between upstream Mach number and plasma beta parameter in perpendicular shocks. Moreover, the results suggested the existence of a critical fast magnetosonic Mach number at which the perpendicular shock becomes unstable. Furthermore, we proposed an AUSM-like expression of HLLD scheme and improved the shock robustness by adding the velocity diffusion term into the pressure flux. In particular, the numerical fast magnetosonic speed varying with plasma beta parameter was introduced for the robustness against the parallel shocks. The matrix stability analysis indicated that the present scheme, HLLD+M, improved the stability of both perpendicular and parallel shocks. We confirmed that the present scheme, compared with existing methods such as HLLD, achieved better shock robustness against the carbuncle instability through numerical tests.

ABSTRACT. Advances in high-performance computing and large-scale data acquisition have significantly expanded our capability to characterize turbulence and, more broadly, chaotic flows across a wide range of regimes. However, the intrinsic multiscale, nonlinear,and intermittent nature of these systems, together with their sensitive dependence on initial conditions, continues to pose formidable challenges for predictive modelling and control. This minisymposium focuses on emerging methodologies that integrate reduced-order modelling (ROM), machine learning (ML), data assimilation (DA), and uncertainty quantification (UQ) to develop interpretable and computationally tractable representations of turbulent and chaotic dynamics. Emphasis will be placed on physics-constrained ROM formulations, such as Galerkin and Petrov–Galerkin projections, operator-inference approaches, and nonlinear manifold-based reductions, which preserve key invariants, modal interactions, and multiscale couplings. Special consideration will be given to their applicability to chaotic attractors, transient dynamics, intermittency, and non-normal growth mechanisms. Complementary data-driven paradigms, including deep learning architectures, kernel methods, and hybrid ML–ROM closures, will be discussed in relation to their ability to capture coherent structures, intermittency, and extreme-event statistics. This minisymposium will also highlight recent developments in data assimilation for chaotic flows, including ensemble-based approaches and forecast analysis cycling methods designed to stabilize ROM and ML surrogates. These tools enable real-time state estimation, online model correction, and adaptive basis refinement in settings where sensitivity to perturbations and model inadequacy can otherwise lead to rapid forecasting performance degradation. Contributions addressing uncertainty quantification, from Bayesian inference and ensemble-based approaches to stochastic reduced models, will underline the importance of robust error characterization and propagation in predictive simulations of turbulence and chaos. A dedicated theme concerns the construction of next-generation digital twins for turbulent and chaotic flows, wherein physics-based ROMs,ML enhanced closures, and efficient DA frameworks are tightly integrated to support prediction, optimal control and design optimization. Applications span from canonical chaotic systems to complex industrial and environmental configurations. By bringing together experts in computational fluid dynamics, scientific machine learning, dynamical systems, and applied mathematics, this minisymposium aims to consolidate the state of the art and stimulate cross-disciplinary collaboration toward high-fidelity, uncertainty-aware, and real-time digital representations of turbulence and chaotic flows.

ABSTRACT. 1. Introduction Computational Fluid Dynamics (CFD) plays a central role in the analysis of complex flow systems. However, high-fidelity Reynolds-Averaged Navier–Stokes (RANS) simulations are often computationally expensive, which limits their applicability in tasks such as design optimization, uncertainty quantification, and flow control. Reduced Order Models (ROMs) address this challenge by constructing low-dimensional surrogates capable of capturing the dominant flow dynamics at a fraction of the computational cost. This work proposes two complementary, data-driven strategies that leverage scientific machine learning to enhance the accuracy and robustness of ROMs for CFD applications.

2. Methodology The first strategy [1] introduces a space-dependent aggregation framework for non-intrusiveROMs. In this approach, multiple surrogate models—built from different combinations of dimensionality reduction and regression techniques—are locally combined through spatially varying convex weights. Two aggregation pipelines are investigated: one based on aggregating RANS field predictions before ROM construction, and another combining ROMs associated with different RANS turbulence models. The aggregation weights are learned using random forest regression or neural networks, leading to improved predictive accuracy in benchmark problems such as transonic airfoil flows and the periodic hills configuration. The second strategy [2, 3] addresses parametric closure modeling for intrusive ROMs. A deep operator network is trained to learn a reduced correction term that compensates for unresolved dynamics across varying parameters, effectively enhancing the ROM closure. This learned operator enables the intrusive ROM to generalize more reliably across the parameter space. 3. Conclusions Together, these approaches advance hybrid reduced-order modeling by combining surrogate aggregation techniques with operator-learning-based closure models. The resulting frameworks provide a scalable and flexible pathway toward fast, accurate, and physically consistent CFD predictions, particularly suited for parametric studies and many-query applications. References [1] Anna Ivagnes, Niccolo Tonicello, Paola Cinnella, and Gianluigi Rozza. Enhancing non-intrusive reduced-order models with space-dependent aggregation methods. Acta Mechanica, pages 1–30, 2024. [2] Anna Ivagnes, Giovanni Stabile, and Gianluigi Rozza. Parametric intrusive reduced order models enhanced with machine learning correction terms. arXiv preprint arXiv:2406.04169, 2024. [3] Anna Ivagnes, Giovanni Stabile, and Gianluigi Rozza. Data-driven closure strategies for parametrized reduced order models via deep operator networks. arXiv preprint arXiv:2505.17305, 2025.

ABSTRACT. A direct application of standard denoising diffusion probabilistic models to fully developed 3D turbulence fails to produce statistically faithful velocity fields, even in the simplest single-snapshot setting. We identify this behavior as a limitation of standard diffusion objectives, which do not explicitly account for global physical constraints and therefore fail to confine the learned distribution to physically admissible configurations. To overcome this major difficulty, we propose a different implementation of the diffusion model that can implement constraints with arbitrary precision in all steps during the generative denoising phase. Using rotating turbulence as a challenge, we demonstrate that our Physics-Constrained Diffusion Model (PCDM) allows to generate high-resolution data for 3D turbulent fields and preserving high accuracy in all multiscale turbulent properties including anisotropic energy spectra, intermittency and extreme events [1]. This step allows to generate turbulent data with the same accuracy of Lagrangian 1D signals [2]. The robustness of the approach across distinct, non-equivalent implementations of the same physical constraints indicates convergence toward a consistent underlying physical representation. These findings point to broader implications for generative modeling of high-dimensional complex physical systems under constraints.

[1] Tianyi Li, Michele Buzzicotti, Fabio Bonaccorso, Luca Biferale. Physics-Constrained Diffusion Models for Synthesis of 3D Turbulent Data. Submitte to PRX (2026) [2] Li, T., Biferale, L., Bonaccorso, F., Scarpolini, M. A., & Buzzicotti, M. (2024). Synthetic Lagrangian turbulence by generative diffusion models. Nature Machine Intelligence, 6(4), 393-403.

ABSTRACT. Convection-dominated flows in under-resolved regimes usually need numerical stabilization to provide accurate simulations. Indeed, numerical oscillations may arise and often require numerical stabilization to improve accuracy and suppress spurious oscillations. The evolve–filter (EF) strategy is an asset in this setting. The discretized Navier–Stokes equations are \emph{evolved} through a time integration step, and, then, filtered through an elliptic operator. The filter removes high-frequency components, making the flow simulation smoother. The EF performance is highly sensitive to the choice of the filter radius, which controls the amount of numerical dissipation. Standard choices for the filter radius, related to the mesh size or to the Kolmogorov scale, may yield over-diffusive solutions, especially in turbulent regimes. To overcome this issue, a time-dependent filter radius that is adaptively tuned during the simulation can be considered. In the literature, the optimization has been carried out by means of data-driven approaches relying on high-fidelity data and lacking predictive capabilities. This talk proposes an innovative reinforcement learning-based (RL-based) EF that addresses these limitations. We propose an automatic process to learn the value of the filter radius $\delta$ based on RL with both data-driven and data-free reward formulations. In the data-driven case, the architecture is trained using reference data over a limited temporal window, corresponding to $\sim 25\%$ of the global simulation. Moreover, we push the learning capabilities by a data-free setting: in this case, the reward is built from physically motivated indicators, which enable the agent to learn the optimized $\delta$ with no access to the reference solutions. We demonstrate the effectiveness of the proposed RL-EF approach for marginally-resolved 2D simulations of the flow past a cylinder at Re = 1000, and decaying homogeneous turbulence at Re = 40 000, achieving significantly improved accuracy compared to standard EF methods.

ABSTRACT. Reinforcement learning (RL) and multi-agent RL (MARL) have recently emerged as promising strategies for active flow control. These control problems, frequently arising in applied sciences and engineering, are typically high-dimensional in state and control variables and often have parametric dependencies, making the learning of the control policies extremely challenging.

In this work, we explore the use of specialized neural network architectures suitable for RL and MARL to learn general policies for parametric control problems. In particular, we propose two hypernetwork-based architectures for RL and multi-agent RL that encode parametric dependencies into the policy weights and biases directly. Hypernetworks are special types of neural networks that learn weights and biases of another network.

Throughout our numerical experiments, spanning from the control of the flow past a cylinder to improve the lift-to-drag ratio to the reduction of the shear stress in a turbulent channel flow, we show that hypernetworks (i) enable effective active flow control, outperforming state-of-the-art RL and MARL algorithms, (ii) can efficiently deal with parametric dependencies, and (iii) require minimal hyperparameter tuning.

ABSTRACT. Inspired by the Equation-Free multiscale modelling approach, we show how the “embed–learn–lift” framework can be exploited for the construction of surrogate “normal-forms” (reduced-order models of minimal dimensionality) of high-fidelity Navier–Stokes simulations. These can in turn be used for computationally efficient and numerically accurate bifurcation and stability analysis tasks. The framework consists of four main steps. First, we apply manifold learning to discover the intrinsic dimension of the emergent latent space of the complex, high-dimensional Navier–Stokes spatio-temporal dynamics over the parameter space. Second, we construct “normal-forms”—that is, reduced-order models (ROMs)—in the sense that we capture the correct dimension for describing the emergent dynamics on the latent space using Gaussian Process Regression (GPR). Third, based on the constructed surrogate “normal-form” models, which dramatically reduce computational complexity, we exploit the toolkit of numerical bifurcation analysis theory to construct bifurcation diagrams in the latent spaces and perform systematic stability analysis. This allows, for example, the continuation of branches of limit cycles emanating from Andronov–Hopf bifurcations; such an analysis is intractable for high-dimensional systems such as the Navier–Stokes equations using standard Jacobian-based computations. Finally, by solving the pre-image problem, we reconstruct the bifurcation diagrams in the original high-dimensional space. The effectiveness of the approach is demonstrated by focusing on three suitable test-bed flow configurations: the wake flow past an infinitely long circular cylinder, the planar sudden-expansion channel flow, and the fluidic pinball up to chaotic conditions. These flow systems exhibit Andronov-Hopf, pitchfork (symmetry-breaking) and Neimark-Sacker bifurcations by increasing the Reynolds number, respectively. The proposed methodology successfully identifies the appropriate dimensionality of the latent space of the Navier-Stokes equations—and based on them constructs ROMs referred to here as "surrogate normal-forms" using Gaussian Process regression (GPR)—that enable the accurate tracing of the bifurcating solutions and their stability analysis. This includes the accurate tracing of limit cycles, their period, and their stability characteristics, as determined through the computation of Floquet multipliers.

ABSTRACT. Finite element-based high-order methods have emerged as a powerful alternative to overcome the weakness of finite volume methods (FVM), yet they remain vulnerable to shock-induced instabilities such as subcell Gibbs–Wilbraham (G–W) oscillations. To eliminate these oscillations, shock-capturing methods have been developed by combining a troubled-cell indicator with a local stabilization, such as a limiting strategy or an artificial viscosity. Nonetheless, conventional methods are not free from a fundamental trade-off between accuracy and robustness, largely due to the coupled nature of (i) where to flag troubled-cells and (ii) how to prescribe the form and strength of the local stabilization. Recently, the shock-capturing PID controller (SPID) achieved excellent performance up to the range of moderately strong shocks via the linear control theory with the data-driven gain optimization, while its robustness and convergence significantly degrade for hypersonic flows due to the strong non-linear mapping from the PID error signal to the required artificial viscosity. To extend the capability of SPID to hypersonic flows, we develop a fully data-driven shock-capturing framework comprising (i) a parameter-free, decision tree-based troubled-cell indicator and (ii) an ANN-assisted dilatation-based artificial viscosity (DB) for strong shocks. By introducing the data-driven DB model into SPID, a hybrid artificial viscosity method (DB_SPID) is proposed, which significantly improves the performance of SPID in hypersonic flows with a target-oriented flow stabilization over a wide range of discontinuity types and strengths. Through numerical experiments for hypersonic inviscid and viscous flows, it is demonstrated that the proposed DB_SPID method significantly improves robustness and convergence characteristics in high Mach number flows.

ABSTRACT. Numerical simulations of hypersonic flows are challenging due to the combination of thermophysical modeling and numerical strategy. Specifically, in the context of atmospheric vehicle entry, the solver has to deal with strong shock waves which provoke energy excitation and molecular dissociation. This leads to a tremendous heat transfer on the vehicle surface, mainly due to conduction. However, in the presence of gas-surface interactions, mass diffusion becomes a relevant contribution to the total heat flux due to catalysis and ablation, the latter rising from the interaction between the chemical species in the boundary layer and the material composing the Thermal Protection System (TPS). Nitridation and/or oxidation processes may take place, pushing the gas far from the wall (blowing effect) and mitigating the heat flux on the surface by material decomposition and mass loss. In some case, the surface recesses with a finite velocity, provoking a shape change of the TPS.

Most of the numerical solvers in the scientific community employ body-fitted (BF) grids, which allow the use of large cell aspect ratios, thus ensuring very high wall normal resolution in the boundary layer. However, when dealing with complex/moving geometries, multiblock structured grids or unstructured grids are necessary. As an alternative, immersed-boundary (IB) techniques exploit Cartesian meshes for the fluid domain, where the geometry is embedded as a Lagrangian (solid) domain. The governing equations are then solved in the Eulerian (fluid) domain with a forcing term which accounts for the forces exchanged between the fluid and the body, such that the equations obey to the desired boundary condition on the surface (e.g., no-slip). Specifically, thanks to a ray-tracing procedure, the cells close to the surface are tag as interface cells: in these cells, an appropriate reconstruction of the flow field variables is needed to satisfy the wall boundary condition.

This work presents numerical results for hypersonic flows obtained with two different IB techniques, compared to reference results to provide an assessment about the accuracy of IB techniques under different conditions.

ABSTRACT. Aleatoric uncertainties are propagated through a computational model representative of a reference experimental hypersonic flow, which exhibits strong nonequilibrium thermochemical effects. Namely, a hypersonic stream over a double-cone axisymmetric geometry for which observations are available. The experiment is reproduced numerically. The inflow uncertainty is first characterized, and then sampled to generate a database. The set of realizations is then postprocessed to build standard Polynomial Chaos Expansions to perform the uncertainty analysis. A global sensitvity analysis is carried out through Sobol's method. Then, the true experimental operating conditions are inferred based on the available observations through a calibration procedure.

ABSTRACT. 1. Introduction The extreme temperatures generated near the surface of a vehicle moving at hypersonic speeds give rise to complex high-enthalpy processes, including vibrational excitation, chemical reactions, thermal relaxation and gas-surface interactions. In fully-developed high Mach number turbulent boundary layers, the coupling of velocity and temperature fields, combined with shock wave interactions and finite-rate chemistry effects, changes turbulence structure and transport properties [1]. For this reason, there is a direct impact on the prediction of wall quantities that are important for the design of a Thermal Protection System (TPS) and for reliable assessment of aerodynamic heating in hypersonic applications. In this context, Direct Numerical Simulation (DNS) represents the most accurate way to analyze turbulent flows, as it resolves all turbulent scales without relying on turbulence models. However, its computational cost restricts applications to canonical configurations at moderate Reynolds numbers, such as turbulent channel flows.In recent years, extensive DNS and WRLES databases for supersonic and hypersonic boundary layers have been generated, spanning a wide range of Mach numbers, wall temperature ratios and thermochemical conditions [2,3]. These databases have proved to be of fundamental importance because they provide key reference parameters for the development of wall models for LES. Recent developments have extended ODE-based wall models to high-enthalpy turbulent boundary layers, solving simplified momentum and energy equations in ordinary differential form and incorporating enthalpy–velocity relations and generalized Reynolds analogy [4]. These approaches have shown promising results for hypersonic WMLES configurations under thermochemical non-equilibrium conditions. The objective of the present study is twofold. First, to generate a high-fidelity DNS database of compressible turbulent channel flows at moderate Mach numbers, both in perfect- and reactive-gas conditions, which will serve as a consistent reference framework for hypersonic wall-model development. Second, to perform a priori and a posteriori studies of wall-model assumptions based on the generated DNS database. DNS-resolved near-wall quantities are used to directly test the equilibrium hypotheses underlying reduced wall-model formulations. The same wall-model equations are then integrated within the numerical framework and their performance is evaluated in a fully coupled simulation environment.

2. Methodology The governing equations correspond to the compressible multi-species Navier–Stokes system with finite-rate chemistry and vibrational relaxation. In conservative form, the formulation includes transport equations for each chemical species, for momentum, for total energy and for vibro-electonic energy. The species equations describe the temporal variation of the partial densities as a balance between convection, molecular diffusion and chemical production rates. The momentum equation accounts for convective transport, pressure forces and viscous stresses. The total energy equation accounts for convective and pressure work terms, viscous dissipation and heat fluxes. An additional equation governs the vibro-electronic energy mode, including its transport and source terms due to thermochemical relaxation. Spatial discretization is performed using a high-order finite-difference formulation in generalized curvilinear coordinates with one-dimensional wall-normal stretching. The convective terms are discretized in a sixth-order skew-symmetric formulation in order to preserve the discrete kinetic energy balance and reduce aliasing errors in turbulent regimes. In the presence of strong gradients, a hybrid strategy is adopted in which the skew-symmetric central discretization is locally replaced by a fifth-order WENO reconstruction to ensure robustness and shock-capturing capability. Viscous terms are discretized using sixth-order central differences. Time integration relies on a third-order accurate, explicit Runge–Kutta scheme combined with operator splitting for stiff thermochemical source terms [5]. The solver used, DAMISO, is implemented within a multi-GPU MPI-CUDA framework and has been previously verified on canonical compressible benchmarks.

3. DNS database generation The baseline DNS database includes turbulent channel flows under both low- and high-enthalpy conditions. All simulations are performed in a canonical plane channel configuration with periodic boundary conditions imposed in the streamwise and spanwise directions and no-slip isothermal walls in the wall-normal direction. A one-dimensional grid stretching is applied toward both walls to achieve a target near-wall resolution of y+ ≈ 0.8. In order to speed up convergence, simulations are first initialized on a coarser grid and subsequently interpolated onto a finer one with final DNS resolution. Time averages are collected when the bulk temperature and wall density in the domain reach a statistically steady state. Two perfect-gas configurations are considered at Mach 1.5 and Mach 3, corresponding to well-established reference results by Modesti & Pirozzoli [6]. The computational domain extends over Lx × Ly × Lz = 4πh × 2h × (4/3)πh, where h denotes the channel half-height. For the Mach 1.5 case (PG-M1.5), the friction Reynolds number is Reτ = 219.7 and the friction Mach number is Maτ = 0.079, with a heat-flux parameter Bq = 0.048. The grid resolution is 256 × 128 × 128 in the streamwise, wall-normal and spanwise directions, respectively. The near-wall resolution is y+wall ≈ 0.8, with streamwise and spanwise spacings of Δx+ ≈ 10.8 and Δz+ ≈ 7.2. For the Mach 3 case (PG-M3), the friction Reynolds number is Reτ = 448.0 and the friction Mach number is Maτ = 0.11, with Bq = 0.14. The grid resolution is 570 × 260 × 250, maintaining y+wall ≈ 0.8, with Δx+ ≈ 10.0 and Δz+ ≈ 7.6. Instantaneous fields of density-gradient magnitude for the PG-M3 case show intense small-scale structures and elongated near-wall streaks, highlighting the presence of fully developed compressible turbulence. Reynolds-averaged wall-normal profiles for velocity, temperature, viscosity and turbulence intensities display excellent agreement with literature data [6], further verifying the numerical setup. The DNS database will be enriched with two high-enthalpy cases, at Mach 3 and with wall temperature values of 1750 K and 3500 K, respectively. For these cases, the flow is initialized with constant pressure, temperature and air composition (i.e., XN2 = 0.79 and XO2 = 0.21).

4. A-priori and A-posteriori assessment of wall modeling The second objective of the present study is the assessment of wall-model assumptions for hypersonic applications. A one-dimensional equilibrium ODE-based wall model is being implemented within the numerical framework to provide a consistent bridge between the DNS database and future WMLES applications in hypersonic regimes. Under local equilibrium and constant total-shear-stress assumptions, the reduced formulation consists of a set of ordinary differential equations describing the wall-normal evolution of the wall-model velocity, the mixture enthalpy and the species mass fractions . The momentum equation reduces to the statement that the wall-normal derivative of the total shear stress, given by the sum of molecular and turbulent contributions, is zero, implying constant total shear stress across the modeled layer. The energy equation accounts for the balance between the work of viscous stresses and the combined molecular and turbulent heat fluxes. The species equations describe the balance between molecular and turbulent diffusion of each species and the corresponding chemical source terms. The mixture thermodynamic properties are evaluated from the local species composition and temperature. In particular, the mixture enthalpy and specific heat at constant pressure are obtained as mass-fraction-weighted sums of the corresponding species properties. Turbulent heat and mass diffusivities are modeled through eddy-diffusivity relations, where the turbulent thermal conductivity is proportional to the turbulent viscosity and inversely proportional to the turbulent Prandtl number, and the turbulent mass diffusivity is proportional to the turbulent viscosity and inversely proportional to the turbulent Schmidt number. The turbulent viscosity μt is obtained using the mixing length model corrected with the Van Driest damping function, while the turbulent Prandtl and Schmidt numbers will be tested either as constants (Prt = 0.9, Sct = 0.7) or through a formulation based on the generalized Reynolds analogy [4]. In the a priori analysis, the DNS database is used to directly extract near-wall gradients, total shear stress distributions and thermodynamic property variations. These DNS-resolved quantities are introduced into the reduced one-dimensional equilibrium wall-model formulation in order to evaluate its assumptions under both perfect-gas and reacting conditions. In this context, the validity of the constant total-stress hypothesis and equilibrium turbulence closure is tested by comparing the wall-model prediction against the exact DNS wall shear stress. In contrast, the a posteriori analysis consists of embedding the wall-model formulation within the simulation framework and evaluating its performance in a fully coupled numerical environment, without directly imposing DNS values. This allows assessment of the model behavior under realistic WMLES-like conditions, whereby the WALE model is employed to compute the subgrid-scale viscosity. The influence of the wall-model exchange location is assessed by varying the matching position within 10 ≤ y+ ≤ 50. For the reacting multi-species Mach 3 cases, turbulence statistics, species mass-fraction distributions and dissociation degree are analyzed. The dissociation degree exhibits a quadratic variation across the channel core, consistent with recent high-enthalpy DNS studies. Compared to perfect-gas conditions, thermochemical non-equilibrium modifies viscosity distributions, temperature gradients and near-wall turbulence structure, directly affecting wall-model assumptions.

5. Conclusions A consistent DNS-based framework for hypersonic wall-model development has been presented. A perfect-gas database of compressible turbulent channel flows has been generated at Mach 1.5 and Mach 3 and both configurations have shown excellent agreement with available reference DNS data, providing verification of the numerical methodology. The extension of the database to multi-species reacting Mach 3 cases is ongoing and will allow a systematic evaluation of thermochemical non-equilibrium effects on near-wall modeling. Based on this DNS reference framework, systematic a priori and a posteriori studies of wall-model assumptions have been performed. The a priori assessment quantifies the validity of equilibrium and constant-stress hypotheses using DNS-resolved near-wall fields, while the a posteriori assessment evaluates the same wall-model formulation when embedded in a fully coupled numerical environment. The present work establishes a coherent link between DNS database generation and physics-based wall-model assessment for future hypersonic Wall-Modeled Large-Eddy Simulations.

References [1] F. Spina, A. J. Smits, and S. K. Robinson. The physics of supersonic turbulent boundary layers. Annual Review of Fluid Mechanics, 1994. https://doi.org/10.1146/annurev.fl.26.010194.001443.

[2] M. Cogo, F. Salvadore, F. Picano, and M. Bernardini. Direct numerical simulation of supersonic and hypersonic turbulent boundary layers at moderate-high reynolds numbers and isothermal wall condition. Journal of Fluid Mechanics, 2022. https://doi.org/10.1017/jfm.2022.574.

[3] D. Passiatore, L. Sciacovelli, P. Cinnella, and G. Pascazio. Thermochemical non-equilibrium effects in turbulent hypersonic boundary layers. Journal of Fluid Mechanics, 2022. https://doi.org/10.1017/jfm.2022.283.

[4] H. Huang, H. Su, Q. Guo, P. Liu, X.Yuan, and J. Cai. Ordinary-differential-equation-based wall model for large eddy simulation of hypersonic high enthalpy turbulent boundary layer. Physics of Fluids, 2025. https://doi.org/10.1063/5.0295932.

[5] G. Pascazio, D. Ninni, F. Bonelli, and G. Colonna. Hypersonic flows with detailed state-to-state kinetics using a GPU cluster. In Plasma Modeling (Second Edition): Methods and applications. 2022. https://doi.org/10.1088/978-0-7503-3559-1ch10.

[6] D. Modesti and S. Pirozzoli. Reynolds and mach number effects in compressible turbulent channel flow. International Journal of Heat and Fluid Flow, 2016. https://doi.org/10.1016/j.ijheatfluidflow.2016.01.007.

ABSTRACT. We analyze the excitation of second Mack-mode disturbances in the boundary layer over the HIFLIER flight vehicle during its ascent. The various sources of disturbances considered relevant to free flight, including free stream turbulence and particulates, are analyzed within the CHAMPS CFD framework. A combination of high-fidelity simulation capabilities and reduced-order modeling are employed to determine which sources were most likely responsible for initiating transition to turbulence over the straight cone flying at approximately Mach 6.

ABSTRACT. Multicomponent real-fluid flows are fundamental to a wide range of engineering and natural systems, such as high-speed vehicles, energy conversion devices, chemical and process engineering, and extreme weather processes. Accurate representation of real fluid flows is critical for advancing predictive simulation tools that support safer, more efficient, and more sustainable engineering designs. This work develops three-dimensional (3D) high-order numerical methods for multicomponent flow simulation with a real-fluid equation of state (EOS) representation. The diffuse interface concept is adopted to approximate flow discontinuities, such as shock waves and material interfaces. As a result, the localized artificial viscosity method is well suited for flow discontinuity capturing as it can spread a discontinuity over a finite thickness, which is controlled by the diffusivity strength. A conservative thermodynamic calculation approach is used to evaluate thermal states of real fluids, which can effectively eliminate the so-called spurious pressure oscillation, even for real fluids with highly nonlinear EOS, although it has to increase the computational cost spent on thermal state evaluation. The real-fluid effects on the Richtmyer–Meshkov instability (RMI) are then studied with high-order numerics. It finds that high-order schemes significantly outperform their low-order counterparts.

ABSTRACT. Numerical discretization of spatial derivatives in partial differential equations inherently introduces errors due to the replacement of the continuous model by its discrete counterpart. At the operator level, discrete derivative approximations can be interpreted as exact derivatives acting on implicitly filtered representations of the solution. Based on this concept of numerically induced filtering, the Approximate Deconvolution Discretization (ADD) method was previously proposed to compensate for discretization-induced filtering effects. However, the original ADD formulation was restricted to periodic boundary conditions.

In this work, we extend the ADD framework to problems with non-periodic boundary conditions. The induced filter is constructed locally using Lagrange interpolation on asymmetric stencils near the boundaries, while maintaining the formal order of accuracy. The ADD procedure is applied to the first-derivative operator, whereas the second derivative is discretized using a stable high-order compact scheme in order to isolate the stability properties of the modified first-derivative operator.

The proposed formulation is assessed using the two-dimensional vorticity–stream function formulation of the Navier–Stokes equations and the regularized shear-driven cavity flow (Burggraf flow), for which an analytical solution is available. Stability is analyzed through eigenvalue spectra of the modified operators. Although direct inversion of the Lagrange-based induced filter leads to operators that are theoretically unstable for certain stencil sizes, numerical simulations demonstrate that stable and highly accurate solutions can still be obtained for low and moderate Reynolds numbers, and for higher Reynolds numbers on sufficiently fine grids. In these cases, viscous dissipation suppresses the growth of unstable numerical modes, leading to conditionally stable behavior despite operator-level instability.

The results show significant accuracy improvements compared to standard second-order finite difference discretization. Ongoing work focuses on identifying stability limits of the approach and on developing strategies to mitigate near-boundary instabilities.

ABSTRACT. High-order numerical methods offer greater accuracy and efficiency than traditional low-order approaches, but their practical adoption can be limited by stability and meshing challenges. Ensuring nonlinear stability is key to making these methods robust and reliable for real-world engineering problems, especially in situations involving moving meshes—such as fluid-structure interactions, moving parts, and adaptive grids. The compressible Navier–Stokes (CNS) equations possess a convex entropy function, related to thermodynamic entropy. This function acts as an admissibility criterion for identifying physically relevant weak solutions and provides a means of demonstrating nonlinear (entropy) stability by bounding in the L2 norm the solution. Maintaining entropy stability for CNS equations has been a major focus of research. Recent advances show entropy stability can be preserved even as computational grids deform, using a discontinuous Galerkin split form of the arbitrary Lagrangian Eulerian (ALE) approach. This work extends these results, proving entropy stability for moving wall boundary conditions and demonstrating the effectiveness of high-order entropy-stable DG methods in large-scale simulations with mesh motion.

ABSTRACT. Adaptive Mesh Refinement (AMR) dynamically adjusts mesh resolution to enhance simulation accuracy and efficiency in regions with shock waves and material interfaces. While most AMR algorithms are developed within the Eulerian framework, they suffer from numerical diffusion, leading to interface smearing and reduced accuracy. Lagrangian methods can directly track material motion and preserve sharp interfaces through adaptive grid deformation. However, most existing Lagrangian-AMR methods are based on non-conservative schemes, and research on conservative schemes remains limited despite the Lax-Wendroff theorem’s emphasis on their necessity for ensuring numerical convergence to weak solutions. This paper proposes a Lagrangian-AMR scheme based on the classical cell-centered Lagrangian scheme that preserves both geometrical and physical conservation laws. At level boundaries with hanging nodes, a unique nodal solver enforces geometric collinearity through a flux correction term, which is distributed to neighboring cells via energy-weighted averaging. The algorithm is implemented using the p4est library with density gradient-based adaptivity, and numerical results show good agreement with analytical solutions.

ABSTRACT. We propose a high-order discontinuous Galerkin (DG) method for solving the eikonal equation with improved stability and efficiency. Conventional high-order solvers rely on nonlinear equations with auxiliary variables, increasing complexity and cost. We extend a linearized vanishing-viscosity algorithm to a high-order DG method, but direct extension leads to projection-induced error amplification and oscillations near singularities. To address this, we terminate vanishing-viscosity iterations at an optimal stage and switch to an inviscid formulation to prevent error amplification in smooth regions. A robust criterion is also proposed to detect singular regions generated by characteristic front collisions. The Singular Gradient Regularization Method (SGRM) then applies viscosity only in singular regions while keeping smooth regions inviscid. A region-decoupled flux prevents viscosity contamination of smooth regions. The method uses a scalar-only formulation, avoiding auxiliary gradient equations and reducing computational cost. Extensive numerical experiments on complex 2D and 3D geometries demonstrate high-order accuracy and stable behavior near singular gradients. Integration into high-order RANS and hybrid RANS–ILES solvers confirms its effectiveness for practical wall distance computations.

ABSTRACT. The SLAU2 is one of Euler fluxes used to compute the numerical flux at cell interfaces in the finite volume method (FVM). This study compared performance of recently proposed methods and SLAU2 in our hypersonic flow test. It has been demonstrated that the recently proposed Euler fluxes, i.e., AUSM+M, AUSM2025p, and AUSM2025u, are as robust as SLAU2 against the notorious carbuncle phenomena or other shock-related anomalous solutions, despite numerically represented thinner shockwaves when the initial location exactly matches the grid lines. At this moment, it is impossible to determine whether this fact (the shock thickness is sensitive to its initial location) is desirable in terms of numerical robustness and theoretical explanations. Moving shock tests, in which the shock experiences any possible locations relative to the gridlines, will hopefully provide some clues. Apart from this, a shock/shock or shock/boundary-layer interaction examples are also of importance from an engineering application point of view. The final paper will include these extended test cases, as well as a matrix stability analysis. It will also cover other flux functions, such as AUSMDV2+ and LDFSS2025u/p/M. These will strongly contribute to further understanding and development of numerical flux functions.

ABSTRACT. This work presents a Streamline-Upwind Petrov–Galerkin (SUPG) stabilized finite element framework for predicting non-ionized hypersonic flows in thermal non-equilibrium. The formulation is enhanced with a residual-based discontinuity-capturing (DC) operator to robustly resolve steep gradients and shock-layer features. The framework solves the compressible reacting Navier–Stokes equations coupled with an additional conservation equation for vibrational–electronic energy, using pressure–primitive variables. The governing system includes one continuity equation for each chemical species, enabling finite-rate chemistry and its coupling to the two-temperature model. The approach is assessed on a set of benchmark cases that validate the chemical kinetics and the thermochemical coupling in the presence of critical flow phenomena. Accuracy is evaluated through comparisons with numerical and experimental data available in the literature, including hollow cylinder extended flare and double cone. The results demonstrate the suitability of the proposed formulation for simulating non-ionized reacting hypersonic flows under thermal non-equilibrium. We will also demonstrate the extension of the formulation to the ionized flows and three-temperature non-equilibrium.

ABSTRACT. Entropy stable discontinuous Galerkin (DG) methods display improved robustness for problems with shocks, turbulence, and under-resolved features by enforcing a mathematical entropy inequality. Such methods have traditionally relied on entropy conservative (EC) fluxes that can be computationally expensive to evaluate [1]. An alternative approach for enforcing an entropy inequality is through a minimally dissipative artificial viscosity called entropy correction artificial viscosity (ECAV). We review how to construct an artificial viscosity formulation using the Bassi-Rebay 1 (BR1) and Local Discontinuous Galerkin (LDG) discretization [2], as well as how to impose an entropy stable no-slip/slip adiabatic boundary condition. Then, we present numerical experiments for the Compressible Navier Stokes equations which suggest that ECAV method is capable of tackling difficult problems with complex shock and boundary layer interaction without requiring additional positivity preserving or shock-capturing methods. Finally, we discuss a flux-corrected transport (FCT) version of ECAV, as well as extensions of the proposed approach to non-ideal equation of states.

ABSTRACT. This work concerns the design and analysis of a limiting technique that allows the preservation of invariant domains for high-order numerical approximations of nonlinear hyperbolic systems of conservation laws. The method can be applied to any conservative discretization method in space and we here focus on time implicit integration for steady and unsteady problems. The method limits the high-order solution around a low-order accurate solution that is known to preserve all the invariant domains and can be seen as a generalization of the flux-corrected transport limiter [Boris and Book, J. Comput. Phys, 11, 1973; Kuzmin and Turek, J. Comput. Phys, 175, 2002] to systems of conservation laws.

We will give details on the derivation and implementation of this limiting technique and analyze the properties of the limited solution. Numerical illustrations will concern finite volume or discontinuous Galerkin space discretizations for nonlinear scalar equations and the compressible Euler equations.

ABSTRACT. Today, more than ever, physical and numerical modeling of atmospheric entry flows is one of the primary research focuses of the CFD community. During the entry phase, a body crosses the atmosphere and encounters a rarefied gas that becomes progressively denser with decreasing altitude. In the present work, the continuum regime is considered such that the gas is dense enough to be treated as a continuum medium and the Navier–Stokes equations apply. Strong bow shocks form in front of the body, and the temperature across them reaches thousands of kelvin, triggering non-equilibrium effects such as dissociation and ionization of the gas molecules. To account for these effects, the most common approach is to use the two-temperature model by Park. Although many research and commercial CFD codes use this model or some of its variants, it has been shown that it cannot reproduce the full spectrum of non-equilibrium gas dynamics. Hence, a reduced state- to-state model is considered for the simulation of atmospheric entry on Uranus and Neptune, and an operator-splitting technique is used to separate the physical processes for improved computational efficiency. The solutions are obtained with the CPU-based code SU2-NEMO and compared to those obtained with the GPU-based code DAMISO.

ABSTRACT. Reduced order modeling (ROM) is of paramount importance for efficiently monitoring and controlling complex systems in multi-scenario contexts, such as fluid flows. However, conventional ROM strategies are typically limited to known and constant parameters, inefficient for nonlinear and chaotic dynamics, and blind to the actual system behavior. In this work, we propose a sensor-driven Reduced Order Modeling strategy based on SHallow REcurrent Decoder networks (SHRED-ROM). Specifically, motivated by the well-known method of separation of variables, we consider the composition of a sequence model, which encodes the temporal dynamics of limited sensor measurements in multiple scenarios, and a decoder model, that reconstructs the corresponding high-dimensional spatio-temporal field. To enhance computational efficiency and memory usage, the snapshots dimensionality is reduced by data- or physics-driven basis expansions, allowing for compressive training of the networks with minimal hyperparameter tuning. By employing SHRED-ROM to control problems, it is possible to efficiently retrieve distributed control actions starting from sparse state readings in multiple scenarios. Through applications in fluid dynamics, we show that SHRED-ROM (i) accurately reconstructs and controls challenging dynamics in multiple scenarios with minimal sensor requirement, independently on sensor placement, (ii) can cope with both physical, geometrical and time-dependent parametric dependencies, while being agnostic to their actual values, (iii) can accurately estimate unknown parameters, and (iv) can deal with different data sources, such as high-fidelity simulations, coupled fields and videos.

ABSTRACT. Autoencoders (AEs) have emerged as powerful tools for non-linear dimensionality reduction, often surpassing traditional linear methods such as Proper Orthogonal Decomposition (POD) in scenarios characterized by slowly decaying Kolmogorov n-widths. In the realm of Reduced-Order Modelling (ROM), these models are increasingly utilized to learn low-dimensional representations of solution manifolds associated with parametric Partial Differential Equations (PDEs). However, the high expressivity of AEs presents a challenge: although trained networks typically minimize reconstruction error, they often struggle to capture the essential properties necessary for building accurate ROM. To overcome this issue, it is crucial to introduce bias during training, ensuring the model learns a pair of encoder and decoder that maintains representation consistency—specifically, constraining the composition of the encoder and decoder to be equal to the identity map E ◦ D= id within the local coordinates of the chart formed by the encoder and the solution manifold. Recent work has tackled this challenge in fully connected AEs by proposing representation-consistent architectures. This study builds upon that concept by extending representation consistency into the convolutional domain. We introduce a novel class of symmetric Convolutional Autoencoders (CAEs) designed to embody the primary properties of manifold parametrization mappings. When integrated into a ROM framework, this architecture demonstrates significantly improved predictive capabilities. Numerical experiments indicate that, in comparison to standard CAEs, our proposed approach achieves reduced reconstruction errors and generates more accurate latent-space trajectories.

ABSTRACT. 1. Introduction Urban fluid flow and air-quality simulations generate high-dimensional spatio-temporal datasets whose size makes storage, analysis, and real-time forecasting computationally demanding [1]. Reduced-order models (ROMs) such as high-order singular value decomposition (HOSVD) have emerged as a viable solution. ROMs provide better computational expense and accuracy [2], [3], as they reduce the dimensionality of flow fields while preserving important spatio-temporal data. HOSVD provides a natural tensor-based framework for extracting coherent multilinear structures while preserving the physical organization of multivariable fields. However, its adoption in large-scale applications is limited by the cost of repeated singular value decompositions of full tensor unfoldings. 2. Methodology This work introduces a low-cost HOSVD framework for efficient reduced-order modeling and prediction of urban flow cases. The method computes modal subspaces from reduced unfoldings obtained through sparse spatial sampling or limited measurements, thereby avoiding the assembly of the full spatio-temporal tensor and the repeated singular value decompositions of large-scale unfoldings. It then reconstructs the full-order factors via a lifting strategy. In this way, the dominant tensor structure is recovered without assembling the complete dataset, leading to substantial reductions in memory footprint and computational time while maintaining reconstruction accuracy. Beyond reconstruction, the resulting compact representations are directly compatible with data-driven forecasting tools. Low-cost HOSVD will be combined with deep learning architectures for temporal forecasting. This hybrid forecasting pipeline enables accurate reconstruction and prediction of multi-component velocity and scalar concentration fields at a fraction of the original computational cost. 3. Conclusions The proposed methodology provides a practical bridge between physics-based tensor decompositions and modern machine learning, enabling scalable reduced-order prediction tools for complex urban environments. Tests will be conducted on representative urban and environmental flow databases to demonstrate that high-fidelity results can be achieved using only a small subset of degrees of freedom, making the approach suitable for real-time and large-scale deployment.

ABSTRACT. Urban areas currently support more than 55% of the global population, a figure projected to reach 68% by 2050. High population density and restricted evacuation options elevate the risk of human exposure to airborne contaminants. Consequently, developing rapid response capabilities for the accidental or intentional release of chemical, biological, or radiological substances is a critical challenge. In these scenarios, the application of physics-based urban wind and dispersion models is essential for emergency management protocols and decision-making frameworks.

Ensuring the rapid turnaround required for emergency response necessitates models that are both highly automated and computationally efficient. Cartesian grid methods with embedded boundaries address the automation requirement by facilitating rapid, fully automated mesh generation for complex geometries. Beyond streamlined meshing, these methods offer flow-based adaptation, straightforward numerical formulations, and low numerical dissipation. Furthermore, the inherent regularity of Cartesian grids enables efficient vectorisation and is well-suited for high-performance parallel architectures such as GPUs. To complement these advantages and ensure rapid convergence, multi-level methods are employed to mitigate the deteriorating convergence rates typically exhibited by iterative solvers, such as the SIMPLE algorithm, as grid density increases. Specifically, the non-linear Full Approximation Scheme (FAS) accelerates convergence by efficiently resolving low-frequency error components on coarser grids. Previous investigations have demonstrated the effectiveness of the FAS multigrid approach for a variety of complex problems, including urban windflows.

To meet these demands, we developed a Cartesian grid solver that combines a cut-cell embedded boundary method with a multilevel algorithm to achieve both efficiency and automation. Developed using the AMReX framework, the method solves the steady-state, incompressible Reynolds-Averaged Navier-Stokes (RANS) equations using the SIMPLE algorithm. Turbulence is modelled with the Wilcox k-omega (2006) model, utilising a high-Reynolds number wall treatment for cut-cells. To automate the mesh generation, we use flow based mesh adaptation using sensors to detect areas that need refinement. In this short abstract, we outline briefly the implementation of the high-Re wall treatment and present some validation results.

ABSTRACT. The mechanical interactions between ultrasound-stimulated microbubbles (USMB) and cell membranes are a critical phenomenon in targeted drug delivery and immunotherapeutic strategies. The acoustically driven oscillation of the microbubbles has been shown to induce sonoporation, a transient increase in cell membrane porosity that enhances site-specific drug delivery. In contrast, inertial cavitation can irreversibly damage and kill cells, which should be avoided during USMB-based therapeutic and diagnostic procedures. Direct observation of ultrasound-stimulated microbubble (USMB) interactions with cells remains challenging due to the fast timescales (microseconds) and small dimensions (micrometers) involved. Nevertheless, optimizing ultrasound-based diagnostic and therapeutic techniques requires a detailed understanding of these interactions at the cellular interface. A comprehensive knowledge of how USMBs behave in proximity to cells is essential for enhancing both their clinical efficacy and safety. To this end, we have developed an interface-resolved axisymmetric computational model to investigate these interactions. In our model, the cell is represented as a droplet of specific viscosity and density suspended within a carrier fluid. A phase-field method is used to model the degradation of the intracellular structure due to the cavitation-induced strain. Our model captures the essential bubble-droplet interactions and interfacial mechanics necessary to characterize USMB-mediated cellular effects. The results show that the stably oscillating microbubbles translate toward the denser and more viscous droplets in the carrier fluid. As the bubble gets closer to the cell, the flow induced by its oscillation periodically deforms the droplet interface. Such deformation of the cell boundaries in cell-bubble interaction may instigate the formation of a temporary pore at the cell membranes, enhancing drug delivery through sonoporation. Furthermore, the acoustically driven bubbles are seen to form high-speed micro jets toward the droplets. The impingement of high-velocity microjets during the inertial cavitation of USMBs can cause irreversible damage, leading to cell lysis—an undesirable outcome in ultrasound-mediated drug delivery. The degradation model will track the damage caused to the intracellular structure due to microjetting. However, jetting collapse can be moderated by adjusting the ultrasound frequency. Our calculations demonstrate that microjet formation is suppressed at higher acoustic frequencies. These results indicate that the current model provides critical insights into cell-bubble interactions, serving as a framework for the development and optimization of USMB-based diagnostic and therapeutic techniques.

ABSTRACT. This work presents advancements in Computational Multi-Fluid Dynamics (CmFD) by proposing a thermally driven Marangoni flow model in a Sharp Interface Level Set Method on a co-located grid (SI-LSMcol). Further, the associated formulation, discretization, and implementation of the interfacial jump condition, essential on a co-located grid, are presented. Finally, validation of the in-house two-phase flow solver is validated for a range of benchmark CmFD problems, including thermocapillary single-fluid Stokes flow in an open cavity, linear shear two-fluid flow, and thermocapillary two-fluid flow in a closed cavity.

ABSTRACT. We present a systematically second-order accurate method for the local evaluation of hydrodynamic loads in immersed boundary simulations for fluid–structure interaction applications. The approach is developed within a finite-difference direct numerical simulation solver for the incompressible Navier–Stokes equations, coupled with a sharp-interface immersed boundary formulation.

The work identifies the discrete mechanism responsible for the well-known degradation of viscous stress accuracy near immersed interfaces. The loss of convergence is traced to a breakdown of local discrete error regularity at boundary-affected grid points, where stencil values no longer share a consistent discretization error structure and second-order cancellations fail.

To prevent this degradation, a probe-based reconstruction strategy is introduced in which pressure and velocity gradients are evaluated exclusively at regular fluid locations and subsequently extrapolated to the immersed surface. Probe placement explicitly avoids boundary-affected points, ensuring consistent second-order accuracy without empirical tuning of probe distances.

The method is verified through manufactured-solution tests and validated against reference numerical and experimental benchmarks, including the flow past a sphere in free stream and sphere sedimentation under gravity. Second-order convergence of reconstructed surface tractions is demonstrated.

ABSTRACT. Rocket-Based Combined Cycle (RBCC) engine integrates a high thrust-to-weight ratio rocket and a high specific impulse dual-mode ramjet within a single flow path, enabling efficient and economical flight across a wide speed range. This makes RBCC a key propulsion solution for future reusable space transportation vehicles. The combustion chamber of RBCC engine is persistently exposed to extreme conditions involving high temperature, high heat flux, and high-speed flow. An efficient Thermal Protection System (TPS) is essential to maintain structural temperatures within material limits, and thermal protection technology has become a critical bottleneck in engine development. Among various active cooling strategies, regenerative cooling has attracted significant research and engineering interest due to its high efficiency, stable wall temperature control, improved energy utilization, and potential to enhance fuel combustion performance, thereby supporting long-duration and reusable engine operation. Conventional regenerative cooling systems typically employ straight rectangular channels, favored for their structural simplicity and low cost. However, under multi-mode operating conditions, strong coupling between supersonic combustion and regenerative cooling results in a non-uniform distribution of thermal environment and heat load. Moreover, hydrocarbon fuels often operate under supercritical conditions, where drastic variations in thermophysical properties can induce flow maldistribution and localized overheating, potentially leading to engine failure. Therefore, optimizing cooling channel design is critical to improve flow distribution and overall cooling performance. Topology optimization has become a key methodology for designing channel configurations that achieve optimal performance under specified constraints, with substantial progress in fluid flow and conjugate heat transfer applications[1][2]. Yu et al. [3]employed the moving morphable components (MMC) method to design heat dissipation channels, achieving improved stability compared to density-based methods and reducing the formation of unmanufacturable micro-channels. In open-source platform development, U. Nilsson et al.[4] implemented a topology optimization framework in OpenFOAM using the continuous adjoint method for sensitivity analysis, systematically investigating adjoint equation derivation and discretization of adjoint boundary conditions. Matsumori et al. [5]optimized thermo-fluid channel configurations with the objective of maximizing heat transfer, proposing a design methodology for heat exchanger channels under constant power input and evaluating Reynolds number effects on topological complexity. Research from the Technical University of Denmark (DTU)[6] explored simplified models for fluid topology optimization in natural convection problems, showing that although boundary layer effects are not captured, such models provide effective initial guesses for high-fidelity optimization, reducing computational cost. Existing studies and our previous study[7] indicate that topology-optimized regenerative cooling channels can overcome key limitations of conventional straight channels, including flow maldistribution and localized heat transfer deficiencies. However, state-of -art researches focus on thermal management of cooling structure with small temperature variations. Large thermophysical property changes due to temperature increase from 300K to 900K in practical engine cooling applications make topology optimization of regenerative cooling channels much more difficult. In this paper, topology optimization model considering large property variation for supercritical hydrocarbon fuel via CoolProp is applied to the regenerative cooling channels of a typical ramjet combustor wall. Three-dimensional numerical simulations are conducted to compare flow and heat transfer characteristics between the optimized and conventional straight channels, validating that topology optimization enables the design of lighter regenerative cooling structures with enhanced heat dissipation efficiency. Figure 1 shows a schematic diagram of a typical structural component of a ramjet combustor. Its wall structure is extracted and simplified into the topology optimization domain for the present study. Owing to the symmetry of the physical problem, half of the full topology optimization domain is adopted as the actual optimization domain. The design domain consists of a 100 mm × 13.75 mm rectangular region and six inlets and outlets with a dimension of 4 mm × 2 mm, where the spacing between adjacent inlets and outlets is 2 mm. Each inlet and outlet is located 1 mm away from the symmetry plane and 2.75 mm away from the top surface of the design domain. A volumetric heat source with a magnitude of 1×10⁸ W/m³ is applied within the design domain. The inlet velocity is set to 0.05 m/s, and the inlet temperature is set to 300K. For topology optimization problems involving thermo-fluid coupling, the density method is commonly adopted for optimization. A uniform volumetric heat source is applied within the design domain, whose magnitude is determined by the average heat flux of the heated wall. With prescribed energy dissipation constraint and volume constraint, the optimal spatial distribution of fluid-solid material elements is solved via numerical calculation, so as to minimize the objective function of the average temperature in the domain. In the present work, the volume constraint is prescribed at 0.6. Based on the density method, the topology optimization formulation for the thermo-fluid problem with the objective of minimizing the average temperature can be expressed as follows: In the present work, the high-temperature alloy GH3128 and n-decane are selected as the solid and fluid materials, respectively. Given the drastic variations of n-decane’s thermophysical properties with pressure and temperature in the channel during regenerative cooling, accurate calculation of its thermophysical properties is required, for which polynomial fitting and tabulation with interpolation are the two common approaches. Herein, the tabulation with interpolation method is employed. The n-decane thermophysical property table is established via CoolProp at 3 K temperature intervals and 1 kPa pressure intervals without considering fuel pyrolysis, and the fluid thermophysical parameters in the numerical simulation are acquired by linear interpolation. The final topology optimization result from the calculation is shown in Fig. 2. The ramjet combustor wall is typically divided into an inner wall, a channel layer, and an outer wall. The channel layer is formed by extruding the topology optimization result, and both the inner and outer walls are fabricated from the high-temperature alloy GH3128, to establish the three-dimensional combustor wall model. Numerical simulations of flow and heat transfer are conducted on the three-dimensional topology-optimized regenerative cooling structure. The standard k-ε model is selected for the Reynolds-Averaged Navier-Stokes (RANS) calculations, which has been widely validated to accurately predict the flow and heat transfer of supercritical-pressure hydrocarbon fuel in horizontal channels. The Coupled algorithm is used for pressure-velocity coupling, with the convection and diffusion terms discretized by the second-order upwind and second-order central differencing schemes, respectively. Given the complex species and reaction processes involved in the endothermic pyrolysis of kerosene (a multi-component hydrocarbon fuel mixture), a surrogate component model is adopted for numerical calculation. A simplified three-step global reaction mechanism with 17 species is employed to simulate the endothermic pyrolysis of kerosene under supercritical pressures, with C₁₁.₈₅H₂₃.₈₂ as the kerosene surrogate. The simulation results of three-dimensional flow and heat transfer considering pyrolysis are shown in Fig. 3. The design of engine regenerative cooling channel via the topology optimization method holds great promise for improving flow distribution uniformity and heat transfer efficiency, thus enabling more efficient wall thermal protection. In the present work, three-dimensional numerical simulations are conducted to evaluate the cooling performance of the topology-optimized cooling panel under the operating heat flux conditions of an RBCC engine, with the main conclusions drawn as follows: (1) In terms of mass, the solid-domain fraction in the channel layer is 43% for the topology-optimized configuration and 58% for the straight-channel configuration, corresponding to a 26% mass reduction in the channel layer. (2) The topology-optimized channel has an average heated-surface temperature of 1000 K and a maximum temperature of 1090 K, both lower than those of the straight channels under different fluid volume fractions. This indicates that the topology-optimized configuration can effectively improve the wall-cooling capability of regenerative cooling channels. (3) Under the same volume fraction, the topology-optimized configuration still has a cooling advantage over the straight-channel configuration. The average and maximum temperatures are reduced by 17.5% and 16.5% under the low-temperature condition, by 5.6% and 6.4% under the medium-temperature condition, and by 3.3% and 4.4% under the high-temperature condition. (4) The cooling effect of the topology-optimized channel mainly originates from coolant redistribution by the branching network and enhanced transverse heat transfer. This advantage is most obvious under the low-temperature condition, and the relative temperature reduction decreases as the temperature level increases.The present study demonstrates that topology optimization serves as a feasible and effective computational design tool for regenerative cooling systems in ramjet engines. Future efforts will be directed toward the additive manufacturing of the topology-optimized cooling structure, followed by ground-based direct-connect thermal qualification tests of the ramjet combustor,with the aim of further facilitating the engineering implementation of topology optimization technology in the aerospace propulsion field.

ABSTRACT. Adjoint-based shape optimization in turbulent flows is developed by explicitly treating turbulent fluxes as forcing terms in the governing equations. Conventional RANS-based adjoint methods rely on eddy-viscosity models, which assume isotropic turbulence and may limit accuracy, particularly near walls. To overcome this limitation, turbulent statistics are first obtained from direct numerical simulation, and the Reynolds stresses and turbulent heat fluxes are incorporated into the adjoint framework as virtual forcing terms under the frozen turbulence hypothesis. The method is applied to the multi-objective optimization of a pin-fin configuration in turbulent channel flow, aiming to reduce drag while enhancing heat transfer. Sensitivity information is derived and used to update the geometry through a level-set method. The optimized design achieves a 13% reduction in drag and a 4% increase in heat transfer, demonstrating the effectiveness of the proposed framework for accurate turbulent shape optimization.

ABSTRACT. 1. Introduction Modern computational fluid dynamics (CFD) workflows increasingly rely on large numbers of parameterized simulation runs to support design exploration, uncertainty quantification (UQ), Non-deterministic Evaluation (NDE) and verification activities. Typical examples include angle-of-attack sweeps, grid refinement studies, and repeated analyses under varying operating conditions. While high-performance computing (HPC) platforms enable these studies to be executed efficiently, ensuring that results remain reproducible, auditable, and reusable over time remains a persistent challenge. A well-defined workflow for one instance of the parameterized study, that may even be tracked through version control, addresses how information is derived through a process since it essentially describes the utilized algorithms. Provenance addresses what information is derived through a process since it describes relationships within a set of files for each instance of the parameterized study including reference to the utilized workflow. Provenance, rather than workflow alone, enables an engineer to trace to specific files to investigate interesting nuances to gain insight into the overall results as well as investigate questionable results that drive CFD model refinement (e.g. flow visualization from volume solution files corresponding to one point on a UQ curve).

The current study addresses the problem of CFD results traceability and reproducibility by treating the set of files utilized in a process and associated relationships between them as a first-class computational object. The management, or data control, of a process’ files occurs through a provenance-aware execution method that explicitly records the relationships among simulation parameters, solver runs, generated artifacts, and post-processing steps as workflows execute. Furthermore, file access is virtualized, simplifying review of the many file artifacts used and created during workflow execution. The approach is demonstrated using a parametric CFD workflow based on the XFOIL panel method, representative of the repeated analyses commonly performed in preliminary aerodynamic design and UQ studies. The value of the provenance-aware method can be clearly understood for efficiently working through the use cases of large parameterized studies outlined below.

1. Traceability: Trace specific process instance files to: a. Find source file(s) for a runtime or logic error b. Understand how a CFD model should be refined c. Investigate specific run instances corresponding to interesting nuances of the overall results 2. Reproducibility: Invert retrospective provenance (what actually occurred) into a precise and detailed reproduction workflow to regenerate results rather than hoping for equivalent execution of a pre-defined workflow, which is prospective provenance (what was planned to occur)

2. Methodology 2.1 Workflow-Centric CFD Execution Model In the proposed approach, a CFD study is represented as an executable workflow composed of well-defined stages and run units, which defines how information will be derived. Each run corresponds to a solver invocation with a specific set of input parameters, while stages group related runs and post-processing steps. Key workflow elements include airfoil geometries, run parameters (e.g., Mach number, angle of attack), solver input and output files, scripts, and derived results such as plots or aggregated quantities of interest. By expressing these elements explicitly, the workflow structure becomes independent of any particular file system layout or scripting convention. 2.2 Provenance Capture During Execution Provenance information is captured automatically as workflows execute, recording both the intended workflow structure and the actual execution history, which defines what information is derived. This includes the generation and consumption of data artifacts, the parameters used for each solver run, the software components involved, and the ordering of workflow activities. Capturing provenance at execution time ensures that the full chain of derivation for any result can be reconstructed without relying on external documentation or manual bookkeeping. The capture process is designed to be open and non-intrusive, allowing engineers to run analyses using familiar tools and HPC environments while provenance is recorded transparently. 2.3 XFOIL-Based Parametric and UQ Workflow The methodology is demonstrated using a panel-method aerodynamic analysis implemented with MIT’s XFOIL code. In the example workflow, the airfoil geometry and Mach number are held fixed while the angle of attack is varied over a prescribed range, producing a sequence of solver runs whose results are subsequently consolidated and post-processed to generate aerodynamic polars. The second example is a UQ workflow, in which XFOIL runs are used to create a surrogate model for statistical sampling. These activities are used to perform sensitivity analysis (SOBEL indices) and compute confidence metrics. These examples are representative of common parametric sweeps and UQ studies in CFD practice, where tens or hundreds of solver runs must be managed consistently. 2.4 Results and Reproducibility Analysis Captured provenance allows the content for any generated artifact to be traced in terms of how and what information was derived or reproduced based on the exact parameters, solver version, and intermediate results involved in its creation. The provenance representation supports direct re-execution of the entire analysis suite or selected subsets with modified parameters, enabling controlled comparisons between runs. Differences between workflow execution instances—such as parameter changes or altered solver settings—can be identified explicitly, providing a concrete basis for in-depth investigations, reproducibility assessment and auditability in CFD studies. Organizing artifacts via their workflow provenance provides an additional benefit as users can access file artifacts directly by clicking on thumbnails in the provenance timeline. This eliminates the distraction of tracking down file names and folder paths when seeking details within complex workflows. 3. Conclusions This work demonstrates that explicit provenance capture integrated with workflow execution provides a practical foundation for efficient in-depth investigations to better understand a broad set of results and reproducible CFD analysis, particularly for parametric and uncertainty-driven studies. No additional burden is on the analyst for reliable rerunning, comparison, and auditing of simulation results when using this approach since the total provenance relationship structure spans the pre-defined workflow, all artifacts recorded during execution, users, computers and software, and execution activities that occur over time. The XFOIL-based examples will illustrate how common CFD workflows can be made reproducible and transparent, with direct relevance to verification, benchmarking, and UQ applications. Future work will extend this approach to more complex CFD solvers and multi-physics workflows.

ABSTRACT. This paper presents a gradient‑based aerodynamic shape optimization framework for minimizing ground‑level sonic boom noise using a partitioned nearfield-farfield analysis approach. High-fidelity Reynolds-averaged Navier–Stokes simulations are performed in the nearfield using a finite-difference discretization on curvilinear overset grids, while farfield acoustic propagation is modeled with ray acoustics and the augmented Burgers' equation solver. A coupled adjoint methodology provides efficient sensitivity calculations across the nearfield and farfield domains, enabling the use of ground‑level noise metrics as optimization objectives. The use of curvilinear overset grids allows the off‑body Mach‑cone‑aligned grid in the nearfield domain to remain fixed during optimization, preserving solution accuracy while deforming only the near-body grid. Three test cases demonstrate the framework’s capabilities: minimization of the ground signature of a two-dimensional diamond airfoil, minimization of the A-weighted sound exposure level for the JAXA wing‑body configuration, and the inverse design of the NASA Concept 25D low‑boom configuration. Results show good reduction of noise metrics, demonstrating the effectiveness of the proposed optimization methodology for low‑boom supersonic aircraft design.