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. 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. 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. 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. Loss-of-primary-coolant accidents in pressurized water nuclear reactors involve the propagation of a rarefaction wave through an extended flow domain containing numerous geometric singularities (diaphragms, perforated plates, etc.). The efficient numerical simulation of such configurations relies on the simplified description of these singularities using an added mass source term in the governing equations. The present work focuses on the Finite Volume discretization of the modified hyperbolic system, namely the barotropic Euler equations with an added-mass term. Two approximate Riemann solvers (Rusanov and HLL) are developed for the modified system and assessed on a range of test cases, from a 1D shock tube to the reproduction of experiments carried out on the MADMAX platform. Comparisons with a 1D Finite Element discretization of the modified system, a 2D axisymmetric FV simulation with a discretization of the geometric singularities and experimental results allow to establish the Rusanov scheme is unable to provide physically accurate results as soon as the added-mass term varies in space. Meanwhile, the Finite Volume HLL scheme yields consistent results and provides a conservative alternative to the existing Finite Element approach.

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. We present the development, verification, and validation of a numerical algorithm for Lagrangian particle tracking in one-way coupled multiphase flows over grids undergoing arbitrary motion. The carrier flow is assumed to be known a priori and computed on a moving mesh within an Arbitrary Lagrangian-Eulerian framework. Particles are advanced in a non-inertial reference frame attached to the moving grid, while their inertial state is restored at discrete CFD time levels. The formulation relies solely on the discrete mesh motion provided by the carrier flow solver and does not require continuous interpolation of the grid kinematics. The overall framework is verified through analytical tests and validated on an unsteady pitching airfoil configuration.

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 engine 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 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 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. Although both frameworks have been independently validated against data from the Advanced Supersonic Parachute Inflation Research Experiment (ASPIRE) flight tests, meaningful cross-verification of modern parachute FSI frameworks remains challenging due to the large number of interacting modeling choices present in the coupled problem. In this work, a systematic methodology for cross-verifying high-fidelity parachute FSI frameworks is developed and demonstrated using the ASPIRE SR03 test article. The proposed approach decomposes the coupled FSI problem into a hierarchy of structural, fluid, and coupled benchmarks, enabling sensitivities associated with individual modeling assumptions to be examined before considering the complete inflation process. Application of the methodology reveals close agreement in several isolated benchmark cases while also identifying sensitivities associated with structural discretization, fluid mesh representation, wall-boundary treatment, and porous-canopy modeling. Preliminary coupled simulations further demonstrate that modeling choices with modest influence in isolated analyses can become amplified in the fully coupled problem. These results establish a practical foundation for future code-to-code comparisons of parachute FSI simulations and highlight key challenges that must be addressed to achieve meaningful agreement between independent multiphysics frameworks.

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.

ABSTRACT. We report direct numerical simulations of subsonic, transonic and supersonic flows past a circular cylinder. The far upstream Mach number Ma and the diameter-based Reynolds number Re in these three simulations are (Ma = 0.2, Re = 10000), (Ma = 0.9, Re = 3000) and (Ma = 1.2, Re = 10000), respectively. We present tracer particle path-line histogram, conditional histogram and zonal histogram data extracted from the direct numerical simulations. We focus on the characteristics of particle instantaneous flow reversal and their dependencies on Mach number and tracer release position. Such tracer path-line statistics have rarely, if at all, been reported in the literature on compressible subsonic, transonic and supersonic flows past a circular cylinder.

In the subsonic flow, histogram profiles of tracer particles released from a region very close to the wall and slightly behind the mean flow separation location display unexpected persistent oscillations. No similar oscillations are detected from histograms of tracer particles released from any other regions of this subsonic flow. The oscillation frequency is 50 times of the periodic vortex shedding frequency, reflecting the effect of a very small magnitude zig-zagging type of particle motion in a region very close to the wall and slightly behind the mean flow separation location.

No similar high-frequency oscillations are observed in the transonic and the supersonic flow tracer particle histograms. This is likely because, at Ma > 0.8, periodic vortex shedding ceases to exist and the near-wake becomes a quasi-laminar recirculation zone bounded by a pair of converging slip-layers with a neck opening to the far-wake. At the instant of 5 subsonic vortex-shedding periods after being released from the cylinder surface, there are 1.67%, 35.2% and 39.7% particles that experience flow reversal in the subsonic, transonic and supersonic flows, respectively.

Given the dominance of periodic vortex shedding in the subsonic flow and the total absence of it in the transonic and supersonic flows, further considering that the histogram oscillation frequency is a harmonic of the shedding frequency in the subsonic flow, and taking into account of the fact that even in the subsonic flow no histogram oscillations exist for tracers released from either the upstream or the upwind surface, it can be deduced that the newly-discovered tracer zig-zagging motion likely results from the impact by the periodic vortex shedding phenomenon on those particles originating from the cylinder rear surface near the shedding location.

The second interesting feature observed in the present study is that, at Ma = 0.9, the turbulent far-wake behind the pair of tail shockwaves experiences a distinct transverse confinement. This is likely caused by a unique distribution of the transverse velocity field. Further downstream the transitional far-wake gives away to a fully turbulent far-wake. Such inhibited wake growth is absent in both the supersonic flow case and the subsonic case. Hopefully the two flow features newly discovered from the present study can be confirmed independently and investigated further by other researchers.

ABSTRACT. The control rod guide tube is a critical precision component in nuclear reactors, responsible for guiding the stepwise movement of control rod assemblies and ensuring that the drop time meets nuclear safety requirements. However, excessive wear on some control rods has been observed during refueling inspections in operating nuclear power plants. Flow-induced vibration (FIV), arising from the hydraulic forces of the coolant, has been identified as the primary contributor to the observed wear. The periodic fluid excitation is primarily induced by lateral flows in the upper nozzle area, the control rod guide tube flange–core upper plate gap, and at the outlet of the guide tube square slot. To systematically investigate the flow behavior in these critical regions and evaluate optimization design of core upper plate to mitigate FIV, this research adopts a combined approach of experimental flow field characterization (using Pitot tube and PIV techniques) and numerical validation via the commercial CFD software Fluent. The numerical framework for the fluid dynamics simulation employs the Reynolds-Averaged Navier-Stokes (RANS) equations, with the realizable k-epsilon model used for turbulence closure. A 2×2 control rod guide tube model (including control rod guide tube, core upper plate, upper nozzle and so on) is employed to compare the prototype design and an optimized design under different core upper plate and support column base configurations. This study investigates changes in the mass flow rate distributions at the four holes of the core upper plate, as well as the magnitude of lateral flow intensity in upper internal structure (UIS), resulting from the design optimization. Given the highly mature state of nuclear component structural design, such as that of the UIS, this study adopts an optimization strategy based on minimal structural adjustments to improve flow conditions at critical locations and mitigate FIV. The optimization strategy involves adjusting both the dimensions of three key hole types (guide tube, unobstructed flow, and support column holes) in the core upper plate and the design of the support column base. Experimental results demonstrate that the optimized configuration significantly improves both the uniformity and intensity of flow distribution across the core upper plate. The measured flow distribution ratios for the three types of holes (guide tube hole : support column hole : unobstructed flow hole) improved from 1.11:0.82:0.96 in the prototype to 0.97:0.98:1.08 in the optimized model [1]. Numerical simulations further validate these findings, showing a consistent trend, with the ratios changing from 1.17:0.75:0.90 to 0.98:0.96:1.07. Meanwhile, the mutual validation of experimental and simulation data provides consistent evidence that the optimized model reduces flow exchange among the three types of holes, decreasing lateral flow intensity at the upper nozzle region. It also suppresses both the mass flow entering into the guide tube and the jet velocity from the gap between the guide tube flange and the core upper plate. Furthermore, flow velocity and fluctuations are reduced in the continuous guide tube section. Consequently, the pressure drop across the optimized model is only 0.83 times that of the prototype. These improvements in the flow field demonstrate that the optimized model can effectively decrease the excitation forces acting on the control rods, resulting in effective FIV suppression.

ABSTRACT. This paper presents a numerical investigation into the interaction between dispersed particles and turbulent mixing layers induced by cylindrical Richtmyer-Meshkov (RM) instability. Simulations were conducted using OpenCFD-Comb, an open-source code developed by the authors, which adopts a high-order finite difference scheme within a Eulerian-Lagrangian framework. Specifically, we simulated a cylindrical interface perturbed by multimode initial conditions subjected to impact by a converging shock wave with a Mach number of 1.5. The study focuses on the distribution characteristics of particles (at loadings of 50 million and 100 million particles) and their feedback effects on flow evolution. Results indicate that particles preferentially accumulate in spike regions, a phenomenon attributed to centrifugal effects and particle inertia. It is observed that while low particle loading exerts negligible influence on the early-stage flow, higher loading (100 million particles) introduces a macroscopic "lagging effect" on flow compression. Compared to the particle-free case, this higher loading also slightly suppresses the asymptotic growth of the mixing layer width. A Eulerian-Lagrangian framework was employed for modeling the two-phase flow. The fluid phase is governed by the three-dimensional compressible Navier-Stokes equations. For spatial discretization, a sixth-order monotonicity-preserving scheme was utilized for convective terms, whereas an eighth-order central difference scheme was applied to viscous terms. Time advancement was achieved via a third-order Runge-Kutta method. The motion of the dispersed phase was simulated using a Lagrangian tracking approach, with the governing equations for particle position and velocity derived from Newton's second law. The computational domain consists of a Cartesian mesh with a grid resolution of 1024³. A cylindrical interface separates the inner gas (85% H₂ and 15% N₂) from the outer gas (79% N₂ and 21% O₂). A multimode perturbation was imposed on the interface to trigger instability, and the flow was driven by a converging shock wave with an initial Mach number of 1.5

ABSTRACT. The Weighted Essentially Non-Oscillatory (WENO) scheme has attracted much attention due to its ability to simultaneously achieve low dissipation characteristics of smooth regions and high-resolution capture of discontinuities. However, how to simplify the reconstruction process, improve computational efficiency, and enhance the robustness of the algorithm is still an issue in its engineering applications. In this paper, based on the core idea of nested center reconstruction stencils, a simplified multi-resolution WENO (MR-WENO) scheme with linear weights that can be artificially valued and are independent of grid types has been developed. In addition, a new extremum properties (EP)-based discontinuous sensor was designed and a corresponding hybrid MR-WENO scheme was developed to improve it computational efficiency. Numerical results of the inviscid/viscous compressible flow problems show that these schemes have less numerical dissipation than the classical fifth-order WENO-JS scheme.

ABSTRACT. Turbulent flows are intrinsically chaotic and multiscale, making long-horizon prediction both computationally expensive and difficult to keep statistically stable. Neural-operator surrogates such as Fourier neural operator (FNO) can accelerate turbulence prediction, but standard one-step-ahead training followed by autoregressive rollouts often accumulates small errors and causes long-term deviation in key statistics. To assess reliability beyond short-term pointwise accuracy, we propose a unified evaluation framework that integrates (i) uncertainty quantification (UQ) of physically meaningful statistics via error distributions across time resolutions and Fourier modes, (ii) statistical stability diagnostics for long-term rollouts and sensitivity to initial perturbations, and (iii) autocorrelation function (ACF) analysis to connect temporal coherence with predictive robustness. Forced three-dimensional homogeneous isotropic turbulence (HIT) is used as a representative chaotic system. Filtered DNS (fDNS) data are generated from a 256^3 periodic box of size (2pi)^3 and spectrally filtered with kc=10 to obtain resolved velocity fields on a 32^3 grid; the learning task is to approximate the one-step evolution mapping U(t) to U(t+DeltaT) and then deploy the surrogate autoregressively. We compare four FNO-based models against a conventional dynamic Smagorinsky model (DSM). We also consider a constrained prediction setting that enforces the target kinetic energy in the two lowest wavenumber shells (k=1,2) at each predicted step to suppress large-scale deviation. Results show that multiple FNO-based models can achieve superior statistical stability compared with DSM, and factorized-implicit Fourier neural operator (F-IFNO) provides the most favorable trade-off among accuracy, long-term robustness, and computational efficiency. Implicit factorization mitigates error accumulation during autoregressive rollouts, while appropriate time interval selection and large-scale constraints further enhance reliability.

ABSTRACT. 1. Introduction Processes involving powder suspensions are widely used in chemical, pharmaceutical, and food industries, yet scaling up these processes has been very challenging due to strong dependence of hydrodynamics on the evolving suspension properties. As powder is incorporated into the liquid, the density, viscosity and shear response of the suspension change significantly. Consequently, the flow patterns, free surface behavior as well as mixing efficiency change. Traditional scale-up rules such as constant tip speed, constant power per volume or geometric similarities are often insufficient because they do not account for the strong coupling between suspension properties and flow behavior across scales. This work combines experimental characterization and multiphase CFD to develop a robust scale-up methodology for industrial powder suspensions. Instead of attempting to directly model powder addition or wetting, suspension properties were experimentally measured at each incremental solids loading, enabling physically consistent CFD simulations across pilot and production scales. A hydrodynamic similarity criterion – based on matching bulk and surface velocities - was used to determine the impeller speed at production scale. This approach provides a robust alternative to traditional scale-up methods. 2. Methodology 2.1. Experimental Methods A series of powder addition experiments were performed at pilot plant scale under controlled mixing conditions. At every specified mass-fraction increment, the suspension was characterized in terms of density, viscosity and surface behavior (by qualitative assessments). These properties, in particular, the density and viscosity were used as inputs to the CFD simulations. 2.2. CFD methodology Pilot scale mixing was simulated using multiphase CFD model with a Volume of Fluid (VOF) free surface formulation to capture the gas-liquid interface. The governing equations [1] are reported as follows: Continuity equation:

