ABSTRACT. Solving a PDE in a complicated domain with a Cartesian mesh leads one to consider "cut-cells"–Cartesian cells that intersect the boundary of the domain. This type of mesh can handle complicated geometry in a robust and automatic way. The difficulties of mesh generation are replaced with those of accuracy and stability at the cut cells. Many interesting ideas have been proposed.

This talk will review a variety of approaches to these problems. I will focus on our newer approach called State Redistribution (SRD), which stabilizes finite volume and Discontinuous Galerkin schemes in a practical post-processing step at every time step. I will introduce a new version we are developing called Delta Redistribution (DRD), to resolve one of the problems with SRD. Computations in two and three space dimensions will be shown. We end with what I see as the current bottlenecks, and a discussion of open problems.

ABSTRACT. Global well-posedness for three-dimensional fluid flow equations remains a profound open problem. Recent efforts have shifted toward statistical solutions as a robust and physically meaningful framework for describing turbulence; yet, efficient computational tools to explore these solutions in three dimensions are scarce. We develop novel probabilistic lattice Boltzmann methods to compute and analyze statistical solutions for three-dimensional incompressible flows. We combine various types of lattice Boltzmann methods with probabilistic approaches, such as Monte Carlo and stochastic Galerkin methods. Through a platform-agnostic implementation on heterogeneous high-performance computing systems, we efficiently target application-relevant flow configurations. We present numerical results demonstrating the computational exploration of statistical solutions and their associated observables, including Wasserstein distances. We validate the convergence of these statistical measures in periodic regimes and demonstrate the computational feasibility of extracting statistical solutions for complex, wall-bounded configurations where deterministic simulations are nonunique. All developed methodologies are implemented in the open-source framework OpenLB to ensure public accessibility and sustainable reusability, enabling applications such as training generative diffusion models.

ABSTRACT. 1. Introduction

Despite the increasing relevance of blunt-body configurations in transonic transport systems, established RANS methods often fail to reliably capture the unsteady shock dynamics that characterise this flow regime. While RANS and URANS can predict integrated quantities, such as total aerodynamic drag, with reasonable accuracy, they generally do not recover the spectral behaviour of shock-induced separation or the resulting unsteady loads that drive structural excitation and passenger discomfort. Previous work by the authors demonstrated that Hybrid Temporal Large Eddy Simulation (HTLES), developed by Duffal, de Laage de Meux, and Manceau [1, 2], can recover these unsteady shock dynamics that URANS cannot capture while remaining computationally cost-effective [3]. However, those results lacked experimental validation. This paper addresses that gap through wind-tunnel measurements, providing a first step toward experimental validation of HTLES for shock-wave/boundary-layer interaction (SWBLI) flows. First results indicate good agreement in both qualitative and quantitative measures, including the spectral behaviour of shock buffet, Fig. 1.

The authors' previous work focuses on evacuated tube train (ETT) concepts (often referred to as Hyperloop) as a case study [3]. There, we compared URANS against HTLES in conditions just above the isentropic limit (choking), assessing both integrated metrics and spectral behaviour. To ground those numerical findings and to begin validating HTLES for SWBLI flows, the present study adds targeted wind-tunnel measurements. Although a conventional transonic indraft wind tunnel cannot reproduce the exact boundary and operating conditions of an ETT system, it can reproduce key SWBLI phenomenology, enabling mechanism-based validation. This includes the occurrence of choking, the overall shock structure (e.g. a lambda-shock), and the associated shock motion, i.e., global unsteadiness modes similar to transonic shock buffet. The validation is performed through direct comparison between experimental high-speed schlieren-like imaging and numerical schlieren visualisations derived from the HTLES and URANS baseline simulations. For the experimental imaging system, we employ an infinite-fringe differential interferometry setup (Wollaston-prism shearing interferometry) with a monochrome high-speed camera. In this configuration, the resulting fringe pattern encodes a finite-difference approximation of the density gradient in the shear direction, which yields a schlieren-like visualisation for small shear angles. More details are provided in Sec. 2.1 (Experimental Methodology). These high-speed measurements enable both qualitative matching of the overall shock structure and quantitative analysis of shock-buffet frequencies. Visual schlieren observations are complemented by wall static-pressure measurements just upstream of the shock to compare pre-shock Mach numbers. The geometric setup and computational domain used in the HTLES and URANS simulations are based on a digital model[4] of the physical transonic and supersonic indraft wind tunnel at FH Aachen. Core sections, including the nozzle inlet, test section, and diffuser, are reproduced with minimal geometric simplification, while the extended ductwork between the diffuser and vacuum tank is represented only via its cross-sectional area rather than its exact geometry. The digital model allows us to control operating conditions and match boundary conditions between experiment and CFD precisely. More details on the numerical setup are given in Sec. 2.2 (Numerical Methodology).

2. Methodology

The validation methodology rests on two pillars: (a) high-speed differential interferometry (Fig. 2) in a transonic indraft wind tunnel and (b) scale-resolving CFD of the same facility. Both approaches yield schlieren-like density-gradient fields that can be compared directly, enabling mechanism-based validation based on non-intrusive measurements. The experimental and numerical setups are summarised below.

2.1 Experimental Methodology

Experiments are conducted at the FH Aachen transonic and supersonic wind-tunnel facility. The indraft tunnel is driven by the pressure differential between a 24 m^3 vacuum tank and a matching-volume, atmospheric-pressure, dehumidified air reservoir. The working-section assembly (nozzle, windowed test section, and diffuser) is 1.6 m long, of which 600 mm corresponds to the nozzle section. The nozzle consists of a fixed upper constant contour and an exchangeable lower insert. In conventional operation, this insert is used to set the test-section Mach number by adjusting the throat-to-test-section cross-sectional area ratio. In the present study, a custom-built insert is employed that does not merely precondition the test-section Mach number; instead, it acts as the test model itself, representing a planar aft section of an ETT pod as seen in Fig. 1. The Mach number upstream of the shock is a key parameter for assessing the operating condition and is determined from the wall static-pressure measurement and the reservoir pressure. The windowed test section, which contains the aft ETT pod geometry, provides optical access for flow visualisation. Here, we employ an infinite-fringe differential interferometry setup with a Photron Mini UX100-M high-speed camera operating at 10 kHz and an exposure time of 50 microseconds. A schematic of the optical arrangement is shown in Fig. 2. As the wind tunnel operates in an intermittent (indraft) mode, the vacuum-chamber pressure downstream of the test section varies during a run. The downstream diffuser is therefore used to control the back pressure acting on the aft end of the pod insert, ensuring controlled and repeatable operating conditions. In the present configuration, it is used primarily to obtain operating conditions just above the isentropic limit [5], i.e., close to choking.

2.2 Numerical Methodology

The numerical setup, including the computational case definition (geometry, mesh, and boundary conditions) and the solver configuration, closely follows the authors' previous paper on HTLES in SWBLI flows [3]. Hybrid Temporal Large Eddy Simulation (HTLES) is a hybrid RANS-LES approach that introduces an explicit temporal filter on top of the (implicit) spatial filtering inherent to LES. By defining a local temporal filter width, HTLES enables a continuous, smooth transition between URANS-like and LES-like behaviour without explicit zoning, which improves robustness and reduces sensitivity to user-defined switching criteria. This is particularly important for applied research and development efforts, where even early design phases require accurate spectral data on potential shock-excitation mechanisms to support rapid prototyping.

All simulations are performed in STAR-CCM+ (v2306), solving the fully coupled compressible Navier-Stokes equations for an ideal gas using an implicit, density-based coupled solver. Roe flux-difference splitting is used to ensure robustness in the presence of strong gradients such as shock waves. The spatial discretisation is finite volume, using a hybrid second-order upwind/bounded-central scheme (hybrid-BCD), which provides stability in RANS-dominated regions while keeping numerical dissipation low in resolved regions. We employ an implicit second-order BDF2 scheme for time stepping. Turbulence is treated using the Scale-Resolving Hybrid (SRH) implementation of HTLES (based on the 2019 implementation by Duffal et al. [1]), coupled with Menter's k-omega SST as the underlying RANS model. Near-wall turbulence is modelled on a low-y+ prism-layer mesh (targeting y+ <= 1), whereas the large-scale unsteadiness in the shock/separation and wake regions is directly resolved. The local URANS-LES blending is governed by the HTLES temporal filter width Delta T_F = 2*pi / omega_c, with omega_c(x_i,t) = min(pi / Delta t, U_s(x_i,t) * pi / Delta x(x_i)), such that the method transitions continuously between URANS-like and LES-like behaviour as local spatial/temporal resolution permits. The coupled system is advanced using an algebraic multigrid (AMG) linear solver (V-cycle, Gauss-Seidel relaxation) with a fixed number of five inner iterations per physical time step, selected after convergence testing to ensure residual reductions of three orders of magnitude per step. The mesh uses a trimmed hexahedral cut-cell/octree topology with targeted refinement in the shock and separation regions, and the time step is chosen from a convective CFL ~ 1 criterion, resulting in time step widths on the order of 10 microseconds to 50 microseconds on a O(10^7)-cell mesh.

3. Results

For the preliminary results, we focus on two operating conditions set by the diffuser. (i) Just above the isentropic limit (choking), the shock does not fill the entire cross-section, is highly unsteady, and repeatedly moves upstream before collapsing and reappearing further downstream. (ii) At further reduced back pressure, where the shock becomes more stable but still oscillates about a mean position.

These preliminary results indicate that HTLES matches the operating conditions well in terms of overall shock structure, mean position, and the spectral content of shock buffet. Figure 3 shows a single operating point above the isentropic limit with an established, stable shock and induced separation. The shock, in both form and position, matches between the experimental quasi-schlieren image (horizontal Wollaston prism), Fig. 3a, and the numerical schlieren, Fig. 3b. The entropy line marking the separated region, starting at the shock foot, is also visible in both images. Note that the numerical schlieren image is taken from an HTLES simulation with a 10 microsecond time step, providing a factor of 10 higher temporal resolution, whereas the experimental image is effectively averaged due to the 10 kHz frame rate and 50 microsecond exposure time. Spectral data further strengthen the comparison, with good agreement in shock-motion spectral behaviour.

--------------- List of Figures ---------------

Figure 1: Instantaneous numerical schlieren image of the shock-induced boundary-layer separation at the rear of a blunt-body ETT pod with a Q-criterion isosurface and Mach number colour scale. Figure 2: Differential interferometry flow visualisation setup with a full-spectrum white LED light source. The Wollaston prism splits the beam into the ordinary (o) and extraordinary (e) components with a shear angle epsilon and projects an interference pattern onto the camera sensor. With both Wollaston prisms placed at their respective focal points, the system operates in infinite-fringe mode, resulting in a schlieren-like image. Figure 3: Comparison of (a) experimental quasi-schlieren and (b) numerical schlieren from HTLES showing the shock structure and separation region at the target operating condition (just above the isentropic limit).

--------------- References ---------------

[1] Duffal, V., de Laage Meux, B., & Manceau, R. (2019). Development and Validation of a Hybrid RANS-LES Approach Based on Temporal Filtering. In Volume 2: Computational Fluid Dynamics. American Society of Mechanical Engineers. https://doi.org/10.1115/AJKFluids2019-4937 [2] Duffal, V., de Laage Meux, B., & Manceau, R. (2022). Development and Validation of a New Formulation of Hybrid Temporal Large Eddy Simulation. Flow, Turbulence and Combustion, 108(1), 1-42. https://doi.org/10.1007/s10494-021-00264-z [3] Geiben, B., Havermann, M., Hale, E., & Bil, C. (2025). Application of hybrid-temporal LES to shock-induced boundary layer separation in hyperloop flows. CEAS Aeronautical Journal. https://doi.org/10.1007/s13272-025-00921-3 [4] Kritzinger, W., Karner, M., Traar, G., Henjes, J., & Sihn, W. (2018). Digital Twin in manufacturing: A categorical literature review and classification. IFAC-PapersOnLine, 51(11), 1016-1022. 16th IFAC Symposium on Information Control Problems in Manufacturing INCOM 2018. https://doi.org/10.1016/j.ifacol.2018.08.474 [5] Lang, A. J., Connolly, D. P., de Boer, G., Shahpar, S., Hinchliffe, B., & Gilkeson, C. A. (2024). A review of Hyperloop aerodynamics. Computers & Fluids, 273, 106202. https://doi.org/10.1016/j.compfluid.2024.106202

ABSTRACT. Turbulent mixing induced by converging spherical shock waves is a canonical yet challenging problem for compressible variable-density turbulence modeling, where rapid distortion, strong dilatation, and pronounced density gradients coexist. In such flows, conventional eddy-viscosity closures often show limited predictive capability for anisotropy and mass-flux-related transport, while full Reynolds-stress transport models may be computationally demanding and sensitive to modeling assumptions. Motivated by the need for a practical closure suitable for CFD applications, this work develops an algebraic-stress-based turbulent mixing model tailored to converging spherical shock-driven configurations.

Starting from a Reynolds stress model, we derive an algebraic stress relation. Instead of expanding the Reynolds stress anisotropy tensor into a large number of tensor basis functions, which is common in many algebraic stress models, we leverage the spherical symmetry of the problem to directly solve for the Reynolds stress anisotropy. The obtained Reynolds stress constitutive equation is then evaluated with a high-fidelity dataset from Implicit Large Eddy Simulations (ILES) of a spherically converging Richtmyer-Meshkov instability, demonstrating good agreement with reference data for Reynolds stress anisotropy.

To implement the proposed algebraic stress closure in a practical turbulence model, we embed it into a modified BHR2 framework and obtain a revised k-l formulation. An algebraic closure is further adopted for the mass-flux term, which eliminates the need to solve the mass-flux and density self-correlation transport equations in the original BHR2 model. In this way, the approach retains the richer physical mechanisms of BHR2 while reducing model complexity. Using the modified k-l model together with the derived algebraic stress constitutive relation, we simulate a spherical converging case, where the evolution of bubble/spike positions and the distribution of volume fraction agree well with reference data, providing initial evidence of reliability of the proposed method.

ABSTRACT. Introduction: Turbulence models are a crucial component of industrial Computational Fluid Dynamics (CFD) due to the widespread usage of Reynolds-Averaged Naiver-Stokes (RANS) simulations to predict various flow phenomena. When introducing a new solver for industrial applications, a robust implementation of these models is crucial. Furthermore, ensuring the accuracy of the software is also paramount as industries require certification of their software to ensure reliability and compliance with regulatory standards. Hence, a new solver must demonstrate both a robust implementation and high accuracy to meet the stringent demands of industrial CFD.

CODA (Leicht et al, 2016) is the CFD software being developed as part of a collaboration between the French Aerospace Lab (ONERA), the German Aerospace Centre (DLR), Airbus and their European research partners. CODA offers the widely used eddy-viscosity models such as the one equation Spalart-Allmaras model (Spalart,1992) and the two equation models like Menter's Shear Stress Transport (SST) (Menter, 1994) model. Apart from them, CODA also offers more complex models such as the SSG/LRR-Reynolds Stress Models (RSM) (Eisfeld, 2005) for the resolution of complex flow features such as three-dimensional separation, vortices, etc. CODA has a combined framework for both Finite Volume (FV) and Discrete Galerkin (DG) discretizations employing strong implicit numerical methods with advanced features like automatic differentiation (AD).

However, obtaining a robust implementation in strongly implicit solvers for complex turbulence models like SST and RSM is a significant challenge (Langer, 2020). One of the major obstacles would be to ensure the positivity preservation of the turbulence variables (Ilinca, 1998). While, several legacy solvers rely on clipping mechanisms (Andre, 1976) to obtain positivity preservation, the same approach may not be suitable for implicit solvers (Langer, 2023). Alternative techniques such as solving the turbulence equations in their natural logarithmic form was provided by (Iilinca, 1998). In the framework of two-equation and Reynolds stress models, many solvers restrict themselves to the usage of the logarithmic form of the length-scale variable (Stefanski, 2020) (Braun, 2025) to avoid the usage of wall functions as the turbulent kinetic energy and Reynolds stresses are zero at the walls, thereby leading to an infinite boundary condition. This approach avoids clipping of the length scale variable which is beneficial for implicit solvers and AD.

However, when transitioning from legacy solvers with the non-logarithmic (non-log) form of implementations to a new solver with logarithmic (log) form of implementations, several new questions become relevant. The effect of the usage of the logarithmic form of the length scale equation in industrial type meshes has not been studied in detail. The studies such as the influence of log form on mesh refinement, discretization settings, far-field decay, effect on blending functions, usage of other limiters such as realizability constraints have not been given attention. This work aims to summarize the recent extensive studies carried out to understand the log form implementation of the SST model and to explore its application to the SSG/LRR RSM model in CODA.

Methodology:

To comprehensively assess the effects of the log form of equations, we consider a range of flow regimes using four different 2D cases and a 3D case - a subsonic 2D NACA4412 airfoil, transonic AC 05-12p6 2D airfoil, high-lift 2D multi-element airfoil and 3D ONERA M6 wing. To ensure a thorough evaluation, we utilize grid families for each of these cases. For the subsonic and high-lift cases, we employ the grid family from the NASA Turbulence Modelling Resource (TMR) website. We use an in-house mesh family for the 2D transonic AC 05-12p6 which is being developed as a part of the DLR-ACTIVATE project. We analyse the mesh convergence for these cases by examining the the boundary integral quantities like coefficient of drag (C_D), coefficient of lift (C_L) and also the surface integral quantities like coefficient of pressure (C_p) and coefficient of skin friction (C_f).

During our studies, we observed a dependence of the log form on discretization settings specifically the use of gradient limiters on the turbulence equations. We investigate why such effects appear only for the log form of the model.

On the finest grid level, we still observed differences in the activation of the blending function and the impact of other limiters like realizability constraints. Although the log form does not theoretically require limiters, we found that employing a realizability limiter, i.e. Durbin's limiter helps to obtain robust convergence. Furthermore, the interaction between the log form and the realizability limiters may differ from that of the non-log form. The outcome of these investigations are on-going and the conclusions for practical application in an industrial environment will be presented in the conference paper.

Conclusions:

Our findings are consistent with the conclusions drawn from studies performed using the DLR TRACE solver (Müller, 2018) (Müller, 2022). We observe that the logarithmic form of the length-scale equation exhibits improved robustness compared to the non-logarithmic form. Furthermore, its interaction with the flow field can have significant effects, which depend on the specific flow regime. The prediction of aerodynamic quantities may also vary, with the magnitude of the differences depending on the particular case and level of mesh refinement. By understanding these effects, we have been able to refine our industrial best practices with confidence, enabling the provision of accurate and reliable predictions.

ABSTRACT. Designing accurate and robust convection schemes for convection‑dominated flows remains a significant challenge. The difficulty arises primarily because sharp discontinuities and large-gradient fields can develop even when the initial conditions are sufficiently smooth. Typical flow features that involve steep gradients include shock waves, material interfaces, thermal and species gradients, shear layers, and reaction fronts. To obtain high‑fidelity simulations of convection‑dominated problems, convection schemes should satisfy at least the following properties:

The high-resolution property. A numerical scheme must accurately capture both small-scale flow structures and sharp discontinuities, which are often smeared by excessive numerical dissipation. The Kelvin–Helmholtz instability is a representative example where small-scale vortical structures coexist with sharp gradients; excessive dissipation can suppress these vortices and significantly degrade solution quality.

The scalar-structure-preserving property. A scheme should produce oscillation-free and physically meaningful solutions. Non-physical oscillations can lead to significant errors or even computational breakdown. Furthermore, the structural integrity of passively transported scalars—such as mass fractions, species concentrations, or volume fractions—should be preserved. Numerical evidence shows that although traditional total-variation-diminishing (TVD) limiters suppress oscillations, they often distort scalar profiles (with the Superbee limiter being a notable example). Preserving discrete physical bounds, such as positivity of density or specific internal energy, is also referred to as maintaining an invariant domain (IDP)

However, constructing a scheme that simultaneously satisfies these properties is nontrivial due to Godunov’s theorem, which states that no linear convection scheme of second order or higher can be monotone. Beyond the two essential properties above, desirable additional features include the use of compact stencils and straightforward extensibility to unstructured grids. In this work, we present our recent efforts toward developing high-resolution, scalar-structure-preserving convection schemes employing compact stencils. The new family of proposed schemes in this work is termed ROUND (Reconstruction Operators on Unified-Normalise-variable Diagram).

This work develops a new family of compact, high-resolution, scalar-structure-preserving schemes using normalised variable diagrams and neural networks. The proposed schemes are further extended to unstructured grids and applied to a finite‑volume subcell limiting framework for the discontinuous Galerkin (DG) method.

ABSTRACT. After numerous efforts aimed at achieving entropy-conserving formulations for the convective terms [1,2], which have been proved to enhance robustness and fidelity of turbulent, compressible fluid-flow computations, the present contribution addresses the issue of entropy stability [3,4] for the diffusive terms. In a discrete framework, the entropy equation induced by the integration of the main system of equations does not necessarily preserve the correct global monotonicity, established by fundamental physical principles. Indeed, depending on the discretization adopted for the mass, momentum, and total energy equations — typically the set being solved — the entropy production over domains with vanishing boundary effects may either not evolve according to its physical trend or even become decreasing. Moreover, the combination of different discrete structures may also compromise the conservation of total energy. In this work, we focus on deriving appropriate split formulations that remain consistent with such physical principles. In particular, by including entropy-conservative formulations for the convective terms, relatively simple and intuitive split formulations are derived to guarantee the consistency of the viscous terms, thereby yielding a fully entropy-consistent numerical framework for spatial discretizations.

The compressible Navier-Stokes equations for viscous, heat-conducting Newtonian fluids are expressed as a system of conservation laws for mass, momentum and total energy for unit mass and volume. Its evolution is driven by the convective and diffusive flux-vectors and the state-variables are naturally conserved in domains with vanishing boundary effects. Discretization of such system may, however, affect how induced quantities - i.e. those ones whose equations are not directly discretized - behave in the discrete framework. Objective of the present work is the correct reproduction of the entropy equation. Specifically, consistency with the continuous physical principles requires that, for periodic domain and in presence of viscosity and heat conduction, volumetric thermodynamic entropy is non-decreasing in time.

To achieve such induced condition without sacrificing total-energy conservation, entropy-conservative convective discretizations - i.e. those giving entropy-conservation in the inviscid limit - are employed for the convective contributions. Subsequent application of summation-by-parts differencing highlights the entropy-stable conditions, to be discretely fulfilled in the integral sense by means of the proposed splittings, tailored on primitive variable gradients and afterwards associated to flux forms with specialized, simple averages.

Finally, representative numerical tests and comparisons will be presented with respect to different existing methods - such as the Laplacian form [5] and non-split schemes based on either primitive or entropy variable gradients - highlighting each approach’s respective performance with respect to resolution and scheme order.

This work presents a fully entropy-consistent numerical framework for the spatial discretization of the compressible Navier–Stokes equations, ensuring physically correct entropy production while being total-energy conserving. By employing entropy-conservative formulations for the convective terms, simple and intuitive split formulations for the diffusive contributions are derived that satisfy the discrete entropy-stability conditions. Numerical experiments will aim at showing that the proposed approach achieves an improved balance between robustness and accuracy when compared with existing methods, making it well suited for the simulation of viscous and turbulent compressible flows.

[1] C. De Michele, A. Edoh, and G. Coppola. Finite-difference compatible entropy-conserving schemes for the compressible Euler equations. Journal of Computational Physics, 540(114262), 2025. [2] A. Aiello, C. De Michele, and G. Coppola. Entropy conservative discretization of compressible Euler equations with an arbitrary equation of state. Journal of Computational Physics, 528(113836), 2025. [3] E. Tadmor. The Numerical Viscosity of Entropy Stable Schemes for Systems of Conservation Laws. Mathematics of Computation, 49(179):91–103, 1987. [4] A. Peyvan, K. Shukla, J. Chan, and G. Karniadakis. High-Order Methods for Hypersonic Flows with Strong Shocks and Real Chemistry. Journal of Computational Physics, 490(112310), 2023. [5] M. Bernardini, D. Modesti, F. Salvadore, and S. Pirozzoli. STREAmS: A high-fidelity accelerated solver for direct numerical simulation of compressible turbulent flows. Computer Physics Communications, 263(107906), 2021