Momentum equation:

Volume fraction equation:

RNG model [2] was used to model turbulence in the system. Simulations were run in ANSYS FLUENT 2025R1 to obtain key hydrodynamic outputs used as scale-up targets. First, we simulated the suspensions at pilot scale to extract key hydrodynamic descriptors: mainly bulk velocity and surface velocity. These descriptors then serve as the basis for hydrodynamic similarity. Then, we simulated the suspension system at the production scale using identical physics and rheology. We then identified conditions where:

For each suspension system we analyzed, the impeller speed which matched both criteria in Eq. 4 was selected as the recommended operating speed for the production scale. 3. Conclusions This study demonstrates a practical, experimentally anchored, and CFD-driven methodology for scaling up powder suspension processes. The use of measured properties at each powder loading ensures realistic representation of the suspension in the CFD model. Steady-state VOF CFD at pilot scale provides hydrodynamic descriptors—particularly bulk and surface velocities—that serve as robust scale-up targets. By matching these descriptors at large scale, the method achieves more reliable hydrodynamic similarity than traditional constant tip speed or power-per-volume approaches. The resulting workflow allows engineers to determine production-scale operating speeds without explicitly modelling powder addition, reducing development time, minimizing pilot trials, and enabling more confident scale-up. The methodology provides a generalizable digital framework for the design and optimization of industrial powder suspension mixing.

References [1]. Hirt, C.W., & Nichols, B.D. (1981). Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries. Journal of Computational Physics, 39, 201–225. [2]. Yakhot, V., & Orszag, S.A. (1986). Renormalization Group Analysis of Turbulence. Journal of Scientific Computing, 1, 1–51

ABSTRACT. A nuclear cooling system configured as a closed loop primarily generates forced convection using reactor coolant pumps. To simulate the internal closed-loop flow within a reactor pressure vessel driven by pumps, conventional computational fluid dynamics (CFD) analyses perform detailed modeling of the pump vanes and the surrounding flow field. This approach directly simulates the interaction between the rotor’s rotational motion and the surrounding fluid, providing the advantage of accurately predicting unsteady pressure fluctuations in the flow field. However, because extremely high computational costs are required to perform steady-state and transient analyses, it is difficult to apply this method to a wide range of operating conditions or to large-scale systems such as nuclear reactors. Thermal–hydraulic analysis techniques for nuclear cooling systems employ the Euler–Euler approach, which simultaneously treats the gas and liquid phases in order to account for the two-phase flow characteristics in the nuclear thermal hydraulics. In this framework, individual special models are developed and applied to the components constituting the reactor fluid system, such as the pressurizer, steam generator, and pump. Unlike the CFD-based pump flow model described above, the Euler–Euler approach treats the pump as a momentum source. Pump models based on momentum sources have the advantage of simplifying the analysis of multiple pump behaviors or large-scale reactor thermal–hydraulic systems. However, difficulties arise in simultaneously matching both pressure and flow rate to their design values, and numerically, the presence of large source terms can cause problems in reducing solver residuals or achieving stable convergence of fluid variables. In this study, as a solution to the convergence issues of the conventional momentum-source-based pump model within the Euler–Euler framework, a Direct Flow Control (DFC) technique that directly imposes mass into the flow field was developed. Unlike existing momentum-source-based flow control methods, the developed flow control approach directly prescribes the mass flow rate, enabling the target flow rate to be established rapidly and accurately. The proposed direct flow control technique directly applies the pump design flow rate at the faces of computational fluid cells; therefore, unlike momentum-source methods that apply a body force at the cell center, it does not require a new or modified mesh. In contrast to inlet velocity boundary conditions that directly impose flow at wall cell faces and require the generation of wall meshes for flow application, the proposed method does not require wall mesh generation and is capable of simulating reverse flow within the pump caused by adverse pressure conditions. Through this study, a method for directly imposing mass flow rates at computational cell faces was developed and applied to forced-convection thermal–hydraulic analyses of an entire reactor system driven by pumps. The currently developed model can be coupled with reactor pump performance curves in the future, enabling extension to simulations of various pump transients, and can be utilized to evaluate diverse three-dimensional local flow variations under reactor accident conditions.

ABSTRACT. Growing global concern over climate change, coupled with continued dependence on fossil fuels, has intensified the demand for sustainable and clean energy solutions. Among the various renewable energy sources, tidal energy is considered a promising alternative due to its high predictability, long-term sustainability, and relatively low environmental impact. Tidal turbines play a critical role in harnessing this resource, and enhancing their performance remains a significant research challenge. This study aims to enhance the computational efficiency of high-fidelity turbine performance prediction through the development of a surrogate modelling framework based on Blade Element Momentum Theory (BEMT). The proposed surrogate model is used to investigate the aerodynamic effects of a morphing trailing-edge design applied to tidal turbine blades, with particular emphasis on achieving an optimal balance between thrust and power output at the turbine level. The model is validated against BEMT calculations and demonstrates good agreement in predicting key turbine performance metrics. The results indicate that the morphing trailing-edge design yields a notable improvement in turbine performance, achieving an average power increase of 34.69% across the tip speed ratio in a range of 2 to 6. Furthermore, the validated surrogate model is employed to efficiently identify optimised design configurations, highlighting its potential for rapid performance evaluation and design optimisation of advanced tidal turbines.

ABSTRACT. Phase change and interfacial heat and mass transfer (solidification, melting, boiling, and related reactive or catalytic processes) are governed by bulk conservation laws coupled through jump conditions posed at interfaces that are extremely thin compared with the flow and mesh scales. Even in reduced settings, closed-form solutions are limited to simplified configurations, while realistic two- and three-dimensional geometries require numerical approaches that preserve interface sharpness, local and global conservation, and robustness in the presence of discontinuous material properties and strongly constrained geometries.

This work targets configurations in which scalar transport (temperature or species) interacts with incompressible flows and moving phase boundaries. Building on the sharp-interface two-fluid paradigm recently advocated for multiphase Navier–Stokes formulations, as an alternative to classical one-fluid treatments, we develop a conservative Cartesian cut-cell discretization that enforces boundary and interface conditions directly in the algebraic system, avoiding interface smearing and body-fitted remeshing. The approach is applicable to equilibrium closures (for example, Henry-type partitioning) and non-equilibrium kinetics, and naturally extends to Stefan-type interface motion. We consider two immiscible phases separated by a sharp interface. In each phase, incompressible flow is described by phase-specific velocity and pressure fields solving the two-fluid Navier–Stokes equations. A scalar field (temperature or species) is transported by the flow velocity and diffuses with phase-dependent properties. At the interface, conservation of the scalar yields a normal flux balance, while a closure relation prescribes the interfacial jump, for instance through a Henry-type partition law. For the two-fluid coupling, we enforce continuity of velocity and balance of interfacial traction. Spatial discretization relies on a Cartesian cut-cell finite-volume construction: control volumes are intersected by embedded boundaries and the interface, and operators are assembled from a minimal geometric description based on phase-restricted volumes, face apertures, and low-order moments. This yields compact stencils and symmetric definite operators in each phase, while supporting mixed boundary conditions and sharp interfacial coupling. In the two-fluid setting, we introduce bulk and interfacial degrees of freedom on each side of the interface, enabling enforcement of the jump relations in a coupled linear system. To accommodate moving geometries, we adopt a space-time viewpoint: each control volume is extruded over a time interval and the interface trajectory defines space-time cut-cells. The transient balance laws are then integrated on these space-time regions, providing a locally conservative update that accounts for the swept phase volumes and the fluxes through moving apertures. Time advancement uses an implicit scheme, and the space-time construction naturally supports repeated creation and removal of small phase volumes (fresh/dead cells) while preserving phase-wise conservation. For flow coupling, incompressible Stokes/Navier-Stokes equations are solved on the staggered embedded grid; the resulting phase velocities drive scalar advection and enter in interfacial balances. The approach is validated on two-phase heat and mass transfer configurations featuring strong interfacial partition effects, large diffusivity contrasts, and moving phase boundaries. The results demonstrate accurate resolution of interfacial fluxes and robust phase-wise conservation in two and three dimensions, including complex embedded geometries and evolving interfaces. These capabilities support predictive simulations of interfacial transport processes relevant to energy systems, materials processing, and multiphase thermal applications.

ABSTRACT. Junction flow is a flow formed by the interaction between a boundary layer and an obstacle, exhibiting pronounced three-dimensional characteristics. It generates horseshoe vortices in front of obstacles, which can have adverse effects in engineering applications, such as scour near structural joints or flow-induced noise. This study numerically investigates the development of junction flow under wave-influenced free-surface conditions. Specifically, this study analyzes the flow field generated by an infinite wave-current passing over a finite-length, finite-width flat plate with a cylinder mounted on it. The flat plate is submerged at a finite depth below the free surface, and both the waves and the incoming flow are parallel to the plate. Numerical simulations are performed using the commercial software FLUENT, employing the URANS approach with the SST k–ω turbulence model. The free surface is captured using the Volume of Fluid (VOF) method. The development of junction flow is examined for different water depths of the flat plate. The characteristic velocity is taken as the uniform inflow velocity, and the characteristic length is the cylinder diameter, yielding a Reynolds number of approximately 5.5×10^5. The results show that the horseshoe vortex in front of the cylinder varies periodically with the incident waves. When a wave crest approaches the cylinder, the primary horseshoe vortex moves closer to the cylinder and exhibits higher vorticity; conversely, when a wave trough approaches, the opposite occurs. This behavior is caused by the variation of water particle velocities within the waves. The submergence depth of the flat plate also influences the wave-induced effects. The closer the plate is to the free surface, the stronger the wave effect, leading to more pronounced changes in the horseshoe vortex position and strength. In addition, waves affect the downstream vorticity distribution of vortices shed from the horseshoe vortex. Different plate submergence depths not only alter the shedding period of the downstream vortices but also create a wake period that differs from the wave encounter period, resulting in an overall non-periodic vortex shedding pattern in the wake.

ABSTRACT. This study investigates the transient growth of perturbations within a laminar rotor-stator cavity following an impulsive acceleration of rotor speed, a critical transient process in aircraft engine air systems. As the limitations of traditional linear stability theory in characterizing short-term dynamics over time-dependent base flows, the research employs nonmodal stability analysis. A time-varying base flow is constructed based on a self-similarity assumption for a scenario where the Reynolds number instantaneously increases from 1000 to 2000. The analysis reveals that transient growth is a significant feature of the flow, with a maximum energy amplification of 5541.4 observed at an azimuthal wavenumber of 18. Physically, the optimal perturbations manifest as circular waves propagating radially within the stationary-disk boundary layer. These findings demonstrate that substantial short-time energy amplification can occur even in linearly stable, low-Reynolds-number regimes, providing a theoretical upper bound for the flow's potential transition to turbulence or complex states under linear dynamics.

ABSTRACT. In practical ramjet engines, combustion instabilities induce significant pressure oscillations, which affect the shock train structure in the isolator and introduce considerable operational risks. This study investigates the modulation of the shock train in an isolator using suction and injection. The results demonstrate that both flow intervention methods are effective but operate through distinct mechanisms. Suction creates a low-pressure zone, removes the near-wall boundary layer, and suppresses flow separation, thereby stabilizing the shock train downstream (7.62%L) and enhancing its stability. In contrast, microjet injection establishes a favorable pressure gradient zone, which shifts the shock train downstream with a smaller displacement (2.70%L). The larger separation bubble on one side within the shock train demonstrates significantly greater streamwise development, whereas the smaller one on the opposite wall expands predominantly in the spanwise direction. While suction simplifies design and improves stability, it incurs flow loss. Conversely, injection preserves the captured airflow but brings flow unsteadiness. Suction and injection both can reduce flow distortion and improve flow quality at the exit when the shock train is effectively controlled. Suction effectively dampens self-excited oscillations by removing part of the leading-edge separation bubble, whereas injection introduces additional unsteadiness. They both concentrate the oscillation intensity in the shock train leading edge. The oscillation frequency of the separation bubble (30-40 Hz) is markedly lower than both the primary (360-393 Hz) and secondary (131-139 Hz) characteristic frequencies of the leading-edge shock wave. Resolving this multi-scale characteristic is fundamental to constructing a comprehensive physical model of the shock train phenomenon.

ABSTRACT. Computational Fluid Dynamics (CFD) is an essential tool for simulating flow and heat transfer in aerospace, energy, and power applications. However, its high fidelity simulations are often hindered by excessive computational time and resource demands, limiting wider engineering adoption. While traditional data driven surrogate models and machine learning methods offer some relief, they still face challenges such as high training costs due to large scale data requirements, strong data dependency, and limited simulation accuracy. To overcome these limitations, this study proposes a physics‑aware domain decomposition reduced order model (DD ROM). First, a K‑means weighted sampling method partitions the physical domain into spatial regions. Data points are selectively sampled by assigning weights according to each region’s contribution to the feature distribution, thereby reducing dataset size while retaining core physical features. Next, feature decomposition is applied to reduce the dimensionality of the sampled data, projecting high‑dimensional data into a low‑dimensional subspace. A deep learning reduced order model, built on a fully connected neural network, is then trained to achieve fast reconstruction of flow and heat transfer states through nonlinear fitting and modal superposition. Simulations under both steady and unsteady conditions demonstrate that the proposed DD ROM significantly reduces storage requirements. Comparative experiments on nozzle flow and heat transfer prediction—benchmarked against a Dynamic Mode Decomposition (DMD) ROM—show that the DD ROM achieves higher predictive accuracy and better generalization. Moreover, its memory footprint during model generation and prediction drops sharply to 0.07 MB and 0.01 MB, respectively, far below the DMD ROM’s 7.94 MB and 3.29 MB. This model effectively alleviates the computational burden of processing high‑dimensional large‑scale data, offering an efficient approach for rapid prediction of flow and heat transfer in complex engineering systems.

ABSTRACT. This study presents a comprehensive numerical investigation of the coupled thermo-fluid and adsorption-desorption dynamics within a desiccant air dryer. To overcome the limitations of conventional experimental approaches and one-dimensional modeling, a transient axisymmetric model was developed in ANSYS Fluent, incorporating user-defined functions (UDFs) to describe vapor mass transfer based on the linear driving force (LDF) approach and the equilibrium adsorption capacity using the Sips isotherm. The numerical framework was validated against experimentally measured breakthrough curves, demonstrating good agreement. The results elucidate the internal transport phenomena by resolving the transient evolution of velocity, molar concentration, and temperature fields. During the adsorption process, the propagation of the adsorption front is characterized by localized temperature peaks resulting from the exothermic heat of adsorption. In contrast, the desorption process, driven by high-temperature purge gas, exhibits pronounced endothermic cooling in regions of active desorption. The analysis further reveals that the discrepancy between the equilibrium capacity and the instantaneous loading served as the primary driving force for mass transfer, guiding the system toward thermodynamic equilibrium. Overall, this research provides a high-fidelity numerical tool for probing internal thermal and flow characteristics that are difficult to access experimentally, offering a robust framework for the design and optimization of industrial air-drying systems.

Acknowledgement This study was conducted with the support of the Korea Institute of Energy Technology Evaluation and Planning (KETEP) with funding from the government (Ministry of Trade, Industry and Energy) in 2026.

ABSTRACT. Accurate prediction of conjugate heat transfer across interfaces among porous media, fluids, and solids under local thermal non-equilibrium (LTNE) is critical for the design of advanced thermal systems. Conventional numerical approaches typically impose empirical interfacial conditions that compromise thermodynamic consistency. We propose a lattice Boltzmann model that inherently satisfies thermal continuity at material boundaries without ad hoc matching. The formulation introduces a generalized volumetric heat capacity ratio to unify the macroscopic energy equation across heterogeneous domains. By appropriately designing the source term, the modified lattice Boltzmann equation recovers the macroscopic energy equation with second-order accuracy through Chapman–Enskog expansion. Validation against benchmark configurations—including forced convection in a circular pipe partially filled with porous media—demonstrates quantitative agreement with established interface-resolved simulations. The model captures the dependence of heat transfer enhancement on porous layer thickness and resolves the transition between LTNE and local thermal equilibrium regimes. This mesoscopic framework enables high-fidelity, multiphysics simulation of conjugate heat transfer in complex heterogeneous media for next-generation thermal engineering applications.

ABSTRACT. This study develops a high-order parallel solver for unstructured grids based on the Flux Reconstruction/Correction Procedure via Reconstruction (FR/CPR) method, implemented in Fortran. By integrating four limiters—TVB, Modal, KXRCF, and CKXRCF—numerical simulations were conducted for two classic test cases: double Mach reflection and shock–vortex interaction. The results demonstrate that all limiters effectively capture discontinuous flow features, clearly resolve shock morphology and their multiple reflections with walls, and exhibit well‑defined hierarchical wave structures. This verifies the reliability and efficiency of the method for high‑accuracy computational fluid dynamics simulations.

ABSTRACT. Thermo-compression bonding (TCB) using non-conductive film (NCF) is a critical technology for achieving fine-pitch interconnections in high bandwidth memory (HBM). However, the complex multiphase behavior involving molten solder and NCF, coupled with the difficulty of real-time experimental observation, poses significant challenges for process optimization. This study proposes a unified numerical framework designed to capture the fully coupled dynamics of the TCB process. Integrating an enthalpy-based multiphase flow formulation with dynamic meshing, the model accounts for governing physical phenomena, including solder phase change, wetting-driven shape evolution, and the curing kinetics-dependent rheology of the NCF during die displacement. This physics-based approach provides a systematic tool for investigating the interplay between material properties and process parameters, offering practical guidance for the design and optimization of high-reliability TCB processes.

ABSTRACT. We present a high-performance, two-way coupled CFD--LPT framework for particle-laden flows spanning dilute suspensions and locally dense regimes. To enable large-scale simulations, the solver employs MPI-based domain decomposition and an FFT-based Poisson solver that accelerates the pressure-projection step. Interphase coupling is implemented with conservative Gaussian particle-to-grid (P2G) projection and second-order grid-to-particle (G2P) interpolation on a staggered grid. A robust PP–Ergun hybrid closure models fluid–particle interactions, using adaptive weighting to ensure smooth transitions from dilute drag to dense packed-bed resistance. Verification is conducted for (i) single-sphere sedimentation, where the present PP prediction closely overlaps recent particle-resolved results, achieving near PR-DNS accuracy at a fraction of the computational cost, and (ii) flow through periodic FCC packed beds over $\phi = 0.15$--0.55 and $Re_{sup} = 50$--1000. In the packed-bed tests, the hybrid model reproduces the expected monotonic increase of $F_d^{*}$ with $\phi$ and follows DNS-based FCC trends; compared with a Gidaspow-type hard switching closure, it yields a smoother response across the dilute--dense transition and is consistently closer to particle-resolved FCC correlations in the moderately dense range. The results confirm rigorous momentum conservation and numerical robustness, establishing a scalable foundation for future four-way coupled simulations of wall-bounded turbulent deposition.

ABSTRACT. High-speed underwater separation is a critical process for cross-media vehicles, yet the underlying unsteady hydrodynamics and multi-body interactions remain unresolved. This study numerically investigates the axial underwater separation process of a vehicle–shell system to clarify the physical mechanisms governing its hydrodynamic characteristics. A multiphase flow model utilizing the Reynolds-averaged Navier–Stokes (RANS) equations, the Volume of Fluid (VOF) method, the Schnerr–Sauer cavitation model, and an overset grid technique is established. Based on the kinematic and flow features, the separation process is systematically divided into a loading phase, an unloading phase, and a cruise phase. The effects of separation force, initial velocity, and initial water depth on the system's separation behavior are comprehensively analyzed. Results indicate that separation force is the dominant factor; a higher force significantly shortens the separation time, while an insufficient force leads to separation failure. Furthermore, a higher initial velocity prolongs the total separation time. Conversely, a lower initial velocity can induce a negative velocity state that triggers severe cavitation inside the shell. Moreover, deep-water environments generate a high ambient backpressure, which creates a sharp pressure difference that exerts a backward drag force on the shell, thereby significantly delaying the overall separation process. These findings provide valuable theoretical guidance for the design and optimization of high-speed underwater separation systems for cross-media vehicles.

ABSTRACT. The Shattered Pellet Injector (SPI) system is a crucial technology for mitigating plasma disruption in large tokamak-type fusion reactors, specifically the International Thermonuclear Experimental Reactor (ITER). The HUN-REN Centre for Energy Research (CER) Fusion Plasma Physics Laboratory built and commissioned a test laboratory to help the development of ITER’s SPI for its Disruption Mitigation System. The 27 injectors placed around the torus will form 28.5 mm diameter cryogenic solid hydrogen pellets and launch them towards a shattering head at several hundred meters per second. To effectively mitigate disruption, it is crucial to minimise the amount of propellant gas that reaches the plasma before the solid pellet fragments. The SPI system employs a gas retention or suppressor system, which has undergone several recent design updates. The suppressor design is assisted by Computational Fluid Dynamics (CFD) simulations; however, accurately modelling the problem is challenging with traditional numerical methods. The volume is initially evacuated to 10e−7 mbar, while the propellant gas pressure is around 80 bar. A significant portion of the domain operates in the rarefied regime, with Knudsen numbers well above unity, making classical finite volume schemes unreliable in these areas, particularly for modelling the expansion front. Mean while, the high-pressure regions and strong temporal dependence make kinetic-theory-based methods computationally inefficient.

We are developing a CFD solver that is based on the analytical solution of the Riemann-vacuum problem. The program solves the discretised Euler equations, our code introduces a new flux definition at the gas-vacuum boundary. This modification allows the exclusion of vacuum cells from the calculations, improving stability and applicability. By extending the model into higher dimensions, we apply the method to model the propellant gas expansion within the SPI system. This process necessitates a technique to geometrically track the evolution of the gas-vacuum boundary in space, achieved through a piecewise-linear tracking method. We define a set of boundary points for each partially filled cell and update their positions over time according to the constitutive equations. This enables interface tracking below the cell resolution. We discuss the applicability of multiple flux schemes in the gas-gas regions and their effect on the solution. The consistency of the two-dimensional solver is assessed by comparing its output with 1D models. The results show excellent agreement with the analytical solution. We compare our data with simulations made with the commercially available ANSYS Fluent software, which is used for the SPI suppressor design. The suitability of our approach for modelling the SPI system and other related problems is also discussed, as well as limitations of the solution and possible enhancements for further applications.

ABSTRACT. 1. Introduction Water-entry problems have been extensively investigated for their critical engineering value in a variety of applications, including lifeboat launching [1] air-dropped torpedoes [2], and AUV release [3]. During the water entry process, limited launch heights or near-surface target locations often result in shallow-angle entry. This results in complex hydrodynamic characteristics and flow field evolution, significantly impacting structural integrity and motion behavior. Therefore, understanding the flow evolution and motion characteristics of vehicles with different geometries during water entry holds substantial guiding importance for practical engineering applications. In the early stages, scholars primarily investigated the water entry problem through experimental and theoretical approaches. Worthington [4] was the first to capture images of splashes and cavities produced by spheres under different conditions using flash photography techniques. Subsequently, May [5] conducted extensive water-entry experiments on objects of various shapes. On the theoretical side, Logvinovich [6] proposed the independent expansion principle to describe the formation and evolution of cavities during water entry. Later, Lundstrom et al. [7] summarized the development patterns of cavities generated by high-speed vehicle water entry and established models for the resulting impact loads. Collectively, these contributions have established a comprehensive supercavitation dynamic framework for water-entry problems. With the rapid development of computational technology and computational fluid dynamics, numerical simulations have been increasingly used to study water-entry problems. In terms of flow-field characteristics, Yi et al. [8] examined how cavity shapes change when small spheres enter rotating water. Kintea et al. [9] performed numerical studies on the water entry of rotating and non-rotating spheres. Kiara et al. [10] studied the effects of contact location, hydrodynamic forces, and rotational motion on cylinder water entry. Sun et al. [11] used fluid-structure interaction (FSI) methods to investigate cavity evolution and load characteristics during vertical cylinder water entry. In summary, current research on water-entry problems has mainly focused on flow-field characteristics of simple geometries, such as spheres and cylinders, as well as the effects of structural parameters on the stability of vehicles entering water at large angles. However, in practical engineering applications, underwater vehicles often have multi-component geometries. Therefore, studying the water entry of vehicles with complex structures has significant practical value. Moreover, systematic studies on shallow-angle water entry remain limited, and the definition of "shallow angle" is often vague. In this work, the relationship between vehicle structure and entry angle is used as a basis to systematically investigate the shallow-angle water entry of a truncated-cone vehicle. Three entry models are defined, and the flow-field characteristics and motion behavior under each scenario are compared. The results provide valuable guidance for engineering applications. 2. Methodology 2.1. Governing equations We consider two-phase incompressible flow in this study due to the vehicle entering the water at the speed of 100 m/s. A three-dimensional simulation model for the high-speed water entry of vehicles is developed using ANSYS Fluent software, combining the finite volume method, Volume of Fluid (VOF) multiphase flow model, and overset grid technique. The compressibility of air and water is considered. The governing equations are expressed as follows: (1) (2) (3) (4) Where is the density of the mixed phase, is the velocity ; is pressure, is gravitational acceleration, is the eddy viscosity, and are the volume fraction of the vapor and air, , and are the density of the water, vapor and air, respectively, is the kinetic viscosity of the mixed phase, and are kinetic viscosity of the vapor and air, The subscript is the direction in the Cartesian coordinate system. Energy equation: (5) Where is the internal energy, , is the specific heat capacity, is the temperature ; is the temperature, is heat transfer coefficient of fluid, is the increase in internal energy caused by viscous forces. 2.2. Turbulence model The Realizable turbulence model is employed to solve complex flows, offering better predictive accuracy for high-Reynolds-number flows. The transport equations for turbulent kinetic energy and its dissipation rate are expressed as follows: (6) Where and are the turbulent kinetic energy generated by mean velocity gradients and buoyancy, is expansion to the total dissipation rate in compressible turbulence. is the turbulent viscosity coefficient, and are the turbulent Prandtl numbers for and , , , is the mean strain rate tensor, , is the kinematic viscosity, , and are the model coefficients, , , , . 2.3. Cavitation model During the water entry process, the vehicle experiences cavitation phenomena, which occur when the water pressure drops below the saturation vapor pressure, causing liquid water to vaporize. This word employs the Schnerr-Sauer cavitation model to solve the water-vapor mass transfer process: (7) (8) Where and are the empirical calibration coefficients of evaporation and condensation, is the cavity radius, is the saturated vapor pressure, and are the mass transfer rate of evaporation and condensation, , , . 3. Conclusions Due to its specific geometric configuration, characterized by the truncated-cone angle c, the truncated-cone vehicle can be classified into three water-entry modes at shallow angles based on the relationship between the truncated-cone angle c and the entry angle (Figure 1). M-a( < c): initial contact occurs at the shoulder. M-b( = c): simultaneous contact at the shoulder and cavitator. M-c( > c): Initial contact occurs at the cavitator. Owing to the differences in the initial contact locations, the vehicle motion and cavity evolution vary among the three entry modes. Therefore, numerical simulations are performed for entry angles = 2 , 3 , 4 , 5 , and 6 at an entry velocity of V0 = 100 m/s. The differences among the entry modes are analyzed mainly in terms of cavity evolution and vehicle motion characteristics. The main conclusions are summarized as follows: (1) Based on the surface wetting of the vehicle and the evolution of the flow field during water entry (Figure 2), the main differences among the three entry modes are analyzed. During the cavity formation stage, M-a exhibits the most extensive wetted area and the longest duration to achieve symmetry between the upper and lower cavity walls. In contrast, M-b and M-c demonstrate high similarity in flow-field evolution, and the difference in the time required for cavity symmetry is relatively small. During the underwater motion stage, the differences are mainly reflected in the surface wetting during tail slapping. M-a induces bilateral wetting at both the shoulder and tail, whereas M-b and M-c experience wetting exclusively at the tail. In addition, due to compression by water vapor, the air pockets at the tail in M-b and M-c detach earlier. (2) From the motion characteristics during water entry, M-a induces the highest hydrodynamic loads to its extensive wetted area. As a result, the vehicle experiences rapid velocity decay and a relatively large trajectory deviation. In M-b, the contact area at the initial impact is the largest, which generates a larger instantaneous pitching moment. In addition, with the increase in entry angle, the wetted area on the vehicle surface decreases rapidly, leading to reduced velocity decay and trajectory deviation. For M-c, the wetted area varies only slightly with increasing entry angle. Consequently, the velocity decay and trajectory deviation also show relatively small differences.