ABSTRACT. This work presents the GPU-porting of a high‑order modal Discontinuous Galerkin (DG) solver equipped with a structure‑preserving entropy‑aware formulation. The baseline scheme employs a conservative discretization of the compressible Navier–Stokes equations, which however exhibits severe robustness limitations as compressible phenomena take place. To obviate this an efficient projection–correction strategy to enforce entropy conservation/stability is embedded in the solver. The GPU porting is implemented using CUDA Fortran, enabling the offloading of all computationally intensive routines to the device, keeping only preprocessing operations on the CPU. The performance and accuracy of the resulting implementation are assessed through the Re=1600 Taylor–Green vortex, simulated in the near‑incompressible regime using three mesh configurations at polynomial order n=3. Numerical experiments demonstrate that the GPU‑accelerated solver achieves substantial speedups—of order x7 over the CPU baseline—while maintaining high fidelity in the predicted kinetic‑energy decay and dissipation profiles. In addition to single‑device benchmarks, parallel scalability tests conducted on the LEONARDO supercomputer reveal near‑optimal strong scaling up to 64 GPUs, confirming the suitability of the proposed DG formulation for large‑scale high‑order CFD simulations on modern many‑core architectures.

ABSTRACT. 1. Introduction

Performing accurate and reliable simulations of compressible flows, especially in the turbulent regime, is an arduous challenge from the numerical point of view, even in the absence of shock waves. Minimization of dispersion errors typically requires high-order methods, whereas minimization of numerical diffusion, which is often necessary in turbulent non-shocked regions of the flow field, is achieved through the use of central schemes. However, high-order central discretizations are subject to several instability problems, typically originating from the discretization of the nonlinear convective terms, whose importance increases as the complexity of the simulation grows. In recent years, structure-preserving, or physically compatible numerical methods have gained great attention from the numerical community for their ability to keep under control some of these sources of instabilities, and are nowadays an active research field.

Kinetic Energy Preserving (KEP) [1] methods are the most well-known and widely used among the class of structure-preserving discretizations. Their ability to nullify the spurious production of kinetic energy coming from the discretized convective terms constitutes a fundamental element for the suppression of nonlinear instabilities. More recently, also “affordable” Entropy Conservative (EC) methods have been developed [2], which add the important structural property that the numerical discretization does not alter the correct induced balance of entropy due to convective terms. Finally, the enforcement of the Pressure Equilibrium Preserving (PEP) property [3] i.e. the ability of the discrete scheme to suppress spurious pressure oscillations and to reproduce the traveling density wave solutions of the compressible flow equations, is another important requirement, especially in the case of variable specific heats, arising in non-ideal and/or multicomponent gas flows.

The simultaneous enforcement of some (or all) of these structural properties, in addition to the discrete conservation of primary invariants, is an important goal for the design of efficient and reliable modern numerical discretizations of compressible Navier–Stokes equations. As regards PEP formulations, which are the topic of the present contribution, several KEP and PEP schemes have been proposed in recent years in the case of single-component calorically perfect gas flows [2, 4, 5], which have shown increased robustness with respect to standard discretizations. However, when turning to non-ideal or multicomponent flows, only partial or approximate solutions have been proposed [6, 7]. The construction of a KEP and PEP discretization able to automatically preserve also the primary invariants mass, momentum, and total energy is still an open problem for non-ideal gases, with a potential great impact on the reliability of the numerical solvers. In this contribution, we discuss the extension of the PEP methods to the case of arbitrary gas models, still retaining full preservation of primary invariants and kinetic energy by convection.

2. Methodology

The analysis is conducted by considering the induced discrete evolution equations for the velocity u and the pressure p for the case of compressible Euler equations with an arbitrary equation of state. Assuming semidiscretization, with negligible errors coming from the temporal integration, the spatial error associated with velocity and pressure evolution is evaluated by analytically manipulating only temporal derivatives [4]. In the case of a multicomponent gas with arbitrary equations of state, the general condition for discrete Pressure Equilibrium is given in [7] and has been derived for a general (but known) functional dependence for internal energy as a function of densities and pressure, resulting in a global constraint for the numerical fluxes for internal energy and partial densities in the case in which u and p are uniform in space. Exact solutions to this constraint (satisfying also primary invariants conservation and Kinetic Energy preservation) have not been obtained so far. Only partial solutions are available at the moment, either satisfying pressure equilibrium approximately, or exactly but without exact conservation of total energy [6].

In this contribution, a formulation which exactly satisfies pressure equilibrium in the case of a single compound with a thermally perfect or real gas model is derived, which also discretely preserves kinetic energy and enforces exact conservation of mass, momentum and total energy in non-viscous, shock-free regions of the flow field. In contrast to usual approaches, the numerical PEP fluxes are obtained by using nonlinear averages induced by the particular functional dependence of the internal energy with respect to density and pressure. The final method still uses a conservative discretization for mass, momentum, and total energy, naturally enforcing global and local conservation of primary invariants, and guarantees kinetic energy preservation by a coordinate discretization of mass and momentum equations. The formulation is tested on controlled nonviscous simulations showing exact theoretical pressure equilibrium.

In Fig. 1 we report a snapshot of the pressure field for a two-dimensional inviscid double-jet flow for a Peng--Robinson gas at transcritical conditions. The computations are performed by using a newly proposed formulation (Fig. 1(a)) and a standard one [3] (Fig.~1(b)), confirming the ability of the proposed scheme in suppressing the pressure oscillations contaminating the solution computed with standard schemes.

References

[1] G. Coppola, F. Capuano, S. Pirozzoli, and L. de Luca. Numerically stable formulations of convective terms for turbulent compressible flows. J. Comput. Phys., 382:86–104, 2019.

[2] H. Ranocha. Entropy conserving and kinetic energy preserving numerical methods for the Euler equations using summation-by-parts operators. In S. J. Sherwin, D. Moxey, J. Peiró, P. E. Vincent, and C. Schwab, editors, Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2018. Lecture Notes in Computational Science and Engineering, volume 134. Springer, Cham, 2020.

[3] N. Shima, Y. Kuya, Y. Tamaki, and S. Kawai. Preventing spurious pressure oscillations in split convective form discretization for compressible flows. J. Comput. Phys., 427:110060, 2021.

[4] C. De Michele and G. Coppola. Novel pressure-equilibrium and kinetic-energy preserving fluxes for compressible flows based on the harmonic mean. J. Comput. Phys., 518:113338, 2024.

[5] C. De Michele, A. K. Edoh, and G. Coppola. Finite-difference compatible entropy-conserving schemes for the compressible Euler equations. J. Comput. Phys., page 114262, July 2025.

[6] M. Bernades, L. Jofre, and F. Capuano. Kinetic-energy- and pressure-equilibrium-preserving schemes for real-gas turbulence in the transcritical regime. J. Comput. Phys., 493:112477, 2023.

[7] H. Terashima, N. Ly, and M. Ihme. Approximately pressure-equilibrium-preserving scheme for fully conservative simulations of compressible multi-species and real-fluid interfacial flows. J. Comput. Phys., 524:113701, 2025.

ABSTRACT. Discontinuities and highly anisotropic flow features represent a considerable numerical challenge for high order discretization techniques. These challenges appear commonly in hypersonic flows, where extremely thin boundary layers and shocks are difficult to resolve, yet extremely important to represent accurately to effectively predict surface heating and shear forces. We discuss the application of Moving Discontinuous Galerkin with Interface Condition Enforcement (MDG-ICE) methods to viscous conservation laws. These methods promise to accurately represent extremely anisotropic features and ultra-thin regions with large gradients. We present a novel formulation of MDG-ICE for viscous conservation laws and contrast it with existing formulations of MDG-ICE. We validate this new approach with three test problems, designed to model the numerical challenges encountered in hypersonic flows, where without additional treatment, standard DG methods succumb to spurious oscillations.

ABSTRACT. Introduction

Continuous progress in adjoint-based optimization and automatic differentiation has made it possible to compute gradients of PDE-based objectives efficiently with respect to input parameters such as surface meshes. Most adjoint-based optimization approaches focus on shape optimization via surface meshes or control points. While the usage of gradient-based optimization algorithms significantly reduces the curse of dimensionality, the optimizations still need many evaluations of the objective function. Therefore, most implementations are either restricted to two dimensional cases, a small number of optimization parameters or coarse discretizations.

Simultaneously, advances in CAD parameterization and differentiable geometry kernels enable the direct differentiation of geometric models. We plan to combine CAD parameterization and adjoint-based optimization into a scalable workflow. Instead of shape optimization of meshes or control points, the CAD parameters should be optimized directly.

This approach promises multiple advantages compared to current implementations. First, the number of optimization parameters drops significantly. Instead of optimizing the coordinates of hundreds of thousands of surface mesh nodes, only the design parameters chosen by the engineer are relevant. Physically relevant parameters like angles of attack, span width, or dihedral angles retain their physical meaning, while also reducing the dimensionality of the design space by orders of magnitude.

Another advantage offered by this approach is the manufacturability of the end product. Since the engineer chooses the design parameters and their effect through the CAD software, the generated designs remain mostly manufacturable. No kinks, buckles, or other non-manufacturable shapes are introduced through the movement of boundary nodes.

In this talk, we focus on novel methods required to realize this concept in practice. In particular, we address the differentiable treatment of mesh deformation and its integration into the adaptive refinement framework used for this project.

Mesh Deformation as a Gradient Propagation Problem

The adjoint solver in this proposed workflow needs a fully differentiable mesh, but the CAD geometry only provides gradients on the surface. Hence, the surface gradients, obtained from the differentiated CAD geometry, which in turn depends on the location on the geometry and the engineering design parameters, have to be propagated into the whole mesh - from the surface nodes to the interior nodes.

A very similar problem is solved by mesh deformation algorithms. In geometry-dependent simulations, mesh deformation is commonly employed to accommodate boundary motion induced by changes in the underlying geometric description. In many existing workflows, mesh deformation is treated as a purely geometric preprocessing step. Its primary objective is to preserve mesh quality by propagating movement of boundary nodes depending on their location to the inner nodes. However, in gradient-based optimization, mesh deformation can play an additional role: it can define how the geometric gradients are propagated from the boundary into the interior of the computational domain.

From this perspective, mesh deformation can be interpreted as a gradient propagation problem. Variations of CAD design parameters induce displacements of boundary representations, which must be consistently extended to all interior mesh entities in order to maintain a differentiable mapping between design space and discretized PDE operators.

Integration with Tree-Based Adaptive Mesh Refinement

The proposed algorithm is implemented within the broader Adaptive Mesh Refinement (AMR) framework provided by the t8code library, which uses the forest-of-trees approach. In this framework, a coarse, conforming input mesh is refined using tree-based algorithms for each input cell, effectively creating a structured, tree-based mesh within every cell. See Fig.1 and [Holke2018] for an illustration.

Within this semi-structured approach, the geometry of the refined elements is derived directly from the coarse input mesh, which means that the mesh deformation algorithm only needs to be applied to this coarse mesh. This represents a potential advantage compared to mesh deformation implementations on fully unstructured meshes, where complexity typically scales with the number of mesh nodes. Preliminary results support this advantage, though a detailed evaluation is still ongoing. The first examples of CAD-based mesh deformation using AMR are shown in Fig.2.

ABSTRACT. Dynamic mesh deformation is essential for unsteady computational fluid dynamics (CFD) with moving boundaries, such as fluid–structure interaction and shape optimization. However, existing methods face a persistent trade-off between computational efficiency and mesh quality. Physics-based analogy methods are robust and preserve element quality, but they are computationally expensive. Conversely, fast algebraic interpolation methods are efficient, yet they can induce severe distortion and may lead to element inversion under distributed motions or large deformation amplitudes. We propose a Transformer-based surrogate that predicts interior displacement fields directly from boundary actuation, trained on high-quality reference deformations. The model is trained with objectives that promote geometric consistency and enforce key constraints. A hard projection step then imposes the boundary conditions exactly. Moving boundaries follow the prescribed displacements, and fixed boundaries remain stationary. This design reduces the risk of element inversion. Numerical studies indicate that the proposed approach substantially reduces per-update cost while preserving deformation quality comparable to robust physics-analogy methods. The method remains stable under distributed motions and large deformation amplitudes, and it shows encouraging generalization across actuation patterns and mesh resolutions. These results support efficient and reliable dynamic meshing for practical moving-boundary simulation workflows.