ABSTRACT. It is well known that, for initial–boundary value problems, boundary closures of high-order finite-difference schemes are often deliberately reduced in order to ensure stability. However, such order reduction can degrade the observed global convergence rate of high-order methods. To enhance the stability of the fifth-order linear Hermite–WCNS (H-WCNS) scheme, we derive high-order accurate boundary closures that satisfy a global conservation property by introducing nonuniform near-boundary solution points together with associated quadrature weights. The resulting globally conservative H-WCNS scheme preserves the designed high-order accuracy on the modified semi-uniform grid and enforces a discrete global conservation relation. Eigenvalue analysis of the semi-discrete operator indicates that the proposed scheme remains strictly stable. Additional numerical tests further demonstrate that the globally conservative formulation delivers improved overall performance.

ABSTRACT. The Ahmed body is a widely used benchmark in vehicle aerodynamics, as it reproduces the main flow features of road vehicles while maintaining a relatively simple geometry. The flow topology around the rear slanted surface is strongly dependent on the slant angle, leading to markedly different separation patterns, wake structures, and aerodynamic loads. In particular, rear slant angles of 25° and 35° are known to generate complex three-dimensional separation and recirculation phenomena, which pose significant challenges for numerical simulations.

In this work, Large Eddy Simulation (LES) combined with an Immersed Boundary Method (IBM) is employed to investigate the flow around the Ahmed body for rear slant angles of 25° and 35°. The simulations are performed using the open-source CFD solver CHAMÁN, which solves the incompressible Navier–Stokes equations using a finite-volume approach on Cartesian grids. Three subgrid-scale turbulence models are evaluated: the classical Smagorinsky model, the Wall-Adapting Local Eddy-viscosity (WALE) model, and the Anisotropic Minimum Dissipation (AMD) model. Two different mesh resolutions are considered, with minimum cell sizes of 2 mm and 1.5 mm in the vicinity of the body.

The numerical results are assessed through comparisons with available experimental and numerical reference data from the literature. Mean velocity profiles along the rear slanted surface show good agreement with experimental measurements, particularly for the 35° configuration, which exhibits low sensitivity to mesh resolution and turbulence modeling. The 25° configuration is more demanding, with improved predictions obtained using the finer mesh, especially when combined with the WALE model. Aerodynamic coefficients derived from pressure integration decrease with mesh refinement and are generally overestimated relative to reference data, although the best agreement is achieved using the finer mesh and the Smagorinsky model. The results demonstrate the capability of the IBM–LES framework to reproduce the main aerodynamic features of the Ahmed body while highlighting the influence of near-wall resolution and turbulence modeling.

ABSTRACT. Cavitation induces severe operational challenges in engineering systems, ranging from material erosion to performance degradation. As the critical precursor to developed cavitating flow, the regulation of incipient cavitation represents a research frontier with profound engineering implications and scientific merit. In this work, a numerical investigation is conducted to systematically explore the modulation mechanism of vortex generators (VGs) on incipient cavitation dynamics. The numerical framework incorporates a modified cavitation model coupled with an advanced turbulence closure, enabling quantitative analysis of incipient cavitation characteristics under varying VG installation angles. The results show that VGs with an appropriate installation angle effectively suppress near-wall flow separation. Excessively large angles increase the trailing-edge spacing between paired VGs, weakening vortex mixing in the intermediate region and thereby compromising separation inhibition. As the cavitation number decreases, the initial cavitation manifests as vortex-induced cavitation generated by the vortex generator, with no attached primary cavitation occurring on the wall. Comparative analysis demonstrates a positive correlation between vortex cavitation scale and the intensity of VG-induced counter-rotating vortices. In single-phase flow, vortex generators (VGs) significantly reduce hydrofoil drag coefficient. At small installation angles, they further increase the lift coefficient. When cavitation occurs, hydrofoil drag increases slightly, while lift variation remains consistent with single-phase flow behavior.

ABSTRACT. Yield-stress fluids are common in engineering and geophysical applications such as cement grouting and fracture-scale flows. Classical Bingham and Herschel–Bulkley models are widely used but neglect thixotropy, which introduces time-dependent structural breakdown and rebuilding. This study examines the steady radial diverging flow of a thixotropic grout between parallel plates using a structural-kinetics model. Despite steady macroscopic flow, spatial variations in the internal structure arise due to advection and shear-rate gradients. Comparisons with an equivalent Bingham solution reveal conditions under which steady thixotropic flows depart from time-independent yield-stress behavior.

ABSTRACT. The Molten Salt Fast Reactor (MSFR), categorized as a Generation IV system with circulating liquid fuel, necessitates high-fidelity multiphysics analysis because of the pronounced interaction between neutronic and thermal-hydraulic phenomena. An advanced steady-state multiphysics coupling framework is developed for MSFR analysis, combining LES-based thermal-hydraulics implemented in the CUPID code with GPU-accelerated Monte Carlo neutronics in PRAGMA through the preCICE coupling interface. Through iterative data exchange, the coupled solvers consistently resolve power distributions and temperature–density feedback until convergence of the steady-state solution. The steady-state results for the MSFR reference configuration show that the proposed framework accurately reproduces essential reactor characteristics, including power–temperature feedback and flow recirculation, thereby validating its applicability to high-resolution MSFR multiphysics studies.

ABSTRACT. This work presents a discrete loss-based optimization solver (FDTO) that integrates conservative finite-difference discretization with a time-stepping optimization strategy for solving nonlinear partial differential equations (PDEs) on several geometries. By formulating the governing equations as a residual-based loss and performing node-wise optimization on body-fitted structured meshes, the proposed FDTO framework preserves the physical consistency and numerical stability of classical computational fluid dynamics (CFD) schemes while maintaining the flexibility of optimization-based solvers. Benchmark simulations of incompressible Navier–Stokes flows demonstrate that the method achieves high accuracy and robust convergence across a wide range of Reynolds numbers. Compared with existing ODIL-based approaches, the proposed framework exhibits improved generalization capabilities and faster error decay, particularly in strongly nonlinear flow regimes. Moreover, the method successfully captures key flow structures, including multi-vortex dynamics at high Reynolds numbers, confirming its predictive capability. Overall, the proposed solver offers a geometry-aware, matrix-free alternative to traditional CFD and surrogate modeling methods, providing potential for efficient simulation of unsteady flows and complex engineering configurations.

ABSTRACT. This study presents a comprehensive transient CFD-based performance evaluation of a large-capacity twin-screw blower specifically designed for pneumatic ash handling systems in coal-fired power plants. To capture the complex, time-varying pressure and mass-flow characteristics dictated by rotor meshing and discharge-port timing, Unsteady RANS (URANS) simulations utilizing the $k–\omega$ SST turbulence model were conducted within the ANSYS CFX framework, facilitated by high-quality rotor-domain grids generated via TwinMesh.Numerical investigations were performed across three discharge-pressure conditions ($0.5, 0.8, \text{ and } 1.1\text{ kgf/cm}^2$) at synchronized rotor speeds (male: 1750 rpm; female: 875 rpm). The results quantify a distinct performance sensitivity: as discharge pressure escalated, the mean outlet mass flow rate exhibited a downward trend, whereas torque, power requirements, and outlet temperatures showed substantial increases, indicating intensified compression work and elevated thermal loading under off-design conditions.Notably, intermittent backflow was identified toward the end of the discharge phase, occurring when the discharge-side static pressure exceeded the internal rotor-chamber pressure—a phenomenon that intensified with increasing backpressure. Flow-field inspections further highlighted a biased discharge pattern characterized by housing wall impingement and local recirculation within the discharge region. These features are significant contributors to additional energy losses and the amplification of unsteady pressure fluctuations. This study provides critical quantitative insights for optimizing operating strategies and refining discharge-port and duct geometries to enhance the efficiency and operational reliability of large-capacity screw blowers in industrial ash handling applications.

ABSTRACT. In a gas-generator-cycle liquid-propellant rocket engine, the hot combustion gases produced in the gas generator are discharged through the turbine exhaust system (TES). Despite its seemingly simple role of expelling hot gases, the TES contributes to engine performance in several ways. However, the fundamental fluid-dynamic behavior of the TES remain relatively less understood, particularly from the perspective of fluid-structure interaction. Motivated that combustion-produced soot can accumulate on internal wall surfaces and thus alter the effective geometry, we numerically investigate the effects of soot deposition on the fluid-dynamic characteristics of the TES. As a fundamental investigation into the fluid-structure interaction of the TES, geometric changes due to the soot accumulation are represented by varying the dimensions of representative internal structural features, and their effects on flow behavior and pressure drop are analyzed comprehensively.

ABSTRACT. 1. Introduction Liquid film cooling via oil-jet impingement is an important thermal management strategy for high-power electric vehicle (EV) motors. The cooling efficiency is fundamentally governed by the spreading dynamics of the thin liquid film on the coil surfaces. However, despite its industrial importance, comprehensive numerical studies that simultaneously resolve the multiphase hydrodynamics and conjugate heat transfer (CHT) specifically for concentrated winding configurations are notably scarce. Most existing works simplify the wetting physics by relying on a constant Static Contact Angle (SCA) assumption [1, 2]. This simplification ignores the velocity-dependent nature of the moving contact line, often leading to inaccuracies in predicting the effective wetted area. To overcome this limitation, precise modeling of the contact line motion is essential. This study aims to numerically investigate the hydrodynamic behavior and cooling performance of the oil film on a simplified coil surface, chosen to isolate the fundamental wetting physics while addressing the lack of coupled multiphase-thermal simulations in this domain.

2 Methodology 2.1 Numerical Framework The simulations are conducted using a custom solver developed within the OpenFOAM environment. The flow field is resolved by solving the transient incompressible Navier-Stokes equations. To capture the sharp gas-liquid interface, the Volume of Fluid (VOF) method is utilized, specifically employing the geometric VOF (isoAdvector) algorithm to minimize numerical diffusion and interface smearing. The phase distribution (α) is governed by the transport equation [3]: ∂α/∂t+∇∙(uα)=0 , where u denotes the velocity field. Furthermore, the temperature-dependent viscosity of the cooling oil is modeled to account for the thermal thinning effect. To accurately model the wetting dynamics at the wall, a Dynamic Contact Angle (DCA) boundary condition is incorporated into the solver framework, which adjusts the local contact angle θ_d as a function of the contact line velocity U_cl. Simultaneously, the energy equation is tightly coupled with the flow solver to evaluate the thermal characteristics of the oil jet. To account for the thermal interaction between the fluid and the solid coil, a conjugate heat transfer (CHT) approach is implemented, strictly enforcing the continuity of temperature and heat flux at the fluid-solid interface.

2.2 Parametric Study Comparative simulations are designed to quantify the deviation in wetting dynamics between the SCA and DCA models. Additionally, the inlet flow rate (Reynolds number) is varied to analyze how jet momentum influences the liquid film footprint. The cooling performance is evaluated using the local and average Heat Transfer Coefficient (HTC).

3. Conclusions In this study, the oil jet cooling process is analyzed to evaluate the impact of wetting dynamics. Current simulation results using the SCA model have established a comprehensive baseline for the wetted surface area and thermal distribution (Fig. 1). However, as the SCA assumption neglects the contact angle hysteresis inherent in rapid film expansion, ongoing work focuses on deriving the corresponding DCA results. The final comparative analysis will quantify the deviation introduced by dynamic wetting physics, highlighting its critical role in accurately predicting the wetted surface area for high-performance EV motor cooling applications.

References [1] B. W. Webb and C.-F. Ma, "Single-phase liquid jet impingement heat transfer," Advances in Heat Transfer, vol. 26, pp. 105–217, 1995. [2] M. M. Rahman, A. J. Bula, and J. E. Leland, "Conjugate heat transfer during free jet impingement of a high Prandtl number fluid," Numerical Heat Transfer, Part B: Fundamentals, vol. 36, no. 2, pp. 139-162, 1999. [3] H. Roenby, H. Bredmose, and H. Jasak, “A computational method for sharp interface advection,” Royal Society Open Science, vol. 3, no. 11, p. 160405, 2016.