ABSTRACT. Accurate simulation of fluid flows with moving boundaries is essential in many engineering applications, yet remains challenging when boundary displacements become large. Arbitrary Lagrangian–Eulerian (ALE) methods provide an accurate representation of the fluid domain and near-wall physics, but their robustness and computational efficiency are often limited by the mesh motion strategy.

This work presents a novel ALE mesh motion formulation based on a fully nonlinear fictitious structural problem recast as a hyperbolic equation. Unlike classical elliptic approaches, the proposed method allows for an explicit time integration of the mesh motion, avoiding any global linear or nonlinear system inversion. This significantly improves computational efficiency and parallel scalability while maintaining robustness under extreme mesh deformations.

The approach is implemented in the open-source CFD software TrioCFD, relying on a flexible description of the fictitious material behavior through the MGIS/MFront interface. This design minimizes additional code development and enables the use of advanced nonlinear constitutive laws for mesh motion.

Validation and demonstration cases include verification of the geometric conservation law, benchmark tests with large boundary displacements, and challenging configurations involving severe mesh distortion. The results demonstrate that the proposed strategy effectively preserves mesh quality and ensures accurate physical solutions, extending the practical applicability of ALE methods to complex moving-boundary problems.

ABSTRACT. We numerically investigate the noise reduction performance of a NACA-0012 airfoil at a low subsonic Mach number and a moderate Reynolds number. Noise suppression at the trailing edge is targeted by minimizing the wall pressure jump, which is stipulated by the unsteady Kutta condition. Guided by this principle, we propose an effective trailing edge configuration capable of significantly suppressing the radiated noise. Numerical results reveal that a pure streamwise-decreasing impedance trailing edge (ITE) fails to provide effective noise reduction across a broad range of characteristic frequencies; in some cases, it even leads to a significant increase in noise radiation levels. Therefore, we develop an extended trailing edge (ETE) configuration, which achieves a notable noise reduction of 3.62 dB. This performance is further improved to 4.22 dB with the incorporation of an impedance design, forming an extended impedance trailing edge (EITE). The reduction in radiated noise is accompanied by modifications in the local flow near the EITE, primarily characterized by a decrease of the streamwise and spanwise Reynolds stresses and a rapid decay of coherent structures. With the aid of high-fidelity near-wall data, the distribution patterns of acoustic sources in the wavenumber/frequency space are analysed in detail using a multi-faceted Kutta condition framework, thereby isolating the distinct role of different elements in achieving the desired noise reduction. For the ETE, it is found that noise suppression is attributed to two primary effects: (\textit{i}) enhanced destructive interference between sound waves scattered by wall sources near the trailing edges, and (\textit{ii}) a significant suppression of the vortex shedding process, particularly as evidenced by the distribution of near-field longitudinal process sources. Furthermore, extra noise attenuation is attained by the use of EITE, as the impedance surface leads to a lower convective velocity of pressure fluctuations and reduced energy within the supersonic phase speed region, consequently contributing to a lower scattering efficiency.

ABSTRACT. This paper introduces OpenCFD-FWH, an open-source implementation of the Ffowcs Williams-Hawkings (FW-H) acoustic analogy for permeable surfaces under incident mean flow. The method is developed for coupling with the OpenCFD-EC solver. Convective effects are accounted for by transforming the problem into a moving reference frame using the Garrick triangle approach. A rotational transformation further corrects angle-of-attack (AoA) effects. This formulation simplifies the FW-H surface integral and improves computational efficiency. Validation is performed on three benchmark cases: a stationary monopole, a stationary dipole (both in uniform flow with angle of incidence), and the 30P30N high-lift airfoil. For the monopole and dipole, numerical results match analytical solutions closely. For the 30P30N airfoil at Re = 1.71e6 and 5.5° AoA, far-field noise is predicted. The flow field is computed using improved delayed detached-eddy simulation (IDDES) in OpenCFD-EC. Acoustic predictions agree well with experimental data. To handle large-scale simulations, a hybrid MPI/OpenMP parallelization is adopted. Speedup reaches 538.5 times in the 30P30N case and overcomes memory constraints. The code is released publicly at https://github.com/Z-K-L/OpenCFD-FWH.

ABSTRACT. Combat aircraft typically carry weapons within concealed bays; however, during weapon release, unsteady acoustic–structure interactions dominate over other aerodynamic phenomena, and this effect becomes increasingly severe with increasing cruise speed. Despite significant advances in understanding the influence of Mach number on cavity acoustics, several key issues remain unresolved, particularly the prediction of acoustic levels, their statistical origin, and the scaling behaviour of both oscillation frequency and sound pressure level. A transitional-open type of cavity with length to depth (L/D) ratio of 6.25 and width to depth (W/D) ratio of 2 has been chosen to gain a deeper understanding of the acoustic noise characteristic. To achieve a trade-off between the computational cost and accurate representation, Detached Eddy Simulations are performed. Initially, the simulation methodology was established in a commercial CFD solver, HiFUN, and validated against an experimental study. From this standpoint, a further analysis is performed to investigate the dependence of the cavity acoustic field on Mach number. A systematic increase of dB level is observed with increasing Mach number while the Strouhal number decreases slightly.

ABSTRACT. Accurate prediction of jet noise remains a major challenge in computational aeroacoustics due to the high computational cost associated with resolving both the jet near- and far-field. In the present work, we perform large-eddy simulations (LES) for subsonic (Mj = 0.9) and supersonic (Mj = 1.5) isothermal round jets at high Reynolds number at inflow, using high-order compact finite difference schemes. The acoustic propagation in the far-field is calculated by solving isentropically linearized Euler equations (ILEE) using the same high-order compact finite difference schemes as the near-field LES. Unlike the conventional one-way coupling approach, in which near-field flow information is passed to the far-field solver without feedback, in the present two-way coupling approach, far-field results are taken into account when calculating near-field variables. Results show smooth acoustic wave propagation across the coupling interface. The predicted potential core lengths, jet spreading rates, and axial variation of overall sound pressure level are consistent with trends reported in earlier experimental and numerical studies. The results demonstrate that two-way coupling is a robust and physically consistent approach for direct jet noise computation for compressible round jets.

ABSTRACT. This study investigates whether an algebraic Volume-of-Fluid (VOF) method remains stable when coupled with a low-dissipative high-order convective discretization for incompressible air–water free-surface flows. Breaking waves, characterized by strong nonlinear transport and severe topological transitions, are employed as stringent validation benchmarks to assess the robustness of the algebraic VOF framework without geometric reconstruction. The incompressible Navier–Stokes equations are solved using a one-fluid formulation within an in-house projection-based solver. Nonlinear convective terms are discretized using a WENO-based high-order scheme, while the interface is captured using an algebraic VOF approach with interface sharpening. Results demonstrate stable interface evolution throughout crest steepening, jet formation, plunging impact, and air entrainment. The normalized wave energy history shows close agreement with reference data, indicating that excessive artificial dissipation is not introduced despite the reduced numerical dissipation of the high-order scheme. Extended three-dimensional simulations further confirm robustness under both spilling and plunging conditions. These results clarify the stability characteristics of algebraic VOF under high-order low-dissipation transport and provide insight into its applicability to severely nonlinear free-surface flows.

ABSTRACT. Eulerian dispersed-phase models are widely used in industrial computational fluid dynamics because they avoid Lagrangian particle tracking and integrate efficiently into implicit finite-volume solvers. However, many dilute dispersed-phase formulations are effectively pressureless, leading to a weakly hyperbolic momentum subsystem. This degeneracy can produce nonphysical concentration spikes and reduced robustness of Godunov-type fluxes. In addition, fixed-diameter closures do not capture the feedback of breakup and coalescence on interphase momentum and heat transfer without resorting to computationally expensive population balance models.

This work presents a pressureless Eulerian dispersed-phase framework augmented by a one-group interfacial area transport equation (IATE) and a packing-activated friction-pressure closure. The IATE provides computationally efficient mean-size evolution, with the particle diameter obtained consistently from the transported interfacial area concentration. Hyperbolicity is restored in dense regions by introducing a friction pressure that activates near packing and yields a bounded characteristic speed. The resulting subsystem admits finite wave speeds along face normals, enabling robust Rusanov, AUSM+up, and HLLC flux formulations while preserving the pressureless limit in dilute regions. A receiver-side packing limiter ensures conservative transport and bounded volume fraction.

All fluxes and source terms are linearized using forward-mode automatic differentiation, allowing fully implicit Newton–Krylov solution strategies. The model is extensively tested and is actively used in coupled multi-dimensional solid rocket motor simulations, where strong acceleration and local packing effects challenge conventional pressureless formulations.

ABSTRACT. In order to accurately simulate compressible two-phase turbulent flows, non-dissipative properties in convective numerical scheme are essential to resolve fine vortex structures without nonphysically smearing by numerical dissipation. However, the simulations of compressible two-phase flows with non-dissipative numerical schemes remain a major challenge, as steep interface gradients and high density ratio trigger numerical instabilities in simulations of two-phase flows. To overcome this challenge, this paper proposes a novel non-dissipative and robust numerical scheme for compressible two-phase flows. We derive the kinetic energy and entropy preserving (KEEP) scheme, which exactly conserves the entropy, for compressible two-phase flows based on the seven-equation model. Our proposed scheme achieves both kinetic energy and entropy conservation at the discrete level and enables robust and non-dissipative simulations of compressible two-phase flows even under high-density ratio conditions, something that existing numerical schemes fail to do robustly. Furthermore, to the best of our knowledge, this is the first scheme to strictly conserve entropy for compressible two-phase flows. The theoretical analysis of entropy conservation error leads to the formulations for half-point numerical fluxes that strictly satisfy entropy conservation. In numerical tests, the entropy conservation property and the robustness of the proposed scheme in two-phase turbulent flows with high density ratios are verified.

ABSTRACT. The pressure wave propagation through confined environments scattered with obstacles occurs in various fields of application: for instance the studies of the impact of ventilation arrangements during the accidental explosion of a hydrogen cell inside a ship, or the mitigation of micro pressure waves associated with the entry of a high-speed train into a tunnel using retaining blocks structures. The present work is focused on the case of loss-of-coolant accident (LOCA) where a rarefaction wave propagates through the vessel of a pressurised water reactor (PWR). The complexity and scale of the entire primary loop of a pressurised water reactor (PWR) in the event of a loss-of-coolant accident (LOCA) make it impossible to use a fine mesh to represent the geometry of the obstacles. The methodology presented relies on a coupled 1D/3D representation where geometric details, such as perforated plates in the vicinity of the reactor core or diaphragms positioned in pipes, are modeled through carefully calibrated local impedance relations. The MADMAX experimental setup allows to perform the measurements needed for such a calibration process and also provides reference data for code validation. The impedance relations used in the efficient, yet accurate, simplified 1D approach are implemented in the Finite Element framework of the EuroPlexus code as an added mass term. The objective of the present study is to develop a similar simplified 1D approach in the Finite Volume framework also available within the EuroPlexus software.

ABSTRACT. This study investigates the energy partitioning during the collapse and rebound of a spark-generated cavitation bubble, explicitly accounting for mechanical, thermal, phase change, and shockwave energy transformations. For the first time, the phase change energy associated with condensation and evaporation, and the resulting mass transfer, is directly computed using a comprehensive cavitation model, with each dominant energy component evaluated individually rather than inferred as a residual. The temporal evolution of potential, kinetic, internal, phase change, and shockwave energies is tracked throughout the bubble dynamics. Although the system is not perfectly conservative due to numerical dissipation and far-field losses, the primary energy pathways exhibit a consistent overall balance. Near-wall bubble collapse is shown to suppress shockwave emission while enhancing kinetic energy dissipation, primarily through jet formation, in agreement with experimental observations. These results establish a validated framework for energy partitioning and provide new insights into the coupled mechanical, thermal, and phase change processes governing cavitating flows and underwater explosions.

ABSTRACT. A spacecraft atmospheric re-entry engages considerable complex phenomena regarding the surrounding fluid, the entering solid, and the aerothermal coupling between the two. During its hypersonic re-entry, the heat shield experiences extreme conditions, due to very high heat flux, and will be greatly damaged, in order to protect the vehicle's interior [1]. By sacrificing itself, the thermal protection system lowers the temperature through in-depth degradation, due to pyrolysis chemical decomposition and through surface degradation, due to ablation. Furthermore, the thermal protection system undergoes deformations, such as thermal expansion, swelling or shrinkage, throughout its heating. As a result of numerous sudden thermal phenomena, the domain undergoes significant shape variations. The aim of this work is to developed a robust numerical formulation accounting for the two different domain variations caused, on the one hand, by ablation, which will be addressed through an Eulerian mesh motion method, and, on the other hand, by deformations, which will be handled through a Lagrangian mesh motion. Therefore, an ALE-formulation is stated. After introducing the ablation-pyrolysis-deformation model employed, numerical schemes preserving mass and energy conservation are presented, using a finite volume formulation. Moreover, since preserving physical properties is required, implementation steps and the Newton resolution method are presented. Also, alongside the two previous mesh motions, a mesh adaptation process is added in order to enhance accuracy, by tracking the pyrolysis front and the temperature gradients. The mesh motion studied in this work is built on the spring analogy and preserves the number of nodes in the domain. In the following section, the governing equations of the thermal and mechanical response are described, stating the assumptions, outlining the boundary conditions and focusing on conservative equations. Subsequently, the discretization is presented. In the next section, the various mesh movements are discussed, focusing on the methods used and their limitations. In the last section, the performance of the numerical methods is analyzed on an ablation test case [2] and validated on CEA experimental arc-jet cases.

The study considers an ablative charring solid during hypersonic re-entry, which may endure deformations. In the first part of this work, the full thermal and mechanical response is established, focusing on mass and energy conservation. The study domain consists of a pyrolyzable solid and its pyrolysis gases, generated by its decomposition, which may vary in time. First, mass conservation in the domain, composed of the solid and the gas, implies that the solid mass loss, due to pyrolysis, is transferred to the pyrolysis gas mass. Therefore, the mass loss modeling is introduced as a sink term into the solid mass conservation equation and as a source term in the gas mass conservation equation. Moreover, in order to ensure mass conservation in the solid enduring deformation, two newly equations are derived on the virgin density and the char density [3]. Thus, the solid mass conservation equation is stated. Besides, the mass conservation equation of the pyrolysis gas, of a bulk density, is stated using the mass flow rate [1]. Moreover, energy conservation equation is stated using conservative form in the domain [4,5], considering the solid and gas energy equations and assuming thermal equilibrium. Furthermore, the solid is subjected to a surface ablation because of extreme heat flux. Complex interactions arise between the incoming flow, species in the viscous boundary layer, pyrolysis gases and radiative phenomenon. Therefore, ablative boundary conditions [6,7] must be applied to the energy conservation equation. Thus, this ablative boundary condition leads to the computation of wall recession. Finally, the pyrolyzable solid can be subjected to various deformations during the heating process. Thermal expansion, shrinkage or swelling can occur. Thus, deformations are defined by an empirical law derived by Henderson [8], assuming that the domain velocity depends on the temperature variation and the material degradation advancement variation [3].

The numerical finite volume schemes employed to tackle this problem are detailed, describing the equation discretization and the resolution steps of the thermal and mechanical material response model. The time discretization is performed using implicit Euler scheme. Each equations are integrated on cells which may evolve over time, due to solid deformation, ablation or mesh adaptation. Therefore, the Discrete Geometric Conservation Law (DGCL) [9] has to be applied because of the mesh motion, introducing a new term into equations, in order to preserve property conservation. Moreover, the discretization of terms linked to pyrolysis gases is performed implicitly using an upwind scheme in the direction of the pyrolysis gases. Furthermore, the conduction term is implicitly discretized with a conductivity on the edges computed using a harmonic mean, in order to preserve the piecewise linear fields [10].

At each time step, two mesh motions are required to capture ablation and deformations. In addition, mesh adaptation is employed to enhance accuracy. As a result, three mesh movements are computed at each time step, preserving the same number of nodes. First, the ablation mesh motion is driven by the ablation velocity determined by the strong fluid-solid interaction at the wall, which is resolved iteratively. An Eulerian mesh motion method is used. Besides, the mesh adaptation also adopted an Eulerian mesh motion method. The mesh adaptation, tracking the pyrolysis front and the strong temperature gradients, is based on spring theory [11]. Each node is connected to its neighbors by a spring and the boundary nodes are fixed. This method can be applied in 2D and 3D by splitting the resolution in each direction [12]. Refining the mesh at strong temperature and density gradients, and coarsing the mesh where variations are less significant, this method is particularly suitable on problems with continuous field evolution. However, the selection of the spring stiffness is a critical aspect of the method. Unlike the two previous mesh motions, the deformation mesh motion is a Lagrangian movement related to the solid. Moreover, given the characteristic times for the solid thermal response and the mechanical response, deformations are computed in steady-state. Each new node position is computed using the deformation equation. When moving the mesh through Lagrangian deformation, the properties associated with the cells must also be adjusted to ensure consistency and preserve mass and energy conservation. Therefore, a conservative projection is performed [3].

The solid thermal and mechanical response model have been validated on several test cases. First, the second ablation workshop test case [2] is analyzed by checking the thermal response, the mass and energy conservation at the discrete level and the performance improvement achieved by the mesh adaptation. Additional code validation was performed through an experimental campaign performed on CEA arc-jet facility, operated by ArianeGroup, in Issac, France. Arc-jet campaign are conducted to evaluate the thermal protection system and characterize the materials used during atmospheric re-entry. While heat flux and pressure experienced by the re-entry body cannot be fully replicated, arc-jet tests give an effective representation of intermediate flight conditions [13], especially the flow interaction at the wall. Wall displacement was measured by laser scans, analyzing the shape deformation under ablation and swelling. Two incoming flow regimes were studied, corresponding to two phases of an atmospheric re-entry: a pyrolysis regime and an ablation regime. First, the wall displacement experimental data was compared to two simulations (one considering deformations and the other excluding deformations). Without taking swelling into account, the wall displacement is far from the experimental results. Including swelling in the simulation provides a reasonable close approximation of wall recession and improves the prediction of ablation rate.

Detailing how three different mesh movements are handled and analyzing numerical methods are the main focus of this study. First, the physical phenomena and the model are described. Numerical schemes and resolution are summarized focusing on mass and energy conservation. Each mesh movement method is discussed, analyzing the numerical methods and the time step limitations. The validation of the model has been performed using the second ablation test case where the results obtained are very satisfactory when compared to the reference values. Furthermore, the model and the numerical resolution preserve conservation laws. Applying these methods on in-house arc-jet tests, a close approximation of wall displacement during hypersonic re-entry simulations is achieved. Furthermore, the solid deformation has a slight effect on the temperature evolution of near-wall thermocouples. To further emphasize the importance of resolving the coupling of the three mesh motions, simulations must be validated using other experimental data, such as other arc-jet experimental campaign and tested on in-flight experiments, using 2D or 3D simulations.

[1] G. Duffa. Ablative thermal protection systems modeling. American Institute of Aeronautics and Astronautics, Inc., 2013. [2] J. Lachaud, A. Martin, T. van Eekelen, and I. Cozmuta. Ablation test-case series #2. numerical simulation of ablative-material response: Code and model comparisons. 2012. [3] A. Cas, C. Baranger, H Beaugendre, and S. Peluchon. Conservative models and numerical methods for pyrolysis-thermal coupling of heat shield degradation and deformations. International Journal of Heat and Mass Transfer, 256:127962, 2026. [4] M. Howard and B. Blackwell. A multi-dimensional finite element based solver for decomposing and non-decomposing thermal protection systems. In 45th AIAA Thermophysics Conference, page 2506, 2015. [5] A. Amar, N. Calvert, and B. Kirk. Development and verification of the charring ablating thermal protection implicit system solver. In 49th AIAA aerospace sciences meeting including the new horizons forum and aerospace exposition, page 144, 2011. [6] C.B. Moyer and R.A. Rindal. An analysis of the coupled chemically reacting boundary layer and charring ablator. part II -finite difference solution for the in-depth response of charring materials considering surface chemical and energy balances. Technical report, NASA, 1968. [7] J. B.E. Meurisse, J. Lachaud, F. Panerai, C. Tang, and N.N. Mansour. Multidimensional material response simulations of a full-scale tiled ablative heatshield. Aerospace Science and Technology, 76:497–511, 2018. [8] J.B. Henderson and T.E. Wiecek. A mathematical model to predict the thermal response of decomposing, expanding polymer composites. Journal of composite materials, 21(4):373–393, 1987. [9] C. Farhat, P. Geuzaine, and C. Grandmont. The discrete geometric conservation law and the nonlinear stability of ale schemes for the solution of flow problems on moving grids. Journal of Computational Physics, 174(2):669–694, 2001. [10] G. Chanteperdrix. Modélisation et simulation numérique d’écoulements diphasiques à interface libre. Application à l’étude des mouvements de liquides dans les réservoirs de véhicules spatiaux. PhD thesis, École nationale supérieur de l’aéronautique et de l’espace, 2004. [11] F.J. Blom. Considerations on the spring analogy. International journal for numerical methods in fluids, 32(6):647–668, 2000. [12] K. Nakahashi and G.S. Deiwert. Three-dimensional adaptive grid method. AIAA journal, 24(6):948–954, 1986. [13] A. Balter-Peterson, F. Nichols, B. Mifsud, and W. Love. Arc jet testing in NASA AMES research center thermophysics facilities. In AlAA 4th International Aerospace Planes Conference, page 5041, 1992.

ABSTRACT. This work presents an interpolation-free strategy for handling grid connectivity changes in multi-step Arbitrary Lagrangian–Eulerian (ALE) simulations of inviscid flows. In unsteady adaptive CFD computations, mesh deformation is typically combined with periodic remeshing to maintain grid quality and accurately capture evolving flow features. However, connectivity changes require transferring the solution from the old mesh to the new one, usually through interpolation. Such interpolations may degrade conservation properties, compromise stability, and complicate the use of multi-step time integrators, which require consistent data at several previous time levels.

The proposed methodology embeds mesh adaptation directly within the ALE formulation by representing local connectivity changes as a sequence of fictitious continuous deformations of the control volumes. Each topological modification is decomposed into collapse and expansion phases, allowing the volume swept by the interfaces to be computed geometrically. These contributions are accumulated to define interface velocities that satisfy the discrete Geometric Conservation Law (GCL) by construction. As a result, grid motion and connectivity updates are incorporated without any explicit interpolation of the solution.

To address the challenges posed by time-dependent index sets in multi-step schemes, the approach introduces ghost entities, enabling a consistent treatment of elements that disappear or are created during adaptation. This cumulative framework allows straightforward integration with backward differentiation formula (BDF) time schemes while preserving conservation and stability.

Numerical results for inviscid test cases demonstrate that the method accurately reproduces steady and unsteady solutions, maintains mass conservation under significant mesh adaptation, and avoids the complexities associated with solution transfer between meshes. The proposed strategy therefore offers a robust and efficient tool for adaptive unsteady CFD simulations based on multi-step ALE formulations.

ABSTRACT. Hypersonic flows are very challenging for numerical methods. In this context, the mesh is as important as the numerical scheme to achieve accurate simulations. It has been shown that Voronoi meshes tend to dampen shock instabilities such as the carbuncle effect, thanks to the greater rotational invariance that their cells provide in numerical dissipation compared to their simplicial counterparts. Furthermore, a specific kind of Voronoi meshes, namely Centroidal Voronoi Tessellations (CVTs), naturally arise as minimisers of the error of cell-centered Finite-Volume (FV) schemes. Several efficient algorithms have been developed to generate triangulations and are implemented in various meshing tools. Leveraging the geometrical duality between Voronoi tessellations and Delaunay triangulations, one can construct the Voronoi diagram from a triangulation thus rely on the same efficient core algorithm. Mesh adaptation has proven to be an efficient tool to increase the accuracy of numerical simulations while keeping a moderate computational cost. Our aim is thus to develop a methodology to construct and adapt Voronoi meshes according to stiff compressible effects that occur when simulating hypersonic flows.