ABSTRACT. This study investigates the aerodynamic effects of incorporating multiple chordwise slits into a single-slotted aircraft flap. At the high deflection angles used during takeoff and landing, boundary-layer separation develops on the flap upper surface, and its extent varies along the span owing to three-dimensional flow. To address this issue, chordwise slits are introduced to segment the flap, with the aim of suppressing separation and enabling spanwise variations in flap geometry. As a first step toward spanwise-adaptive flap designs, we assess the aerodynamic impact of the slits and the sensitivity to geometric parameters. Numerical simulations are conducted for several configurations with different slit lengths, numbers, and widths, and the resulting lift characteristics are compared with those of a conventional flap. The results show that an optimal configuration increases the lift coefficient by up to 4.7% relative to the non-slit baseline. This improvement, however, involves a tradeoff: the slits promote earlier upstream pressure recovery, increasing the pressure near the flap leading edge.

ABSTRACT. Uncertainty quantification (UQ) in computational fluid dynamics for aerodynamics, particularly involving imprecise boundary conditions and defects in airfoil geometry, incurs a significant computational cost. Addressing this challenge demands a scalable and efficiently parallelized method for approximating solutions to the Navier--Stokes equations, for which the lattice Boltzmann method (LBM) is well suited. High-resolution UQ for two-dimensional incompressible flow past airfoils at $Re \geq 10^4$ remains unexplored since direct numerical simulations (DNS) are expensive. Furthermore, traditional numerical methods, such as finite-volume methods (FVM), require mesh deformation algorithms to account for geometrical defects, which introduce additional implementation overhead and a source of simulation failure. A solution to these problems is to use LBM in the efficiently parallelized \textbf{OpenLB-UQ} C++ library~\cite{openlb-uq}.

In this work, a framework for simulating incompressible flow past arbitrarily shaped airfoils using LBM is developed and validated on a uniform Cartesian lattice for Reynolds numbers up to $Re = 2 \times 10^4$ and across different angles of attack. Furthermore, a detailed approach for computing the skin-friction coefficient on a structured lattice, along with a grid-convergence study of all aerodynamic coefficients, is presented for the first time. The framework is subsequently used to carry out UQ experiments in a four-dimensional uncertainty space using the quasi-Monte Carlo (QMC) method. Statistical convergence results for the mean, standard deviation, and the $1$-Wasserstein distance indicate both expected convergence orders and an accurate approximation of the \textit{pushforward} measure of the force coefficients for Reynolds numbers up to $Re = 10^4$ at two angles of attack $\alpha = 0^\circ$ and $\alpha = 5^\circ$.

\section{Methodology} The standard LBM-BGK (Bhatnagar--Gross--Krook) formulation is employed to ensure efficient parallelization on multi-CPU/GPU configurations. Moreover, the class/shape function transformation (CST) is utilized for airfoil shape parameterization~\cite{cst}, since CST parameters directly influence airfoil properties, making it a versatile approach for optimization and UQ. The drag coefficient $C_D$ and the lift coefficient $C_L$ are computed using the momentum exchange algorithm (MEA) directly from the particle distribution functions. The pressure coefficient $C_p$ is computed simply by time-averaging the wall pressure, which is measured one lattice node away from the wall. Lastly, a combination of analytical normals obtained from the CST parameterization, a higher-order interpolated bounce-back boundary condition, and linear extrapolation/fitting is required to accurately compute the skin-friction coefficient $C_f$.

The non-intrusive UQ framework detailed in~\cite{openlb-uq} is used to conduct the UQ experiment at $Re = 10^4$ and $\alpha =0^\circ, 5^\circ$. The uncertainties comprise the boundary conditions: $Re, \alpha$; and the geometry: upper surface $S_u$, lower surface $S_l$; parameterized as \begin{equation}\label{eqn:parameterized-inputs-experiment} \begin{gathered} Re(\boldsymbol{Z}) = Re(1 + Z_1), \quad \alpha(\boldsymbol{Z}) = \alpha(1 + Z_2), \\[0.5em] S_u(\boldsymbol{Z}; A_i) = S_u(A_i)(1 + Z_3), \quad S_l(\boldsymbol{Z}; B_i) = S_l(B_i)(1 + Z_4), \end{gathered} \end{equation} respectively. For $Re = 10^4$ and $\alpha = 5^\circ$, the QMC algorithm is used to sample $M = 1500$ samples from the input uncertainty space $\boldsymbol{Z} = [Z_1, Z_2, Z_3, Z_4]^\text{T}$, where $Z_1\sim\mathcal{U}[-0.1, 0.1]$ and $Z_j\sim\mathcal{U}[-0.4, 0.4]$ for $j \geq 2$. Next, the OpenLB solver is used to propagate these uncertainties to determine the discrete probability measures of the quantities of interest $C_D$ and $C_L$.

\section{Results} \subsection{Validation and Convergence of Aerodynamic Coefficients} The validation of force coefficients (not presented here) in comparison to FVM and XFOIL is performed for flow past the NACA0012 airfoil up to $Re = 10^4$ for different angles of attack. Moreover, the surface coefficients are also validated under different flow conditions for flow past symmetric NACA airfoils. Fig.~\ref{fig:validation-cp-cf} represents the comparison between LBM and FVM~\cite{elimelechdns} for one such flow condition involving the flow past a NACA0009 airfoil at $Re = 2 \times 10^4$ and $\alpha = 3^\circ$. For both the upper surface and the lower surface, OpenLB simulated $C_p$ and $C_{fx}$ are in \textbf{strong agreement} with the reference. Only minor deviations are present at the upper surface leading edge for both $C_p$ and $C_{fx}$.

Convergence under grid refinement (\textit{diffusive scaling}) is performed for both the force and surface coefficients at $Re = 10^4$. A previous study has shown such convergence results for $C_D$ at $Re = 500$ in terms of accuracy. However, to the best of the author's knowledge, no studies have been published that demonstrate a rigorous convergence study of force coefficients for flow past airfoils at $Re = 10^4$. Furthermore, convergence studies for surface coefficients have not been reported at any Reynolds number. Fig.~\ref{fig:scaling-surface-coefficients} indicates that $C_p$ and $C_f$ converge at approximately \textbf{first-order} in the $L^1$-error norm (in terms of both accuracy and consistency). This result is significant for the present method of computing $C_f$, as it yields mesh-independent results with a healthy convergence rate under grid refinement.

\subsection{Uncertainty Quantification of Force Coefficients} Fig.~\ref{fig:pdf-uq} indicates the probability density functions of the force coefficients obtained in UQ experiments for $Re = 10^4, \alpha = 5^\circ$ with uncertainties according to Equation~\eqref{eqn:parameterized-inputs-experiment}. It can be observed qualitatively that there is a \textbf{substantial statistical spread} in the force coefficients. To assess the statistical convergence of the QMC-approximated \textit{pushforward} measure $\hat{\mu}_{\mathrm{QMC}}$ toward the reference measure $\hat{\mu}^{\Delta}$, the $1$-Wasserstein distance $W_1(\hat{\mu}^{\Delta}, \hat{\mu}_{\mathrm{QMC}})$ is plotted as a function of $M$ for the force coefficients in Fig.~\ref{fig:wasserstein-convergence}. Convergence of the $1$-Wasserstein distance directly implies convergence in the first and second-order moments. For QMC, it is expected that $W_1(\hat \mu^\Delta, \hat \mu_{\text{QMC}}) \sim \mathcal{O}(M^{-\lambda})$, where $\lambda \in [0.5,1]$. In Fig.~\ref{fig:wasserstein-convergence}, this is successfully achieved since $\lambda \approx 0.8$ for both force coefficients.

\section{Conclusions} \vspace{-0.5em} The three major novel contributions of this work are summarized below. \begin{enumerate}[topsep=0pt, partopsep=0pt, itemsep=0pt] \item \textbf{Detailed and accurate computation} of $C_f$ in LBM with a uniform Cartesian grid. \item \textbf{Convergence} of all aerodynamic coefficients under grid refinement. \item \textbf{Statistically convergent UQ} for flow past airfoils at $Re = 10^4$, accounting for imprecise boundary conditions and geometrical defects. \end{enumerate} These novel results pave the way for the following future research directions. \begin{enumerate}[topsep=0pt, partopsep=0pt, itemsep=0pt] \item \textbf{PDE foundation models}~\cite{poseidon-fm} have recently emerged as a promising approach for learning solution operators of PDEs. These models have demonstrated success in the fast and accurate inference of complex fluid flow phenomena. This motivates the use of foundation models in combination with LBM and QMC sampling to allow for UQ of airfoil flows at higher $Re$ and higher-dimensional input uncertainty spaces. \item Due to the turbulent nature of the flow past bluff bodies at high $Re$, increasingly smaller vortices are resolved as one approaches the inviscid limit $\nu \searrow 0$, i.e., with finer resolution and an increasing $Re$. The development of smaller vortices as the flow field is better resolved is a clear indication of \textbf{poor sample convergence} in the inviscid limit. The theory of \textbf{statistical solutions to the Euler equation}~\cite{statistical-solutions-iee} develops a concrete notion of convergence for these solutions in the sense of time-parameterized probability measures. The LBM-UQ framework developed in the present work shows capabilities for computing statistical solutions of incompressible flow past bluff bodies. \end{enumerate}

ABSTRACT. This study employs particle-resolved direct numerical simulation (PR-DNS) to investigate the complex physics of particle-laden impinging jet flows. Unlike conventional point-particle models, this approach resolves the full scale of turbulence near particles, providing reliable insights into surface sweeping and deposition patterns.The numerical framework utilizes the immersed boundary method (IBM) with direct forcing to solve the incompressible Navier-Stokes equations on a staggered Cartesian grid. This method incorporates a 4-point delta function for Lagrangian markers and has been enhanced to maintain stability for neutrally-buoyant and low-density-ratio particles. The computational domain maintains a fixed jet-to-plate distance of H/D = 2. A systematic parametric study was conducted by varying the particle Reynolds number, normalized particle size, density ratio, and volume fraction. Results indicate that increasing the volume fraction intensifies turbulence modulation, leading to nonlinear changes in particle transport that simpler models cannot capture. This research highlights the critical roles of particle inertia and size in shaping near-wall flow structures, offering a high-fidelity dataset for validating simplified industrial models.

ABSTRACT. In numerical simulations of fluid mechanics, achieving higher accuracy can be approached in several ways: refining the physical modeling of flows, increasing the resolution of discretization, or improving the precision of numerical methods. To balance computational cost and accuracy, anisotropic metric-based Adaptive Mesh has demonstrated its effectiveness for steady problems. However, extending this approach to unsteady flows remains a significant challenge due to the complexities introduced by time dependence. One key issue is the lagging of the mesh relative to the evolving solution. When a mesh is optimized for a solution at specified time, it quickly becomes suboptimal as the solution advances. Several methods in the literature address this issue, but they require introducing complex interpolation techniques to transfer the solution between meshes, as well as intricate metric definitions to ensure accuracy. In this work, we introduce a novel numerical method designed for unstructured meshes, ensuring conservation in both space and time, while mitigating these challenges in an efficient way.

Unlike classical time-marching approaches for unsteady problems, which typically rely on a single global time step imposed over the entire spatial domain, the present space-time setting naturally accommodates locally varying time increments.

ABSTRACT. High-fidelity CFD simulations provide accurate and predictive flow representations, but their computational cost limits parametric studies, control design, and real-time forecasting. Reduced-order models (ROMs) alleviate this by evolving the flow in a reduced representation identified from data, for instance via POD, DMD/Koopman-type operator-based approaches, or learned embeddings. Despite substantial progress, reliable long-horizon prediction and robust generalization across operating conditions and parameters (e.g., Reynolds-number variations) remain challenging, particularly when only sparse sensor measurements are available. We propose an uncertainty-aware, parameterized ROM tailored to sparse sensing that couples (i) a probabilistic encoder mapping sparse measurements and the parameter to latent coordinates, (ii) a parameter-conditioned latent dynamics model enabling fast multi-step rollouts, and (iii) a decoder reconstructing full-field quantities of interest from the latent state. Conditioning on the parameter is applied throughout the architecture via feature-wise linear modulation (FiLM) using per-feature scale and shift. The probabilistic formulation provides uncertainty-aware latent representations when sparse measurements do not fully determine the underlying state, enabling confidence estimates for temporal forecasting and parameter interpolation, and for assessing reliability under potential extrapolation.

ABSTRACT. Lithium-ion batteries (LIBs) in electric vehicles require effective thermal management, particularly under high C-rate operations in hot environments. While solid-liquid phase change materials (SL-PCMs) such as n-docosane, when integrated with cooling channels, offer effective thermal regulation, their requirement for complete resolidification between aggressive discharge-charging cycles poses a critical operational constraint. This study investigates neopentyl glycol (NPG), a solid-solid phase change material (SS-PCM), as an alternative to n-docosane for thermal management of a Samsung ICR 18650-26J lithium-ion battery. A three-dimensional numerical model based on the Multi-Scale Multi-Domain (MSMD) approach with the Newman-Tiedemann-Gu-Kim (NTGK) electrochemical model is employed to simulate the thermal behaviour during 3C discharging followed by 1C charging at 40 °C. The effects of NPG thickness (3–4 mm) and ambient temperature (30–40 °C) on cell temperature and phase fraction are systematically evaluated. A direct comparison between NPG (latent heat: 126 kJ/kg) and n-docosane (latent heat: 257 kJ/kg) reveals that while docosane achieves lower peak cell temperatures in a single cycle, NPG eliminates the resolidification requirement entirely by undergoing a reversible solid-solid phase transition. This characteristic makes NPG particularly advantageous for sustained multi-cycle operation in hot climates, where incomplete SL-PCM recovery progressively degrades thermal performance.