Considering a second-order accurate cell-centered FV scheme, the pointwise truncation error that stems from the reconstruction of a given C2-smooth scalar field can be estimated from its Taylor expansion. This estimation depends on the Hessian of the scalar field that is considered. Integrating over all the cells of a mesh, one can get the total error estimate in Lp norm. Now the mesh minimising this error is anisotropic and is described by a Riemannian metric field. For simplicity of mesh construction and robustness of the numerical scheme, an isotropic approximation can be derived by replacing the metric with an isotropic approximation, i.e. a multiple of the identity matrix. This approximation can be chosen to be in some sense close to the original metric by minimising their distance under a given matrix norm. Solving the problem under the Frobenius norm, one gets that the scalar field can be chosen as the arithmetic mean of the eigenvalues of the Hessian matrix. Note that the requirement on the scalar field being C2-smooth is only theorical as constructing the Hessian matrix estimation of a discontinuous field will lead to an implicit smoothing of the mesh density near discontinuities, leading to smooth mesh gradation which is beneficial for numerical schemes. In fact, after the construction of the mesh density field from the Hessian approximation, a smoothing step has been added to control the variation of the field over the whole computational domain.

One would want to build the mesh that minimises the isotropic error estimation. This is a well-known optimal quantisation problem for which local minima belongs to a class of Voronoi diagram called stable CVTs, which can be obtained by gradient descent. Using an asymptotic result which relates the mesh density field of a CVT to the size of the cells as their number tends to infinity, an asymptotically equivalent error estimator can be derived and makes a direct link with the concept of CVT energy.

In the Euclidean plane, the Voronoi diagram of a discrete set of points called generators is a tessellation which is defined as the union of the set of cells, each defined as the set of all points that are closer to their associated generator than to all the others. This diagram has the property of being the dual of the Delaunay triangulation, widely used to generate simplicial meshes for numerical computations. A particular kind of Voronoi diagram is obtained when the centroid of each cell lies exactly on its corresponding generator. This type of diagram is called a Centroidal Voronoi Tessellation and has several mathematical properties which makes it a proper mesh for the FV method.

The Voronoi diagram is constructed from a Delaunay triangulation by leveraging the duality properties that links the two objects. Given a set of generators, their Delaunay triangulation is first constructed, then the corresponding Voronoi diagram is built upon it. The Delaunay triangulation is generated from an efficient implementation of Bowyer-Watson algorithm. The Voronoi diagram is unbounded and needs to be clipped with the geometry in order to form a proper mesh. For this task, we implemented a clipping strategy that relies on Greiner-Hormann polygon clipping algorithm.

The meshing tool is coupled with a hypersonic flow solver in an iterative loop. This loop is initialised from a uniform mesh on which a first flow computation is performed. From this first computation, the value of the Hessian matrix of the pressure field is approximated on each cell. This approximation serves for computing the mesh density field to create a new mesh. For computational efficiency, primitive variables are interpolated from one mesh to the next so that the hypersonic computation is initialised with a good guess on each new mesh. For the same reason, each adaptated mesh is constructed from the previous mesh. The adaptation loop runs until convergence is reached, meaning that the mesh does not change drastically between two adaptation loops. The methodology is applied to an inviscid hypersonic flow at Mach 8.15 around the double ellipse geometry. The incident flow at angle 30° creates a bow shock. This result has been computed on the adapted mesh obtained as the output of the adaptation loop, performed with a fixed number of cells N=20,000 from an initial uniform mesh. The choice of a specific Lp norm affects the number of cells that are created near steeper shocks. The adaptation is therefore done in L1 norm to keep a reasonable amount of cells near wall as well as to resolve adequatly the second shock due to the bump in the geometry. The pressure coefficient obtained with the adapted mesh is compared to the pressure coefficient obtained with the uniform mesh in order to assess the efficiency of the mesh adaptation process. The curves are compared to a reference FV computation by P. Vankeirsbilck and H. Deconinck. The pressure coefficient obtained on the adaptated mesh matches the reference much more closely over the whole geometry than the one obtained on a uniform mesh. The interest of the mesh adaptation process is especially highlighted when looking at the pressure coefficient around the stagnation point and near the geometrical bump where the amplitude of the variations is systematically underestimated on the uniform mesh.

Assuming that the largest error of the numerical scheme comes from the truncation error of the reconstruction process, this work casts the problem of finding the mesh minimising the error of a FV scheme into the optimal quantisation framework. While the error estimation is derived in an analogous manner as the usual Finite-Element-oriented one, the cell shapes are not assumed a priori. This shows that from this perspective, a CVT is actually the best mesh that one can create given a fixed cell budget. This allows for a drastical increase on the quality of the results while keeping the same number of cells in the mesh. The specific CVT to construct is dictated by the choice of the error norm and depends on the tradeoff that the user wants between capturing the sharper features accurately and capturing several flow features adequately. This choice is all the more important when considering discontinuous fields, i.e. for hypersonic flows.

ABSTRACT. The numerical simulation of incompressible flows remains challenging in computational fluid dynamics due to the wide range of spatial and temporal scales encountered in practical applications. Standard discretization strategies often require fine spatial meshes and small time steps to accurately capture relevant flow features, resulting in a high computational cost. To address this issue, several modeling and numerical approaches have been proposed, including Variational Multiscale (VMS) methods, reduced order methods and mesh adaptation techniques.

In this work, we propose a space-time adaptive strategy for the numerical approximation of incompressible flows, based on recovery-based a posteriori error estimators driving anisotropic mesh adaptation in space and adaptive time-step selection. The proposed approach increases the spatial mesh resolution in regions of complex flow behavior while allowing for coarsening elsewhere. In addition, the time step is adaptively reduced when the flow dynamics evolve rapidly and increased when the flow features vary slowly.

The developed space-time adaptive algorithm is validated on several benchmark test cases, including the fluidic oscillator. The numerical results show significant computational gains, with time savings of up to 75% for the unsteady fluidic oscillator and even larger improvements for stationary flow configurations. Overall, this work represents a first step toward fully space-time adaptive strategies in fluid dynamics and lays the groundwork for extending the proposed approach to fully three-dimensional adaptive mesh techniques.

ABSTRACT. This work presents a predictive mesh adaptation strategy for unsteady CFD simulations, aimed at reducing the frequency of remeshing operations while maintaining accuracy. In time-dependent flows, key structures such as shocks and discontinuities move across the domain, requiring dynamic grid refinement to capture steep gradients accurately. Traditional adaptive methods update the mesh frequently to track these evolving features, which increases computational cost. The proposed approach instead predicts the evolution of the target grid spacing (metric field), allowing mesh adaptation to anticipate feature motion over multiple time steps.

The method starts from a computed flow solution, obtained here by solving the compressible Euler equations with a node-centered finite-volume formulation, although the framework is general and not restricted to this setup. A metric field prescribing the desired element size is constructed. Diverse techniques are available in the literature: as an example, in this work, the Projected Hessian-based Mesh Adaptation method (PHMA) is applied to the solution. To predict the metric evolution, the target grid spacing is first interpolated from the unstructured CFD mesh onto a structured Cartesian grid. This enables the use of template matching algorithms from Computer Vision. By comparing metric “images” at two different time levels, a motion field is estimated by minimizing the Sum of Squared Differences. The resulting displacement vectors represent the advection velocity of flow features between the two time steps. Assuming locally constant velocity over short time intervals, the target grid spacing is then extrapolated forward in time to predict the metric for future steps. The predicted metric is finally interpolated back onto the original unstructured mesh and used to guide adaptation.

Preliminary results on a two-dimensional cylindrical converging shock test case show that the predicted motion closely matches the analytical propagation speeds of the main flow features. The predicted metric field agrees well with the exact one, demonstrating that the approach can effectively anticipate the transport of existing structures. The methodology is currently under development and is being assessed on additional test cases. As expected, the method cannot anticipate the emergence of new flow features. Nevertheless, in simulations where the flow is characterized primarily by the transport of already developed structures, the proposed strategy can significantly reduce the frequency of remeshing operations, thereby reducing the overall computational cost of unsteady adaptive simulations.