ABSTRACT. The pseudo high-speed (SHS) method effectively accelerates convergence in low-speed compressible flow simulations by elevating a pseudo Mach number while preserving physical fidelity. However, the original method exhibits poor density convergence due to its reliance on a linear density transformation, which fails to modify the density gradient. This deficiency propagates into the energy equation and undermines solution reliability—particularly in cross-velocity regimes such as transonic or supersonic flows where accurate prediction of aerodynamic heating is critical. To address this, we propose an improved SHS method based on a nonlinear density function transformation, which restores consistent and accurate convergence of both density and pressure fields without sacrificing computational efficiency.

ABSTRACT. 1. Introduction

Computational fluid dynamics (CFD) is a cornerstone of modern aerospace engineering, serving as an indispensable tool for design and analysis. However, the utilization of CFD is hindered by its steep learning curve. Recent advances in large language models (LLMs) offer a transformative opportunity to automate these complex CFD workflows through natural language interaction, potentially democratizing access to this powerful tool. However, their application to CFD is constrained by the need for precise physical understanding and software-specific syntax. This domain knowledge gap hinders effective automation of the complex CFD workflow, which demands a deep understanding of numerical schemes, turbulence models, boundary conditions, and solver-specific implementation details. Fine-tuning offers a direct approach for embedding domain expertise into LLMs, incorporating CFD knowledge into the model's parameters rather than relying on external retrieval [1].

We introduce CFD-copilot [2], an automated framework designed to manage the CFD workflow from case setup to post-processing through natural language. We developed NL2FOAM, a custom dataset comprising 49205 pairs of natural language descriptions and corresponding OpenFOAM configurations, augmented with chain-of-thought (CoT) annotations to capture expert reasoning. Using NL2FOAM, we fine-tuned Qwen3-8B model [3], enabling it to translate high-level natural language problem descriptions into executable CFD setups. A multi-agent system manages the workflow, handling input verification, configuration generation, simulation execution and error correction autonomously. We also integrated a model context protocol (MCP)-enabled post-processing agent, offering scalability and flexibility.The framework was evaluated on canonical aerodynamic benchmarks, including the NACA 0012 airfoil and the three-element 30P-30N high-lift configuration, to demonstrate its applicability to practical engineering problems.

2. Methodology

2.1. Fine-tuning LLM

LLMs, pre-trained on vast corpora of text, encode broad knowledge and language capabilities. They can be specialized for domain-specific applications through fine-tuning, particularly when sufficient labeled training data exists and high-precision is essential. For the CFD-specific adaptation, we fine-tuned Qwen3-8B [3] via Low-rank adaptation on the NL2FOAM dataset. Starting with 36 OpenFOAM cases spanning laminar and turbulent flows, we modified configuration files to generate over 200k variations. An LLM enhanced linguistic diversity by rephrasing the problem descriptions. Test simulations filtered out cases with runtime errors, solution divergence, or excessive runtime. The final dataset contains 49205 cases that link natural language descriptions to executable OpenFOAM configurations, each including a problem description, mesh files, OpenFOAM input files, an execution script, and an LLM-generated CoT reasoning trace.

The fine-tuned LLM serves as an interface between users and OpenFOAM, translating natural language specifications into executable CFD configurations without requiring users to master OpenFOAM's syntax and parameter structures.

2.2. Multi-agent system

We propose CFD-copilot, an integrated and multi-agent framework designed to automate the CFD workflow from simulation setup to post-processing. The initial workflow employs a self-correcting loop of four agents: pre-checker, generator, runner, and corrector to translate a user's natural language prompt into a viable simulation. The pre-checker validates input and mesh files, enabling the LLM-powered generator to produce OpenFOAM configurations. Subsequently, the runner executes the simulation, while the corrector analyzes error logs and iteratively prompts for adjustments until successful completion. This robust loop ensures a high success rate for the automated simulation setup.

Following simulation, the framework transitions to the post-processing stage, which is orchestrated by the post-processor agent. This post-processor interprets natural language queries from the user, such as requests for data analysis or visualization. The agent leverages MCP, a model-agnostic architecture designed to decouple LLM reasoning from software execution. This approach enables a flexible and extensible framework for post-processing.

3. Results

3.1. Experimental setup

The framework's performance was evaluated on two canonical aerodynamic benchmarks: the NACA 0012 airfoil [4] and the three-element 30P-30N high-lift configuration [5]. To ensure a valid assessment of the model's generalization capabilities, these cases were outside the fine-tuning dataset.

3.2. Result analysis

Our validations confirms the efficacy of the proposed workflow. For the NACA 0012 airfoil case, the system achieved an average success rate of 52.86%, with velocity and pressure accuracies of 96.41% and 93.22%, respectively. The detailed automation performance across different angles of attack (AoA), including success rate, correction iteration counts and token consumption. For the more complex multi-element 30P-30N airfoil, the system demonstrated scalability to realistic high-lift configurations. It achieved a 10% success rate, with velocity and pressure accuracies of 93.14% and 88.44%. In direct contrast, larger, general-purpose models failed to produce a single converged solution, highlighting that domain adaptation is critical for handling complex simulations.

We also evaluated the framework's post-processing capabilities. This process is orchestrated by MCP, which interprets natural language prompts to first select and execute the appropriate OpenFOAM utility for data generation.

We adopt the following prompt for the qualitative analysis of flow fields based on streamline visualization: "Please generate streamline by sampling p and U fields along x axis. Take 1000 sampling points along the line from (-5, -20, 0) to (-5, 20, 0).". The resulting visualization correctly depicts the flow accelerating over the suction surface and forming a wake downstream, providing qualitative confirmation of the underlying flow physics.

We adopt the following prompt to yield the pressure distribution on the main surface of the three-element 30P-30N airfoil: "Please sample field p on the 'wall_airfoil' patches.". After extracting data, a single, general command, "Please write a Python script to draw a scatter plot of normalized chord length and pressure coefficient.", instructed the LLM to generate a script that automatically located, processed, and plotted the data. The resulting distributions show the characteristic features of a high-lift system, including the pressure recovery on the main surface. The results match well with experiments [5]. This demonstrates the ability to manage and visualize data from multiple geometric components in one workflow.

4. Conclusions

We presented CFD-copilot, a natural language-driven framework for end-to-end automation of CFD simulations. The implementation for the initial case setup is managed by a multi-agent, self-correcting system organized around a central generator agent. For the subsequent post-processing phase, the framework employs a modular client-server architecture built on MCP, which allows an LLM-driven client to orchestrate a server hosting a comprehensive suite of independent OpenFOAM utilities. The performance of this integrated, language-driven system was evaluated on two canonical aerodynamic benchmarks, the 2D NACA 0012 airfoil and the three-element 30P-30N airfoil.

References

[1] Z. Dong, Z. Lu, and Y. Yang. Fine-tuning a large language model for automating computational fluid dynamics simulations. Theor. Appl. Mech. Lett., 15:100594, 2025. [2] Z. Dong, S. Du, Z. Lu, and Y. Yang. CFD-copilot: Leveraging domain-adapted large language model and model context protocol to enhance simulation automation, 2025. arXiv:2512.07917. [3] A. Yang, A. Li, B. Yang, et al. Qwen3 technical report, 2025. [4] C. L. Ladson. Effects of independent variation of Mach and Reynolds numbers on the low-speed aerodynamic characteristics of the NACA 0012 airfoil section. National Aeronautics and Space Administration, Scientific and Technical Information Division, Washington, DC, 1988. [5] S. M. Klausmeyer and J. C. Lin. Comparative results from a CFD challenge over a 2D three-element high-lift airfoil. Technical report, 1997.

ABSTRACT. High-speed air-breathing propulsion systems, such as scramjets and ramjets, rely heavily on efficient fuel-air mixing within extremely short residence times. Introducing swirl into the fuel injection stream is a promising technique to enhance mixing efficiency, but it inherently alters the shock structure, leading to complex three-dimensional flow phenomena. In particular, underexpanded jets with swirl exhibit highly distorted shock patterns that are difficult to predict using conventional theories. Therefore, establishing a reliable computational framework to capture these physics is crucial for developing advanced active flow control strategies.This study numerically investigates the detailed shock structure and flow characteristics of underexpanded supersonic nozzle flows under the influence of swirl. The experimental configuration of Abdelhafez is adopted as a reference case. Three-dimensional Reynolds-Averaged Navier-Stokes (RANS) simulations with the Shear Stress Transport (SST) k-omega turbulence model are employed to analyze the complex flow physics, focusing on the modification of the Mach disk location and shock cell structure.The computational results successfully capture the key features of the flow field. It is found that the introduction of swirl induces a strong centrifugal force, causing a radial expansion of the jet plume. This expansion leads to an upstream shift and enlargement of the Mach disk compared to the non-swirling baseline case. Furthermore, the centerline pressure analysis reveals that swirl accelerates the decay of pressure oscillations downstream of the Mach disk, indicating enhanced turbulent mixing and faster dissipation of shock strength. The established computational framework serves as a robust baseline for future investigations into plasma-assisted shock vector control and mixing enhancement in supersonic propulsion systems.

ABSTRACT. 1. Introduction The integration of Unmanned Aerial Vehicles (UAVs) into the civil airspace represents one of the most significant paradigm shifts in modern aviation [1]. The acoustic footprint of small-scale propellers remains a major obstacle to their widespread use around humans and animals. Drones have several noise generation mechanisms : propeller noise, motor noise, and rotor-rotor or rotor-frame interaction. Among these, propeller noise is often the dominant source [2] and thus a viable way to reduce it is to directly optimize the shape of a propeller blade. This aeroacoustic shape optimization problem requires multiple elements: a parametrization of propeller blade shapes, either deformative where a reference blade is deformed to improve performance [3], or generative where new shapes are generated using a model, [4]; a model or simulation to evaluate performance, e.g. Blade Element Momentum Theory (BEMT) [5] or a coupled CFD-CAA (Computational Aeroacoustics) simulation; an optimization method that can handle constraints, multiple objectives, and a high-dimensional design space. In this work, we propose a complete blade optimization pipeline to generate mechanically viable, aerodynamically efficient, low-noise drone propeller blades by using a Generative Adversarial Network-based (GAN) parametrizationcoupled with a fast analytical aeroacoustic simulation based on Farassat’s compact F1A formulation [6] for tonal noise and the Brooks-Pope-Marcolini (BPM) model [7] for broadband noise. The aeroacoustic simulation is validated against experimental data measured in-house on select APC propellers [8].

2. Methodology This section briefly introduces all the components of our shape optimization pipeline, namely the BEMT, noise generation mechanisms, geometry parametrization and optimization problem definition.

2.1 Blade Element Momentum Theory The Blade Element Momentum Theory (BEMT) is a widely used analytical approach for modeling the aerodynamic performance of propellers. The method builds on the fundamental assumption that the momentum lost by the fluid passing through the blade disk must be equal to the loading on the blades. Each blade of a propeller is split into radial sections (i.e. airfoils), and the inflow angle at radius r phi follows tan phi = V_a / wr(1-a') with V_a the section’s axial inflow velocity, w the rotational velocity of the propeller and a' the tangential induction factor (Fig. 1). Once phi is known, the local thrust and torque per unit span are computed and integrated to get the total thrust, torque and efficiency of the propeller [5].

2.2 Aeroacoustic Models The generation of noise by a propeller moving through a fluid is an extremely complex phenomenon comprising two types of noise generation mechanisms. Tonal noise is noise produced at multiples of the blade passing frequency, resulting from the movement of the blade through the air. A stationary observer sees each blade of the propeller get close and then move away repeatedly, producing a distinct acoustic signature. Due to the prohibitive cost of high-resolution CFD simulations, a compact F1A formulation of the Ffowcs-Williams and Hawkings acoustic analogy [6] is used to compute the tonal noise generated by the propeller. This reduces the integration of the pressure and velocity fields over the blade’s surface into a simple line integration over the blade’s span, where each section acts as an acoustic point source. Broadband noise does not have clear harmonics, but instead adds acoustic power to a wide frequency range. It is generated by surface pressure fluctuations caused by turbulent inflow, separation over the blade surface, and trailing edge and tip vortices [7, 9]. Due to the complexity of these phenomena, they cannot be modeled from first principles and require empirical data. In this work, we use the semi-empirical BPM model [7] corrected for source-observer time lag, in conjunction with Amiet’s turbulent inflow noise model [9]. The computed spectrum is validated against measurements of the APC 10x blades (10 inch with variable pitch) [8] carried out in an anechoic chamber.

2.3 Blade Geometry Parametrization Similarly to [10], we parametrize airfoils using a BézierGan with a 4-dimensional latent space, implemented as a WGAN-GP [11] to improve training stability. Two airfoils (root and tip) are sampled from the latent space, and linearly interpolated spanwise for all intermediate sections. A twist distribution theta(r; a, b, c) = a(1 − r)^b + cr^b and a chord distribution c(r; a, b, c, d) = a(1 − r)^b + cr^b − dr ln(r) as a function of radius r are sampled to complete our 15-dimensional blade geometry parametrization.

2.4 Aerodynamic Shape Optimization Aerodynamic Shape Optimization (ASO) is an approach to aerodynamic design that frames the problem of designing a shape to meet certain performance criteria as an optimization problem [12]. In our case, assume that a propeller with a fixed number of blades and a span is parametrized by x. Finding an optimal blade geometry consists of finding x* = min_x f(x) s.t g(x) <= 0, where f is an objective function to be minimized, and g is a constraint function. In this work, the problem is multi objective: • Maximum Figure of Merit FoM = 1/sqrt(2) C_T^1.5/C_Q [13], with C_T, C_Q the coefficients of thrust and torque • Minimum area encompassed by the 60 dB(A) A-weighted Overall Sound Pressure Level (OASPL) contour line around the propeller (Fig. 2). Below Z = 0, the area is weighted 5 times. To be feasible, the thrust should exceed a minimum value T_min = 8 N, and the maximum stress sigma_max cannot exceed sigma_y/K where sigma_y = 150 MPa is the material yield strength and K is a safety factor of 5.

3. Results and Discussion In this section, we validate our aerodynamic and aeroacoustic models against experimental thrust and noise data and present the main results of our optimization.

3.1 Aeroacoustic Model Validation The BEMT, F1A and BPM models are validated in Fig. 3 and Fig. 4. The thrust and first harmonic noise level are captured by our model, but the acoustic prediction is poor at low frequencies and not complex enough at high frequencies. This can be explained by the BPM model assumptions being violated (wing with a known NACA profile in a straight line motion). Nevertheless, most of the acoustic power is generated at the first harmonic, so the OASPL values agree despite error throughout the spectrum.

3.2 Optimization Results To solve the ASO problem, we use SMS-EMOA [14]. The main idea of this evolutionary algorithm is to repeatedly prune points that least contribute to the hypervolume. Given the 1-1.5s problem evaluation time (Intel i7-13700K, Nvidia RTX A4000), sample efficiency is secondary, and an evolutionary algorithm is preferred to Bayesian optimization due to minimal computational overhead between rounds. SMS-EMOA has a population of 20 points that evolve over 100 generations. Figures 5 and 6 show the Pareto fronts and hypervolumes over 5 separate runs. Figure 7 shows a Pareto-optimal blade found using this algorithm. The corresponding 60 dB(A) contour plot is shown in Fig. 8. This propeller is predicted to produce 8.001 N of thrust and experience a maximum stress of 28.5 MPa, which is within our feasible set. The very thin cambered airfoils and large chord help reduce noise generation while ensuring sufficient thrust.

ABSTRACT. Voronoi Tessellation-assisted CNN–PINNs with End-to-End Sensor Placement Optimization for Robust Flow Reconstruction Renjie Xiaoa,1, Xiaoping Tiana,1, Bingteng Suna,*, Hao Yina, Qiang Dua,b,c,d,*, Junqiang Zhua,b,c a. Advanced Gas Turbine Laboratory, Institute of Engineering Thermophysics, Chinese Academy of Sciences, Beijing, 100190, China b. Qingdao Institute of Aeronautical Technology c. Nanjing Future Energy System Research Institute, Nanjing, Jiangsu 211135, China d. National Key Laboratory of Science and Technology on Advanced Light-duty Gas-turbine, Beijing 100190, China e. University of Chinese Academy of Sciences, Beijing, 100190, China

Keywords: Physics-informed neural networks, Flow reconstruction, Sensor placement optimization, Voronoi tessellation, Irregular domains.

Introduction Fluid physics is traditionally investigated by high-fidelity simulation and laboratory measurements. Yet direct numerical simulation remains prohibitively expensive for many practical Reynolds numbers [1], and time-resolved flow visualization can demand substantial experimental resources [2]. Flow field reconstruction therefore aims to infer the full state from a limited set of measurements, enabling digital twins and closed-loop control when dense sensing is infeasible [3]. Recent progress in data-driven reconstruction can be viewed through three methodological families: direct supervised mappings, regression-based methods, and data assimilation, as summarized by Dubois et al. [4]. With the rise of modern machine learning for fluid mechanics [5], super-resolution and sparse-to-dense recovery have become a focal direction, including Voronoi tessellation-assisted reconstruction that improves robustness to irregular sensor layouts [6]. However, most pipelines still assume a fixed number of sensors and a stable layout, which limits their deploy ability under sensor dropouts and changing geometries. Physics-informed neural networks (PINNs) complement purely data-driven models by embedding governing equations into the learning objective. Representative variants such as NSFnets [7] and PINN-SR [8] have demonstrated improved data efficiency for incompressible flows. For irregular domains, geometry-adaptive operators such as PhyGeoNet [9] alleviate boundary discretization artifacts, while physics-informed sensor placement has also been explored in related PINN settings [10]. Building on these threads, we propose VSOPINN, a Voronoi tessellation-assisted CNN–PINN framework that jointly reconstructs flow fields and optimizes sensor locations end-to-end.

* Correspondingauthor Email address: sunbingteng@iet.cn (Bingteng Sun) Methodology Our approach, termed VSOPINN, couples three ideas: Voronoi image lifting of sparse data, geometry-adaptive coordinate transformation, and end-to-end sensor placement optimization. Figure 1 summarizes the overall workflow. First, given sensor locations and their measurements, we construct a Voronoi tessellation and rasterize it on a structured reference grid. Each grid point inherits the value of its nearest sensor, producing an image-like representation that preserves the measurement topology while remaining compatible with convolutional feature extraction. To improve smoothness and reduce degeneracy when sensors move, we further regularize the tessellation via centroidal Voronoi tessellation (CVT), where sensor sites are encouraged to coincide with weighted cell centroids. The weight (density) is derived from the local physics residual so that high-error regions attract more sensing capacity. Second, to support irregular physical geometries without staircase boundary artifacts, we adopt a geometry-adaptive mapping that transforms an irregular domain into a regular computational domain, enabling PDE derivatives to be evaluated through standard CNN operations. This idea follows geometry-adaptive convolutional operators as in PhyGeoNet [9], but is embedded here within a physics-informed encoder–decoder network. Third, VSOPINN enforces the incompressible Navier–Stokes equations through a composite loss combining sensor data misfit, interior PDE residuals, and boundary-condition residuals, while treating sensor coordinates as trainable variables. During training, sensor locations are updated jointly with network weights using gradient-based optimization, augmented by CVT regularization to maintain well-distributed sensors and to stabilize the Voronoi lifting. For incompressible flows, the PDE residuals are constructed from the continuity and momentum equations: \mathrm{\nabla}\centerdot u\ =\ 0 \rho(\partial u/\partial t\ +\ u\centerdot\mathrm{\nabla u})\ +\ \mathrm{\nabla p}\ -\ \mu\nabla²u = 0 The overall training objective is a weighted sum of data, physics, boundary, and CVT regularization terms: L\ =\ L_{data}\ +\ \lambda_f\ L_{PDE}+\ \lambda_{bc}\ L_{bc}\ +\ \lambda_{cvt}\ L_{cvt}

Figure 1: Diagram of the VSOPINN model. (a) Encoder/Decoder modules, (b) Attention- enhanced convolutional decoder, (c) Finite difference convolution kernel. Conclusions Across diverse benchmarks, VSOPINN consistently improves reconstruction accuracy over geometry-adaptive baselines under tight sensing budgets. Figure 2 visualizes the reconstructed flow fields compared against the CFD reference. For the lid-driven cavity at Re=100 using only four sensors, the relative L2 error of the velocity magnitude decreases from 4.89×10−1 (physics-only baseline) to 2.66×10−1 with Voronoi lifting, and further to 3.89×10−2 with CVT-guided end-to-end sensor optimization. In the annulus driving flow, residual-guided sensor migration reduces the error from 6.39×10−1 to 7.22×10−2, indicating that the learned layout concentrates sensing capacity near steep shear layers. Overall, the Voronoi-regularized CNN–PINN framework offers a practical pipeline for reconstructing velocity and pressure fields from sparse, irregular measurements in complex geometries while simultaneously designing an informative sensor layout. By coupling Voronoi image lifting, geometry-adaptive mapping, and CVT-guided end-to-end optimization, the method improves accuracy and robustness under sensor-layout changes. Future work will target scalable 3D formulations and integration with real-time sensing for closed-loop experiments.

Figure 2: Visual comparison of the velocity magnitude and pressure fields for the lid-driven cavity flow at Re = 100. The columns correspond to the CFD reference and the predictions from the progressive ablation models. The optimized sensor locations in VSOPINN with CVT enable a more accurate reconstruction of the primary vortex structure and corner singularities. References [1] P. G. Huang, G. N. Coleman, and P. Bradshaw, Compressible turbulent channel flows: DNS results and modelling, Journal of Fluid Mechanics, 305, 185-218, 1995. [2] Zhiwen Deng, Yujia Chen, Yingzheng Liu, and Kyung Chun Kim, Time-resolved turbulent velocity field reconstruction using a long short-term memory (LSTM)-based artificial intelligence framework, Physics of Fluids, 31(7), 2019. [3] Jared L. Callaham, Kazuki Maeda, and Steven L. Brunton, Robust flow reconstruction from limited measurements via sparse representation, Physical Review Fluids, 4(10), 2019. [4] Pierre Dubois, Thomas Gomez, Laurent Planckaert, and Laurent Perret, Machine learning for fluid flow reconstruction from limited measurements, Journal of Computational Physics, 448, 110733, 2022. [5] Steven L. Brunton, Bernd R. Noack, and Petros Koumoutsakos, Machine Learning for Fluid Mechanics, Annual Review of Fluid Mechanics, 52(1), 477-508, 2020. [6] Kai Fukami, Romit Maulik, Nesar Ramachandra, Koji Fukagata, and Kunihiko Taira, Global field reconstruction from sparse sensors with Voronoi tessellation-assisted deep learning, Nature Machine Intelligence, 3(11), 945-951, 2021. [7] Xiaowei Jin, Shengze Cai, Hui Li, and George Em Karniadakis, NSFnets (Navier-Stokes flow nets): Physics-informed neural networks for the incompressible Navier-Stokes equations, Journal of Computational Physics, 426, 109951, 2021. [8] Zhao Chen, Yang Liu, and Hao Sun, Physics-informed learning of governing equations from scarce data, Nature Communications, 12(1), 2021. [9] Han Gao, Luning Sun, and Jian-Xun Wang, PhyGeoNet: Physics-informed geometry-adaptive convolutional neural networks for solving parameterized steady-state PDEs on irregular domain, Journal of Computational Physics, 428, 110079, 2021. [10] Shengze Cai, Zhicheng Wang, Sifan Wang, Paris Perdikaris, and George Em Karniadakis, Physics-Informed Neural Networks for Heat Transfer Problems, Journal of Heat Transfer, 143(6), 2021.

ABSTRACT. Thermal runaway in lithium-ion batteries poses critical safety challenges, while conventional CFD-based analyses remain computationally expensive and unsuitable for real-time prediction. This study proposes an extended physics-aware recurrent convolutional neural network (PARC) for rapid and accurate prediction of battery thermal runaway. The original PARC architecture is reformulated to model diffusion–reaction dynamics coupled with Joule and reversible heating, incorporating a generalized Arrhenius-based reaction term to represent exothermic decomposition processes. The model predicts the spatiotemporal evolution of temperature, electrolyte concentration, and voltage using only initial and boundary conditions. Validation against CFD datasets under diverse operating conditions demonstrates that the proposed framework accurately captures both normal thermal behavior and the onset of thermal runaway. The trained model achieves several orders-of-magnitude speedup compared to CFD, highlighting its potential for real-time battery safety assessment.

ABSTRACT. When a high-speed boundary layer encounters a rectangular cavity, the flow separates at the leading edge and forms a shear layer spanning the cavity opening that reattaches near the trailing edge, resulting in an open cavity flow. Open cavity flows can occur in a variety of situations, such as weapon bays, scramjet combustors, etc. The shear layer associated with the separated flow in the open cavity is inherently unstable and establishes a self-sustained oscillatory mechanism in which vortical structures shed from the upstream cavity edge convect downstream, impinge on the trailing edge, and generate upstream-travelling acoustic waves that trigger new vortex formation at the leading edge. This feedback process produces discrete resonant frequencies known as Rossiter modes. The characteristics of these oscillations are strongly influenced by geometric parameters such as cavity length, depth, and wall heights. Variation in rear-wall height can significantly alter the shear-layer dynamics and the associated feedback mechanism. Configurations with relatively higher rear walls exhibit low-frequency, high-amplitude oscillations, whereas configurations with lower rear walls tend to be dominated by higher-frequency modes. The influence of rear-wall height in high-supersonic regimes remains insufficiently understood. The present work investigates the effect of rear-wall height on the unsteady behaviour of rectangular cavity flows at a freestream Mach number of 4.0. Four configurations with different rear-wall heights are considered for a fixed cavity length-to-depth ratio of L/D = 6. High-order two-dimensional simulations of the compressible Navier-Stokes equations are performed using the GPU-accelerated OpenSBLI framework. The resulting unsteady flow fields are analysed using a spectral proper orthogonal decomposition to identify dominant oscillation frequencies and associated coherent structures. The results indicate that increasing rear-wall height shifts the dominant response toward lower-frequency Rossiter modes and eventually leads to a new low-frequency mode that is different from Rossiter modes, that remains relatively unexplored. Understanding the effect of the rear-wall height can shed light on the unsteady dynamics that can arise when the cavity design and geometry are varied in applications (e.g., weapons bays, scramjet combustors, etc.) and also pave the way for new flow control approaches for mitigating flow oscillations in open cavity flows.

ABSTRACT. This study develops a skeletal chemical kinetic mechanism for Rocket Propellant-3 (RP-3) adopting a three-component RP-3 surrogate, which consists of n-dodecane, 1,3,5-trimethylcyclohexane and n-propylbenzene, thus integrating it into a computational fluid dynamics (CFD) model for scramjet combustor simulation. The detailed RP-3 mechanism with 446 species and 2,942 reactions is reduced by a hybrid reduction strategy combining path-flux analysis (PFA) and an artificial neural network (ANN), which yields a skeletal RP-3 mechanism containing 113 species and 798 reactions. The reduced mechanism is validated against ignition delay times, species profiles in a jet-stirred reactor and laminar flame speeds over relevant thermochemical conditions. The skeletal mechanism is then applied in a two-dimensional (2-D) scramjet engine model with single- and dual-cavity configurations to investigate the combustion characteristics in kerosene-fuelled supersonic flow using hydrogen as pilot. The simulation with shear stress transport (SST) k-ω turbulence model reproduces experimentally wall-pressure profiles for various inlet kerosene temperature, capturing both the magnitudes and the ordering of pressure peaks for both the single- and dual-cavity cases. This CFD model gives insights into the influences of inlet fuel temperature and cavity geometry on combustion behaviors and flow structures, and a pathway analysis is presented to reveal the routes by which RP-3 forms its principal thermal-cracking products. Overall, the proposed skeletal mechanism balances fidelity and computational cost, enabling tractable CFD studies of kerosene-fuelled scramjet combustion and supporting design-oriented exploration of operating conditions and cavity geometries.

ABSTRACT. Hydrogen micromix combustion offers rapid fuel-air mixing and flashback mitigation for low-carbon gas turbines, yet the influence of air-hole staggering on flow organization, mixing development, reaction-zone structure and NOx response remains insufficiently clarified. This study numerically investigates a confined hydrogen micromix unit at a lean equivalence ratio of 0.381, considering four circumferential offsets between two rows of air holes: 0°, 10°, 20° and 30°. Time-averaged velocity fields, mixing uniformity, instantaneous temperature evolution, reaction-zone statistics and radial NOx distributions are analyzed. Results show that air-hole staggering modifies the near-field high-speed passage, shear layer and local recirculation, thereby enhancing circumferential transport and dilution of the hydrogen jet. Compared with the aligned case, staggered configurations improve near-field mixing uniformity; Case 10° provides faster mixing within a limited passage length, whereas Case 20° gives a more balanced performance after further downstream development. In the combustion zone, the flame remains globally anchored near the micromixer outlet, while the flame root exhibits local unsteady wrinkling due to shear-layer entrainment and hot-product recirculation. The stagger angle further changes the high-temperature area fraction, the extent of strongly reactive regions and the radial NOx distribution. Overall, moderate air-hole staggering provides a practical compromise among rapid mixing, stable heat release and spatial control of NOx formation.

ABSTRACT. This study presents a numerical investigation of unsteady flow fields surrounding the Caradonna-Tung rotor, serving as a critical validation baseline for rotorcraft aerodynamic simulations. Utilizing an experimental test case at 1500 RPM and an 8 degrees collective pitch, the research employs a multiblock grid structure framework integrated with a sliding interface to resolve the unsteady transient rotational flow physics. The numerical results demonstrate a strong correlation with experimental surface pressure distributions (Cp) and thrust coefficient (Ct), confirming the accuracy of the unsteady flow solver and the effectiveness of the grid structure. The study provides the necessary foundation for future high-fidelity studies involving fully coupled fluid-structure interaction (FSI) and the implementation of deforming mesh blocks for comprehensive aeroelastic analysis.

ABSTRACT. A scalable numerical framework for the interaction of incompressible fluids and structures (FSI) involving rigid bodies has been developed by coupling a sharp-interface immersed boundary formulation with a multigrid-accelerated fractional-step projection method. The solver is implemented on distributed GPU platforms to enable large-scale simulations, while the underlying numerical formulation remains general and independent of hardware-specific assumptions.

The accuracy and robustness of the proposed framework were verified through a hierarchy of benchmark problems. In the three-dimensional lid-driven cavity flow with an embedded rigid sphere, the solver accurately reproduced the modification of the primary recirculation structure and its systematic evolution with increasing Reynolds number. The framework was also applied to large-scale gravitational settling simulations involving up to 20,000 fully resolved particles, revealing complex flow instabilities and coherent vortex structures in dense particulate systems.

ABSTRACT. In nuclear reactors, the flow of coolant plays an important role in the heat transfer of fuel rod bundle, and the characteristics of the turbulence have a significant impact on the heat transfer and the safety of the fuel rods. In this study, the large eddy simulation (LES) data obtained previously were first employed to perform Proper Orthogonal Decomposition (POD) analysis on the fuel rod bundle flow field, encompassing modal energy analysis and flow field reconstruction to enable rapid flow field reconstruction. Subsequently, a Temporal Convolutional Network (TCN) was integrated to predict the evolution of POD modal coefficients, thereby achieving fast flow field prediction. This hybrid POD-TCN framework facilitates efficient fluid dynamic analysis of the flow field while significantly reducing computational costs. The modal energy spectrum indicates a relatively low energy contribution for each individual mode, while the cumulative energy of the first six modes exceeds 65%. At Re is 11164 condition, using only the first 20 modes achieves a reconstruction accuracy greater than 99%, with errors mainly concentrated in near-wall regions and around rod surfaces where small-scale vortices dominate. In addition, the predicted-reconstructed flow field closely matches the true field in terms of vortex structures and velocity gradients.

ABSTRACT. 1. Introduction Floating offshore wind turbines (FOWTs) operate under platform motions induced by wind, waves, and coupled aero-hydro effects. In a wind turbine, the root vortex generates counter-rotating motion in the wake, while tip vorticity shed from the blade boundary layer forms a coherent helical structure that propagates downstream [1]. Unlike fixed-bottom turbines, floating rotors undergo rigid-body motion relative to their wake, fundamentally altering the stability and structure of the blade-generated vortex system. Surge motion periodically displaces the rotor through its own wake, promoting vortex interaction, pairing, and breakdown, and thereby modifying turbulence characteristics in both the near and far wake. Instability becomes significant when adjacent helical vortex filaments interact strongly and undergo vortex pairing, preceding distortion and breakdown and enhancing turbulent mixing and kinetic-energy transport across the shear layer [2]. The dynamics of vortex ring pairs are governed by collective convection and self-induced velocities arising from their circulation [3]. Opposite-signed circulations induce separating velocities, while mutual interaction promotes paired motion; the balance between these mechanisms determines vortex stability until diffusion or instability disrupts the structure [3]. High-fidelity blade-resolved CFD can capture these unsteady dynamics but becomes computationally prohibitive under six-degree-of-freedom motion. The Actuator Line Model (ALM) provides an efficient alternative by representing blade forces as distributed momentum sources, capturing large-scale wake structures and tip vortex evolution, although its accuracy under motion-induced inflow remains sensitive to force projection and filtering strategies. In this work, the unsteady wake dynamics of a floating wind turbine undergoing prescribed surge motion are investigated using a three-dimensional Actuator Line Model (ALM). The actuator line model is developed by implementing a source term in the momentum equations. Dedicated UDFs (User Defined Functions) are developed to introduce the actuator line model for blade and rotor motions. A tip loss correction technique is applied to account for the effect of the blade tip vortex. Preliminary results show that ALM is capable of capturing vortex trajectories and wake meandering in FOWTs during surge motion. The simulations are carried out using the Ansys software package through implementation of UDFs that introduce the actuator line model for blade and rotor motions. 2. Methodology The flow is governed by the incompressible, unsteady Navier-Stokes equations with additional body-force source terms representing the actuator. Flow field simulations are carried out using incompressible URANS with a pressure-based solution approach. All simulations are performed in Ansys Fluent using a pressure-based transient solver with second-order temporal and spatial discretization. 2.1 Actuator Line Model Formulation The actuator line model (ALM) represents turbine blades as rotating lines of momentum sources instead of resolving each blade geometry inside the flow domain, applied over a projection width. The source term S_ALM on the right-hand side of Equation (1) represents the body force that artificially represents the blade section in the actuator line model. The aerodynamic forces are calculated and distributed onto the computational mesh. These body forces are obtained from local flow conditions and airfoil characteristics and are spatially smeared using a regularization kernel. Although the blade geometry and boundary layers are not explicitly resolved, the ALM allows generation of blade-induced vorticity as well as tip vortices, enabling near and far wake structures to be captured at reduced computational cost compared to blade-resolved simulations [5]. In the ALM technique, each blade is represented by a rotating line discretized into spanwise segments. Each segment has a constant airfoil shape, twist, chord, and incoming wind kinematics. Velocities are calculated locally for each blade element. Blade pitch and sectional lift and drag forces are evaluated using pre-tabulated airfoil polars. The sectional lift and drag forces are projected onto the Cartesian momentum equations and distributed to the surrounding flow field using a three-dimensional Gaussian kernel. The aerodynamic forces are smoothly projected onto the flow field. The projection width of the Gaussian operation is determined by epsilon, chosen proportional to the local grid spacing to ensure numerical stability while preserving vortex coherence. In the current results, epsilon is set to 3 times the mesh spacing in the streamwise direction, that is, epsilon approximately equals 3 delta x. The Gaussian smearing prevents singular forcing and enables smooth coupling between actuator lines and the resolved flow field. Airfoil polars are introduced in the spanwise direction of the rotor, and aerodynamic forces generated by the turbine are computed and distributed as body forces along the line elements representing the blades. This approach allows implementation of aerodynamic forces and corresponding changes without fully resolving blade geometry, resulting in significant reduction in computational cost while still enabling instantaneous vortical structures in the wake to be resolved and detailed investigation of vortex interactions during surge motion of a floating wind turbine [2]. 2.2 Platform Surge Motion Modelling Platform surge motion is prescribed as a harmonic displacement in the streamwise direction of the rotor. The surge motion is incorporated directly into actuator-line kinematics by updating the rotor center position and relative velocity at each time step. Since blade geometry is replaced by the actuator line model, no physical blade geometry exists in the computational domain. Therefore, blade motion is represented by modifying actuator line location and inflow velocities through a UDF. 3. Results These preliminary results are presented for a generic rotor undergoing surge motion. The purpose is to investigate the feasibility of implementing a simplified actuator line model as a low-fidelity approach for capturing formation and development of vortex structures in the wake of a floating wind turbine under surge motion. The computational domain is a cuboid of dimensions 10D, 3D, and 3D in streamwise, spanwise, and normal directions. The distances from inlet and outlet to the rotor center are 3D and 10D, respectively. The mesh is refined around the actuator disk and wake region by introducing two concentric cylinders at the domain center. The total number of cells is 626572. Simulations are carried out at a tip speed ratio lambda equal to 1.66 and zero pitch angle. The harmonic surge motion is defined with reduced frequency F_star equal to 0.9 and reduced amplitude A_star equal to 0.1. Reduced frequency and amplitude are normalized using rotor rotational frequency and rotor diameter, respectively. The time step is delta t equal to 1 divided by (40 times Omega), where Omega is the rotational speed of the rotor. Within the UDF, lift and drag coefficients are defined using thin airfoil theory [6], and source terms are calculated using local aerodynamic forces. 3.1 Wake Structure and Vortex Interaction Mechanisms In the absence of platform motion, the wake is characterized by a stable helical vortex system dominated by coherent tip vortices. Under surge motion, relative displacement between rotor and wake introduces periodic modulation of vortex spacing and strength. As the rotor advances into regions of reduced axial velocity, adjacent helical filaments experience enhanced mutual induction, leading to vortex collisions and pairing. In the present case, surge-induced unsteadiness accelerates vortex distortion and breakdown, promoting enhanced turbulent kinetic energy production. The onset and intensity of vortex interaction depend strongly on surge amplitude and frequency, indicating critical coupling between platform dynamics and wake stability. In the final version, results will be improved in terms of mesh resolution and numerical accuracy and validated through comparison with results for the NREL 5 MW turbine under surge motion. Additional vortex shedding and collision simulations will be performed to induce vortex displacement and collisions, enabling analysis of instability development, structural changes, and their role in wake instability growth. Conclusions The actuator line model has been implemented for simulation of a floating wind turbine undergoing harmonic surge motion. Results demonstrate that platform surge motion fundamentally alters vortex dynamics and wake stability in floating wind turbines. Further simulations in the final submission will expand this framework, allowing visualization of wake breakdown and mixing processes under varying surge amplitudes and frequencies to reveal how platform oscillations interact with different turbulent regimes. 4. References [1] S. Ivanell, R. Mikkelsen, J. N. Sorensen, and D. Henningson, "Stability analysis of the tip vortices of a wind turbine," Wind Energy, vol. 13, no. 8, pp. 705-715, 2010. [2] T. Li, Y. Zhang, Q. Yang, X. Zhou, Z. Zhang, and T. Wang, "Unsteady aerodynamic characteristics of a floating offshore wind turbine in propeller state," Renewable Energy, vol. 246, p. 122861, 2025. [3] P. Weidman and N. Riley, "Vortex ring pairs: numerical simulation and experiment," Journal of Fluid Mechanics, vol. 257, pp. 311-337, 1993. [4] L. Tophoj and H. Aref, "Instability of vortex pair leapfrogging," Physics of Fluids, vol. 25, no. 1, 2013. [5] D. Zhao, N. Han, E. Goh, J. Cater, and A. Reinecke, Wind turbines and aerodynamics energy harvesters. Academic Press, 2019. [6] W. Zhong, W. Z. Shen, T. Wang, and Y. Li, "A tip loss correction model for wind turbine aerodynamic performance prediction," Renewable Energy, vol. 147, pp. 223-238, 2020.

ABSTRACT. Wall-resolved large-eddy simulation (WRLES) has traditionally been considered computationally prohibitive for engineering applications. This study demonstrates affordable aeroacoustic prediction of high-lift device slat noise using high-order WRLES coupled with the Ffowcs Williams-Hawkings (FW-H) acoustic analogy. The space-time expansion of high-order kinetic preserving flux reconstruction scheme (STE-KEP-FR) enables efficient computation, with only a few GPUs are required to perform WRLES simulations. Two high Reynolds number benchmarks are examined: the 30P30N high-lift airfoil and the Trap Wing semi-span model. These configurations represent 2.5D and 3D slat noise problems with degrees of freedom (DoFs) ranging from 31.3 million to 0.32 billion. Results demonstrate excellent agreement with experimental data. The 30P30N benchmark case demonstrates that our high-order solver Dimaxer achieves two orders of magnitude higher accuracy and computational efficiency compared to the CPU-based second-order finite volume solver OpenCFD-EC.