ABSTRACT. 1. Introduction Due to the intense aerodynamic heating encountered under high–Mach-number flight conditions, the thermal protection system (TPS) has become one of the most critical subsystems in vehicle design. In current engineering practice, ablative thermal protection is widely adopted owing to its robustness and high heat-load capability, making it one of the predominant TPS solutions. Unsteady Computational Fluid Dynamics (CFD) simulations of in-flight ablation can characterize, in considerable detail, the evolution of the material temperature field and material loss. However, high-fidelity CFD modeling that rigorously accounts for thermo-chemical nonequilibrium in the flow as well as surface ablation recession typically incurs prohibitive computational cost. This burden is particularly severe in trajectory analyses, where the flow field, ablation, and thermal response must be solved repeatedly across many flight states, making the total cost of a single flight trajectory exceedingly large. In recent years, using machine learning (ML) to replace part of the CFD workload and rapidly predict aerothermal quantities to reduce overall cost has become an active research focus. However, ML studies dedicated to ablation remain relatively scarce; even within the broader aerothermal-ML literature, the emphasis has largely been on generalization across freestream conditions and geometry parameters. Existing datasets are often generated under isothermal-wall boundary conditions, with wall temperatures typically fixed at 300K or near-ambient values [1]. In contrast, the ablation mass flux and the associated aerothermal environment are highly sensitive to wall temperature and can vary markedly with wall-temperature changes. Nevertheless, systematic modeling that treats wall temperature as a key conditioning variable remains limited, and studies that explicitly target material ablation mass flux as a learning objective are even rarer. To address these gaps, this work develops a machine-learning surrogate to replace CFD computations, aiming to rapidly predict the wall-distributed aerothermal heating and ablation mass flux under varying freestream conditions and wall-temperature distributions for carbon-based ablative materials. We design a multi-layer convolutional neural network to capture the nonlocal influence of wall information on local quantities, and employ feature-wise linear modulation (FiLM) to inject freestream parameters into the network in a conditional manner, thereby improving cross-condition generalization. The training dataset is generated using an in-house coupled aerothermo–ablation–thermal-response simulation platform. The approach is evaluated on a representative flight trajectory reported in the literature. Across 21 trajectory points whose freestream conditions and wall-temperature distributions are unseen during training, the model achieves a mean relative error of 3.63% for both aerothermal heating and ablation mass flux in the stagnation-region vicinity. These results demonstrate that the proposed surrogate generalizes well with respect to both freestream conditions and wall-temperature distributions, providing an efficient predictive capability for TPS design and trajectory-level rapid assessment. 2. Numerical methods and datasets In the CFD platform employed in this study, the flow fields are governed by Navier–Stokes equations with thermo-chemical nonequilibrium effects, while the solid domain is modeled by the heat-conduction equation. Fluid–solid coupling is implemented using a partitioned, loosely coupled strategy: the fluid and solid domains are solved with separate governing equations and solvers, and information is exchanged at the fluid–solid interface. At each prescribed time instant, the flow field is solved in a steady manner, whereas the solid domain is advanced unsteadily in time. In addition, after each flow field solution, the computational mesh is deformed to account for geometry changes induced by ablation recession.This study employs the Park model [2] to describe surface oxidation reactions and uses the Knudsen–Langmuir formulation to compute the sublimation mass flux. On the flow field side, an 11-species, 24-reaction chemical-kinetics model including carbon-related species is adopted [3]. During fluid–solid coupling, mass and energy conservation at the coupling interface are enforced. The net ablation mass flux produced by all surface reactions is defined as one output channel of the machine-learning model; for the energy balance, the net sum of the remaining energy terms—excluding the heat flux conducted into the solid and the radiative heat flux—is defined as the other output channel. Based on the coupled CFD framework described above, a representative blunt-body configuration reported in the literature is selected as the test case, with the corresponding geometry and solid material properties adopted [4]. Coupled simulations are conducted at several fixed freestream conditions while the solid progressively heats up and undergoes ablation, thereby forming the training dataset. The freestream cases cover 65km (Ma=23,25), 60km (Ma=20,25), 50km (Ma=15,20,25), 40km (Ma=15,20,25), and 45km (Ma=22). The wall-temperature range in the training dataset spans 300K–3800K. 3. Neural network model and prediction results To address the fact that the target quantities may span multiple orders of magnitude due to variations in freestream conditions or wall-temperature distributions, we apply a logarithmic transform to both the input and output channels. In addition, the predicted quantity is decomposed into the sum of three components: (i) a baseline term determined solely by the freestream conditions to establish the overall magnitude; (ii) a temperature-modulation term induced by the wall-temperature distribution to further adjust the magnitude; and (iii) a shape term that describes the along-surface spatial distribution at the prescribed magnitude. To ensure that these components fulfill their intended roles, the baseline term depends only on the freestream conditions, while the shape term is stabilized via a self-anchoring treatment; the temperature-modulation and shape terms are both generated by the same six-layer convolutional neural network. Within this convolutional network, feature-wise linear modulation (FiLM) is introduced to condition the internal features on the freestream parameters, enabling a deeper representation of freestream effects and improving cross-condition generalization. A representative flight trajectory from the literature [5] is selected to validate the proposed machine-learning prediction method. Along this trajectory, unsteady coupled CFD simulations accounting for ablation recession are performed at 21 trajectory points, covering freestream conditions from 66km down to 39km in altitude and Ma=20–22. For all 21 points, the combinations of freestream parameters and wall-temperature distributions are unseen in the training dataset. The stagnation-point temperature on the vehicle surface increases from 300K to above 3500K. The freestream parameters and wall-temperature distribution at each trajectory point are provided as inputs to the neural network to predict wall aerothermal heating and ablation mass flux, and the predictions are compared pointwise against the CFD results. Focusing on the engineering-relevant regions of high heat flux and high ablation mass flux, the mean relative error aggregated over the two output channels has a mean value of 3.63% and a median of 3.33% across the 21 trajectory points, indicating that most cases fall within an error band of approximately 3%–4%; the 90th percentile of the error distribution is 5.18%. Among the 21 points, the minimum-error case corresponds to 47.3km and Ma=20.69, with a mean two-channel error of 2.38% in the stagnation-region vicinity. The maximum-error case corresponds to 39.4km and Ma=21.35; this altitude is a mildly extrapolative condition relative to the training set, and the mean two-channel error near the stagnation region is 6.27%. Overall, the model achieves high accuracy for interpolative conditions within the training coverage and maintains acceptable errors for mildly extrapolative conditions, demonstrating good generalization with respect to both freestream conditions and wall-temperature distributions. 4. Conclusions Based on coupled CFD simulations of carbon-based material ablation, this study develops a FiLM-CNN surrogate model that treats wall temperature as a key conditioning variable, enabling rapid prediction of wall-distributed aerothermal heating and net ablation mass flux. Validation on a representative flight trajectory with unseen freestream conditions yields a mean two-channel error of 3.63% in the stagnation-region vicinity, with a 90th-percentile error of 5.18%. The model maintains acceptable accuracy for both interpolative and mildly extrapolative conditions, demonstrating good generalization with respect to freestream conditions and wall-temperature distributions. The proposed approach has the potential to substantially improve the efficiency of rapid TPS design and engineering-level trajectory assessment. References [1] L. He, F. Chen, and Y. Qin, “A novel data-driven method for predicting heat flux of hypersonic aircraft based on Fourier neural operator,” Aerospace Science and Technology, vol. 169, p. 111397, 2026, doi: 10.1016/j.ast.2025.111397. [2] C. Park, R. L. Jaffe, H. Partridge, Chemical-kinetic parameters of hyperbolic earth entry, Journal of Thermophysics and Heat Transfer 15 (1) (2001) 76–90. doi:10.2514/2.6582. [3] J. A. McQuaid, A. L. Zibitsker, A. Martin, C. Brehm, Simulation of Graphite Ablation using an Overset Near Body Solver on an Adaptive Block-Structured Cartesian Off-Body Grid, 2022. doi:10.2514/6.2022-4088. [4] Zibitsker A L , Mcquaid J A , Stern E C ,et al.Finite-rate and equilibrium study of graphite ablation under arc-jet conditions[J].Computers & Fluids, 2023:267.DOI:10.1016/j.compfluid.2023.106069. [5] D. W. Kuntz, B. Hassan, and D. L. Potter, “Predictions of ablating hypersonic vehicles using an iterative coupled fluid/thermal approach,” Journal of Thermophysics and Heat Transfer, vol. 15, no. 2, pp. 129–139, 2001, doi: 10.2514/2.6594.

ABSTRACT. Introduction Combined-cycle propulsion systems, particularly Rocket-Based Combined Cycle (RBCC) engines, are considered key technologies for reusable space transportation and high-speed atmospheric flight [1]. During engine operation, intense gas compression and combustion processes subject engine wall structures to extreme thermal loads. For instance, at Mach 8, the total inlet temperature approaches 2700K, while the total exhaust gas temperature reaches 3000K [2]—far exceeding the thermal limits of current material systems. To maintain structural integrity and reliable performance, regenerative cooling using hydrocarbon fuels is widely adopted in RBCC engines. In such systems, fuel absorbs heat through internal cooling channels before being injected into the combustion chamber, serving concurrently as both coolant and propellant. To maximize cooling effectiveness and fuel preheating efficiency, hydrocarbon fuels in RBCC regenerative cooling channels typically operate under supercritical pressure conditions (3–6 MPa) and elevated temperatures exceeding 1000 K. Under these conditions, fluid behavior departs significantly from ideal-gas assumptions. Supercritical fluids exhibit strong thermodynamic non-idealities, including sharp density gradients, specific heat peaks near pseudo-critical regions, and transport property anomalies. In addition, endothermic pyrolysis reactions modify mixture composition and introduce further nonlinear coupling between thermodynamic and transport properties [3]. These complex real-fluid and chemical effects strongly influence local heat transfer coefficients, wall temperature distribution, and flow stability within cooling passages. Accurate numerical simulation of these phenomena requires the use of real gas equations of state (EOS). Models such as the Pen-Robinson (PR) equation or Extended Corresponding State (ECS) are commonly employed to capture supercritical behavior. However, these highly nonlinear EOS models are computationally expensive, particularly in multicomponent reaction systems. Previous studies indicate that in Large Eddy Simulation (LES) of supercritical reaction flows, calculating the thermodynamic and transport properties of real fluids can account for up to 50% of the total computational time [4]. CFD simulations of regenerative cooling channels at the RBCC propulsion-scale require tens of thousands of iterations across tens of millions of grid cells. Therefore, high-fidelity modeling of supercritical pyrolysis hydrocarbon fuels has become the primary computational cost in hypersonic propulsion simulations, limiting parametric design studies and multiphysics optimization. To address this challenge, this study developed a deep learning-based thermophysical substitution framework that replaces iterative nonlinear EOS evaluations with data-driven mapping. By integrating the trained model into large-scale CFD simulations, this approach significantly reduces computational costs while maintaining engineering-level accuracy for RBCC regenerative cooling analysis.

Methodology This study employs n-decane (C₁₀H₂₂) as an alternative fuel and utilizes a simplified RP-3 kerosene cracking mechanism to describe the dominant endothermic decomposition behavior under supercritical conditions. The overall reaction can be represented as: C_{10}H_{22}\rightarrow0.1766H_2+0.7104CH_4+0.748C_2H_4+0.6068C_2H_6 +0.4367C_3H_6+0.2001C_3H_8+0.0482C_4H_6 +0.6148C_4H_8+0.1242C_4H_{10}+0.5612C_6H_6 To efficiently approximate the nonlinear thermophysical behavior of supercritical cracking hydrocarbon fuel, a multi-input, multi-output Deep Neural Network (DNN) surrogate model was constructed. The input vector consists of temperature T, pressure P, and cracking degree α, while the output vector simultaneously predicts density ρ, dynamic viscosity μ, specific heat Cp, thermal conductivity k, and specific enthalpy h. The training dataset was generated using an extended generalized corresponding-state equation of state (ECS-EOS) combined with a simplified cracking mechanism, covering temperature ranges of 280–1300 K, pressures of 2.5–6.5 MPa, and cracking degrees of 0–100%. The training process converged smoothly, and the final loss value was 0.001239. The test set prediction performance showed high accuracy, with an average absolute error (MAE) of 2.328, a mean squared error (MSE) of 22.830, and a determination coefficient (R²) of 0.999751, with an average error of less than 3% for each variable, as shown in Table 1 and Figure 1. These results indicate that the model has excellent regression accuracy and high fidelity over the entire thermodynamic range. Table 1: Prediction accuracy of the DNN model on the test dataset Output MAE MSE R2 MAPE density ρ 0.805 2.562 0.999926 0.3708% viscosity μ 1.617 3.364 0.999821 2.9272% specific heat Cp 6.575 50.729 0.994307 0.1402% conductivity k 0.216 1.001 0.996618 0.2839% enthalpy h 2.425 3.689 0.999990 1.5207%

(a) density ρ (b) dynamic viscosity μ

(c) specific heat Cp (d) thermal conductivity k Figure 1: Prediction accuracy of the DNN model on the test dataset After training, the DNN model is converted to ONNX format and deployed using ONNX Runtime in C language to achieve efficient inference and cross-platform deployment. The alternative model is integrated into the CFD solver using user-defined functions (UDFs), replacing the iterative calculation of the real gas state equation. To ensure scalability in large three-dimensional simulations, we implemented a strategy of tensor memory pre-allocation, batch inference, and cell-level caching to minimize runtime overhead and accelerate simulation speed. The proposed approach was first validated against a typical single-tube experimental case, demonstrating excellent agreement between simulation and measured data, shown as Figure 2(a). Subsequently, the framework was applied to Large Eddy Simulation (LES) of an RBCC regenerative cooling flat-plate channel under supercritical conditions, shown as Figure 2(b). By replacing the conventional ECS-based real-gas property evaluation coupled with a detailed chemical mechanism, the DNN-based framework substantially reduces computational cost. Specifically, the average computational time over ten iterations decreases from 675 s using the ECS + detailed mechanism approach to 38 s with the surrogate model. The DNN-accelerated LES simulations exhibit stable convergence behavior and achieve significant computational acceleration while maintaining physical fidelity compared to traditional real-gas thermodynamic evaluation methods.

(a) The temperature chart of a single-tube experiment (b) The temperature chart of an RBCC Regeneration Cooling Plate Figure 2: The simulation results of the DNN model

Conclusions This paper develops and validates a deep learning–based thermophysical surrogate framework for efficient simulation of supercritical cracking hydrocarbon fuel in RBCC regenerative cooling channels. Through database construction, neural network training, and propulsion-scale CFD deployment, the following primary conclusions are drawn: (1) The proposed model is constructed upon supercritical thermophysical data generated using an extended generalized corresponding-state equation of state (ECS-EOS), covering a wide operating range of temperature (280–1500 K), pressure (2.5–6.5 MPa), and cracking degree (0–100%). By learning the nonlinear mapping between thermodynamic states and fluid properties within this regime, the surrogate model can effectively replace repetitive and computationally expensive real-gas thermal property evaluations in RBCC regenerative cooling simulations, to meet the requirements for equation of state iterative inversion during CFD computations. (2) The proposed multi-input, multi-output Deep Neural Network (DNN) surrogate model can accurately approximate the nonlinear mapping between thermodynamic state variables (temperature, pressure, and cracking ratio) and key thermophysical properties (density, viscosity, specific heat, thermal conductivity, and enthalpy). The trained model achieves an overall coefficient of determination R2=0.99975, with average prediction errors below 3% and maximum errors under 11.5% across the full supercritical operating range. The model remains stable near pseudo-critical regions, where strong property gradients are typically challenging for both numerical solvers and surrogate approximations. (3) Through the ONNX framework, cross-platform deployment is achievable. Leveraging memory pre-allocation, batch inference, and parallel-compatible implementation strategies, this proxy model significantly reduces runtime overhead associated with thermal-physical property evaluation. In a three-dimensional RBCC regenerative cooling simulation involving approximately 1.3 million grid cells, this AI acceleration framework achieved an overall computational acceleration of approximately 6.6 times. This framework establishes a practical engineering pathway for simulating hypersonic air-breathing propulsion systems. By replacing iterative nonlinear EOS evaluations with Deep Neural Network inference, the method significantly enhances computational efficiency while maintaining physical consistency. This enables high-resolution multiphysics simulations, rapid parameter studies, and scalability to other advanced air-breathing hypersonic engine design processes.

References [1] P. Czysz, M. Little, Rocket Based Combined Cycle engine (RBCC)-A propulsion system for the 21st century, in: Proceedings of the 5th International Aerospace Planes and Hypersonics Technologies Conference, 1993. [2] S.L. Zhang, X. Li, J.Y. Zuo, et al., Research progress on active thermal protection for hypersonic vehicles, Prog. Aerosp. Sci., vol. 119.p. 100646. 2020. [3] R. Jiang, G. Liu, and X. Zhang, Thermal cracking of hydrocarbon aviation fuels in regenerative cooling microchannels, Energy & Fuels, vol. 27, no. 5. pp. 2563–2577, 2013. [4] P. J. Milan, J.-P. Hickey, X. Wang, and V. Yang, “Deep-learning accelerated calculation of real-fluid properties in numerical simulation of complex flowfields,” J. Comput. Phys., vol. 444, p. 110567, Nov. 2021.

ABSTRACT. The application of Deep Reinforcement Learning (DRL) to active control (AFC) of the flow over a backward-facing step (BFS) at moderate Reynolds numbers (Re_h=5000–36000) is presented. Two-dimensional numerical simulations are undertaken with flow predictions for the uncontrolled case validated against existing experimental data. The numerical predictions show good agreement, providing a basis for proceeding with the control study. Flow control was achieved using fluid injection or suction through narrow slots, with a DRL agent trained to optimize jet actuation for increasing mean base pressure. Multiple jet positions were evaluated, and both single- and multi-objective reward functions investigated. The results indicate that jet placement upstream of the step yields the best performance, increasing the base pressure coefficient from approximately -0.2 to -0.1 at Re_h=5000, and still maintaining good control at Re_h=36000 with a similar increase. These results demonstrate the potential of DRL to identify effective control strategies for turbulent separated flows, highlighting the importance of actuator placement and reward design. The present work provides a basis for the further extension of DRL-based AFC methods towards more complex flow configurations and practical applications.

ABSTRACT. Thermal phase change phenomena, such as boiling and condensation, are inherently multiphase processes and remain challenging to model due to the strong coupling of heat and mass transfer across phase interfaces. Numerically, the presence of additional source terms in the governing equations strongly influences the stability and convergence of these simulations. Conventional solvers in OpenFOAM, such as interThermalPhaseChangeFoam, assume constant thermodynamic properties and compressibility effects only due to phase change, neglecting the fluid compressibility. These assumptions are introduced to avoid the instabilities due to the non-linear behavior of the phase change process, but often limit the accurate prediction of heat and mass transfer rates in condensation and boiling simulations. To address these limitations and improve the predictive accuracy of heat and mass transfer in compressible boiling and condensing flows, we have extended the current capabilities of compressibleInterFoam to incorporate thermal phase change. Furthermore, the classical assumptions of constant material properties are replaced by a temperature and pressure dependent thermophysical properties in the form of polynomial expressions obtained by fitting the data from CoolProp, an open-source thermodynamic data library. An enthalpy-based formulation of the energy equation is adopted to enhance numerical robustness. However, when thermophysical properties vary with both temperature and pressure, the iterative evaluation of temperature from enthalpy using Newton’s method becomes computationally demanding. To overcome this bottleneck, machine learning (ML) technique is used in modelling this complex sub-process of temperature inversion from enthalpy in the solution procedure. Machine learning techniques in CFD modelling is of growing interest due to its potential to enhance the predictive capabilities of the conventional CFD solvers. Artificial neural network (ANN) models are developed and trained using PyTorch to learn the mapping between enthalpy, pressure and temperature. The ANN models are integrated into OpenFOAM using ONNX Runtime, an open-source, cross-platform inference engine that enables fast and efficient data exchange during simulations. The solver is numerically validated against benchmark cases, Stefan problem and bubble condensation. The results confirm the successful implementation of the framework and demonstrate accurate prediction of phase-change dynamics. This hybrid ML-CFD framework offers a promising pathway in improving accuracy and robustness in condensation and evaporation simulations.

ABSTRACT. We present an integrated framework for uncertainty-aware inverse heat transfer tailored to thermal protection system characterization. The approach targets the recovery of spatially varying thermal conductivity from noisy transient temperature measurements on a one-dimensional slab with mixed boundary conditions, representative of certification-by-analysis workflows for thermal protection materials. To address the cost and fragility of traditional gradient- and adjoint-based calibration in large-scale settings, we employ a physics-informed neural operator with a Fourier neural operator backbone, trained on synthetic ensembles generated from high-fidelity solvers. This learned operator is coupled with Karhunen–Loève-based representations of stochastic input fields to enable efficient uncertainty quantification in the inferred conductivity. A quantum-computing-based Monte Carlo (QCMC) estimator is then used to accelerate sampling-based uncertainty propagation while remaining compatible with near-term quantum hardware. Together, these components deliver near-real-time, uncertainty-aware inverse characterization for thermal protection system–like problems, bridging the gap between high-fidelity Bayesian methods and deployable certification workflows.

ABSTRACT. In the past decade, the numerical simulation of supersonic parachute fluid–structure interaction (FSI) has advanced substantially through the development of Stanford’s AERO Suite and, more recently, NASA’s LAVA framework. Both frameworks have been independently validated against experimental data from the Advanced Supersonic Parachute Inflation Research Experiment (ASPIRE) flight tests, but have not yet been directly compared. In this work, a cross-verification study is performed using the ASPIRE SR03 test article with equivalent physical and numerical settings in both frameworks to evaluate consistency between the frameworks, identify modeling sensitivities, and establish reference benchmarks for future high-fidelity simulations of supersonic parachutes. A sequence of case studies is conducted, including computational fluid dynamics analyses of the payload and payload–parachute systems, structural dynamics analyses of pressure-driven canopy inflation, and fully coupled FSI simulations from line-stretch through full inflation. Preliminary results demonstrate strong agreement between the two frameworks for isolated fluid and structural cases, with fully coupled FSI results for the parachute inflation phase to be presented in the final paper and conference presentation.

ABSTRACT. Driven by the increasing demand for high thrust-to-weight ratios in modern rocket engines, the nozzle design presents two significant trends: large expansion ratios and lightweight structures. However, during low-altitude flight, large expansion ratio nozzles are prone to inducing asymmetric flow separation, which generates unsteady side loads. For lightweight thin-walled nozzles, the coupling effect between the unsteady side loads and the flexible wall is significantly enhanced, which may lead to structural instability or even failure. To reveal the unsteady aeroelastic response mechanisms under low-altitude over-expansion conditions, this study employs a two-way fluid-structure interaction (FSI) method to investigate the transient operation of a thin-walled nozzle with large expansion ratio, focusing on capturing the dynamic feedback mechanisms between fluid and structure. The results indicate that, driven by coupling effects, the separation shock wave exhibits a wave-shaped distortion with a circumferential wave number of 4 during its oscillatory retreat. The spatial distribution of this specific aerodynamic load is consistent with the n = 4 natural mode shape of the structure. Consequently, the divergent section exhibits a square-like deformation consistent with this mode, demonstrating that the shock flow topology is highly sensitive to structural deformation. This study reveals the potential strong coupling aeroelastic characteristics within the free shock separation (FSS) flow field, providing a new theoretical basis for the structural integrity design of thin-walled nozzles with large expansion ratio.

ABSTRACT. Conjugate heat transfer (CHT) in gas-cooled turbine vanes represents a critical multi-physics challenge in turbomachinery design, where accurate prediction of coupled fluid-thermal behavior directly impacts cooling efficiency and component lifetime. Conventional numerical approaches typically address these challenges through partitioned solution strategies, wherein separate solvers for the fluid and solid domains are iteratively coupled via boundary condition exchange. While conceptually straightforward, such loosely coupled methods often impose severe limitations on time-step sizes due to interfacial stability constraints, resulting in slow convergence and computational inefficiency. Alternatively, tightly coupled formulations, which solve the combined system simultaneously, present significant implementation challenges arising from disparate governing equation structures and widely separated characteristic time scales between fluid convection and solid conduction. This paper presents a novel unified high-order CHT solver and its application to efficient simulation of gas-cooled turbine configurations. The methodology employs a double-time-scale nondimensionalization of the compressible Navier–Stokes equations, enabling a single formulation to describe both fluid flow and solid heat conduction domains. A dimensionless time scaling parameter is introduced, which systematically reconciles the inherent temporal disparity between fluid convective and solid conductive transport. This scaling enables the solid domain to evolve at an artificially accelerated rate, thereby dramatically reducing the number of computational steps required to attain a global steady-state solution. To ensure numerical accuracy across material discontinuities, high-order consistent formulations for interface temperature and heat flux are established and enforced. The complete numerical framework is implemented within an existing parallel adaptive high-order discontinuous Galerkin (DG) flow solver, preserving its inherent advantages of high spatial accuracy, geometric flexibility on curved elements, and robust stability at large Courant–Friedrichs–Lewy (CFL) numbers. Following the established validation of the CHT solver [1, 2], this study advances to its application in engineering problems. The capabilities of the proposed unified CHT solver are demonstrated through numerical simulations of representative gas- cooled turbine vane configurations, showcasing its potential for industrial application in turbomachinery thermal analysis and design optimization.

ABSTRACT. Water entry is a critical technical issue in nature, daily life, and engineering fields. In this study, a numerical model for the water entry of large-scale vehicle (with a diameter of 533 mm) was established based on the structured arbitrary Lagrangian-Eulerian (S-ALE) method. The cavity evolution and time-domain characteristics of acceleration during the water entry process were obtained. Furthermore, the frequency-domain characteristics of acceleration were derived by combining modal analysis and shock response spectrum analysis. The relevant research conclusions can provide supplements to the existing numerical simulation results of vehicle water entry.

ABSTRACT. Magnetohydrodynamics (MHD) is a well established description of electrically conducting fluids, with applications in space physics, geophysics, and nuclear fusion. A central challenge in computational MHD modelling is the enforcement of the divergence-free condition for the magnetic field. Although not an explicit transport equation of the MHD system, it appears as a constraint that must always be satisfied given compatible initial and boundary conditions. At the discrete level, however, the solenoidal property may be severely violated if discretization schemes are not constructed with care. Large non-zero magnetic divergence errors can result in unphysical behavior and may compromise numerical stability. Consequently, a variety of strategies have been developed to control divergence errors in computational MHD, including Powell's formulation, projection methods, and constrained transport schemes.

Recently, considerable attention has been directed toward the application of variational data assimilation (DA) to plasma flows governed by ideal MHD descriptions. The authors were the first to apply variational DA to the latter. Data assimilation is the process of systematically combining observational data with a physical model to produce an improved estimate of the system state. In variational DA, this is formulated as a partial-differential-equation-constrained optimization problem, in which a cost functional measuring the misfit between model predictions and observations is minimized subject to the constraints of the governing equations. In the context of solar wind forecasting, for example, MHD-based variational DA can yield significant predictive improvements compared to that of the forward modelling alone.

Analogous to the challenges encountered in forward MHD modelling, the enforcement of the solenoidal condition emerges as a critical problem in DA as well. The distinction, however, is that while forward modelling techniques must treat divergence errors introduced by numerical discretization, DA must additionally handle errors that arise from the data because of sparse spatial coverage, noisy measurements, and issues inherent to the DA algorithm. The present study proposes a novel formulation of ideal MHD-based variational DA that eliminates the introduction of magnetic divergence errors by construction. The central idea is to infer an underlying magnetic vector potential, rather than the magnetic field directly, thereby guaranteeing satisfaction of the divergence-free condition. Details of the proposed vector-potential-based variational DA approach its numerical implementation within a standard upwind finite-volume procedure for the ideal MHD equations will be provided, along with numerical results for several relevant benchmark problems to demonstrate the effectiveness of the approach and its impact on assimilation accuracy and solution robustness.