ABSTRACT. The equations of continuum mechanics encompass both fluids and solids, covering a broad spectrum of applications in engineering and physics. Within the framework of hyperbolic thermodynamically compatible (HTC) models, these equations can be unified into a single mathematical formulation. This model admits an additional conservation law for the total energy, in accordance with the first principle of thermodynamics, and satisfies the entropy inequality dictated by the second principle.

Moreover, in the absence of algebraic source terms, both the distortion field and the specific thermal impulse remain curl-free, provided that the initial data satisfy this condition. The model is further constrained by a nonlinear algebraic relation linking the determinant of the distortion field to the density of the medium. In the stiff relaxation limit, the system asymptotically recovers the compressible Navier–Stokes equations, thereby enabling a unified description of nonlinear solid and fluid dynamics.

In this talk, I will present and analyze innovative semi-discrete thermodynamically compatible finite volume schemes on unstructured meshes, designed to preserve most of these structural properties exactly at the discrete level. I will also briefly discuss asymptotic preservation and introduce a novel approach to ensure geometric compatibility on both fixed and moving unstructured meshes.

ABSTRACT. The phase-change, e.g. condensation, phenomena have been a hot issue in the field of gas-liquid two-phase flow researches. To predict the occurrence of phase-change phenomena accurately in gas-liquid two-phase flows, not only the advanced measurement techniques, but the credible CFD methods should be established. In this study, the authors propose a modified condensation simulation procedure based on an interface-tracking method, in consideration of interface advancement thorough a cell-boundary. The proposed simulation method is verified by solving the well-known Stefan problem to show the superiority to conventional methods. Then, the proposed method is employed to the simulation of a simple experiment of direct contact condensation phenomena to show the applicability to practical complicated condensation flows.

ABSTRACT. Mass transfer across resolved interfaces in two-phase flows is a key mechanism in many industrial and environmental applications, such as gas–liquid reactors, spray systems, and pollutant removal processes. Accurately capturing interfacial mass transfer remains a major numerical challenge due to thin diffusion layers and the stiffness induced by thermodynamic discontinuities of concentration at the interface. Classical one-fluid formulations rely on mixture variables and additional source terms, which may limit their accuracy when strong interfacial jumps are present. In this context, the present work investigates the ability of a discontinuous two-field formulation combined with the Ghost Fluid Method to faithfully represent sharp interfacial physics on fixed structured grids.

The proposed approach is based on a two-field framework for the resolution of mass transfer with explicitly resolved interfaces. In this framework, grid cut-cells are fully populated with two fields: one is directly resolved, while the other is reconstructed using different strategies. In purely single-phase cells, separate advection–diffusion equations are solved in each phase, whereas interfacial diffusive fluxes and jump conditions are handled using the Ghost Fluid Method. The interface dynamics are captured using a Volume-of-Fluid method, and a simultaneous and synchronized transport of the concentration field and the volume fraction is handled though Weymouth’s geometric VOF scheme. The method is implemented in an in-house finite-volume solver and systematically compared to a one-fluid model incorporating an interfacial source term. Benchmark configurations with analytical solutions, as well as more realistic configurations involving moving and deformable interfaces, are used to assess accuracy, convergence, and numerical robustness.

The results demonstrate the clear superiority of the discontinuous two-field/Ghost Fluid approach for sharp mass transfer problems. This assessment is supported by systematic convergence studies conducted on both the concentration field and the Sherwood number, using analytical reference solutions encompassing static and moving interface scenarios. Various configurations representative of the three mass transfer regimes (internal, external and conjugate) are investigated. While a slightly higher order of convergence is observed, the most significant improvement lies in the accuracy, with errors often reduced by up to one order of magnitude compared to the one-fluid formulation. In addition, an analysis of temporal fluctuations of key physical quantities during interface crossing highlights the improved stability and consistency of the proposed advection strategy. Finally, the approach is successfully applied to a three-dimensional convection–diffusion problem involving a deformable bubble, illustrating its potential for realistic simulations where no analytical solutions are available and paving the way for future extensions to complex multiphase configurations.

ABSTRACT. Numerical simulations of multiphase flows are of great importance in practical applications, including ship hydrodynamics, naval architecture and ocean engineering. These simulations typically involve strong topological deformations of free surfaces such as merging and fragmentation, sharp changes in fluid properties such as density and viscosity, and body-fitted unstructured grids to comply with the complex geometry shape of ships and offshore structures, which pose numerous challenges to the application of computational fluid dynamics (CFD). On the one hand, the large fluid density ratio will generate severe spurious currents on the non-orthogonal unstructured grids if an imbalanced discretization is employed between the pressure gradient and the external forces. On the other hand, the accuracy of velocity field will deteriorate or even lead to divergence if a merely first-order accurate scheme is employed to reconstruct the velocity variation. To address the above issues, a unified formulation of the generic balanced-force algorithm and the second-order velocity reconstruction scheme, named as UGBF-VR, is proposed in the finite volume framework to develop a high-fidelity numerical model for multiphase flows simulations on polyhedral unstructured grids. Compared with the conventional formulation, it simultaneously increases the accuracy of velocity field in the entire domain and suppresses the spurious velocity generated near the fluid interface, which offers a promising platform to provide more accurate and robust predictions for incompressible multiphase flows with large density ratio.

ABSTRACT. Solid-liquid phase change phenomena are encountered in a variety of science and engineering problems. For instance, latent heat thermal energy storage, transient thermal management of electronic equipment, and additive manufacturing. The enthalpy-porosity method is one of the commonly used methods for simulation of solid-liquid phase change problems. However, existing formulations of the enthalpy- porosity method rely on the assumption of incompressible flow and typically utilize semi-implicit or implicit time integration methods. In this work, we propose a novel method based on operator splitting between the advection and the diffusion parts, where we use the Super-Time-Stepping method (STS) for integration of the viscous steps, explicit integration of the advection and enthalpy method for the phase transition with improved accuracy and larger time steps. We compare our preliminary results to analytical solutions of the Stefan problem with advection.

ABSTRACT. Topological defects profoundly influence flow behavior and orientational order in active nematics, making the control of such defects important for shaping the dynamics of active nematic systems in both biological and technological contexts. Applied activity patterns can induce self-propulsion of active nematic defects, but methods for achieving general-purpose control over these defects by exploiting this phenomenon have been largely unexplored. We introduce a systematic framework leveraging deep reinforcement learning (DRL) to develop flexible, performant controllers that can apply spatiotemporally varying activity patterns to steer active nematic defects through arbitrarily complex microchannel geometries. We demonstrate the power of our method by training DRL controllers on minimum-time position control tasks in hybrid lattice Boltzmann simulations of active nematodynamics in different microchannel configurations, outperforming static activity patterns and even our own rule-based controller. Moreover, with no additional fine-tuning, the resulting trained controllers can be combined into a single meta-controller capable of successfully steering an individual defect through a test maze of microchannel junctions. Our framework and results unlock new classes of challenging control problems and their solutions in soft and active matter systems, accelerating potential applications in microfluidics, programmable materials, and biophysics.

ABSTRACT. Laminar to turbulent boundary layer transition presents a monumental challenge for the practical imple- mentation of hypersonic flight. Relative to their laminar counterparts, turbulent boundary layers result in an increase in wall heat flux and shear loading. To this end, it is of great importance to the high-speed community that the mechanisms by which the laminar region of the boundary layer may be extended are fully understood so that we may try to delay the transition process. Under low levels of background noise, the driving force of boundary layer transition is the exponential growth of unstable eigenmodes within the laminar region as described by [1], with the dominant instability for a given flow varying based on a surface geometry and Mach number. For hypersonic flows (M > 4.5) with a relatively cold wall boundary condition the dominant instability is the 2D second Mack mode. This disturbance behaves as a trapped acoustic wave between the wall boundary and the sonic line (M = 1) within the boundary layer, propagating downstream with a phase velocity greater than the speed of sound. Laminar flow control (LFC) methodologies generally fall into one of three categories: active, passive, or hybrid control. Active flow control methods utilize an external energy source to impact disturbance growth rate which, while promising, are mechanically complex and difficult to implement. Passive flow control methods on the other hand require no energy input and rely on the design and surface characteristics of the vehicle to delay transition. Hybrid flow control utilizes a combination of both active and passive components. One such proposed method of passive flow control is the use of acoustic metasurfaces which interact with boundary layer disturbances through acoustic waves. This umbrella term ’acoustic metasurfaces’ covers a variety of engineered surfaces including porous media, microstructures, absorptive coatings, among others. Of interest are three flavors: porous media, in which the acoustic energy of the second mode is absorbed/dissipated by the wall boundary, reflection-controlled, in which the reflection angle is altered for an incident wave, and phononic subsurfaces, which we will discuss here shortly. Early analysis of porous media was performed on high-speed boundary layer disturbances in which a model was developed for the acoustic and thermal admittance characteristics of a porous wall boundary condition [2]. These results, along with experimental comparisons, demonstrated strong reduction in second mode disturbance growth rates. The model was subsequently extended to account for higher-order refracted waves and mode interaction between pores, which further demonstrated second mode amplitude reduction [3]. The previous three decades have given rise to a strong interest in the the field of phononic materials, which exist as an engineered periodic medium composed of one or more compounds whose frequency Thirteenth International Conference on Computational Fluid Dynamics (ICCFD13) Milan, Italy, July 06-10, 2026 2 bands exhibit specific bandgaps, allowing for attenuation of mechanical waves. These bandgaps are quite large and highly tailorable, as demonstrated in [4], making phononic materials highly applicable for a variety of applications across a broad range of disciplines. A preliminary study by [5] examined the effect of subsurface phononic materials (termed ’PSubs’) on subsonic Tollmien-Schlichting (TS) waves in a parallel channel flow, with results indicating that for a monochromatic disturbance forcing (single- frequency) the presence of a PSub worked to reduce disturbance kinetic energy and stabilize regions of the flow. The PSub application was extended to spatially-devloping 2D subsonic boundary layer flows [6]. A multi-input multi-output (MIMO) PSub was designed which exhibited a broad range of achievable phase relations. This MIMO PSub implementation resulted in a decrease in the turbulent kinetic energy of the disturbance and a delay in transition onset based on skin friction. The use of PSubs, with their highly-tunable frequency-response properties, as a potential method of disturbance stabilization is an exiting avenue to pursue and shows early potential. This abstract demonstrates successful reduction of boundary layer disturbances (both broadband and narrowband) through the placement of optimized subsurface phononic materials.

2. Methodology The PSub implementation is examined using linear stability theory (LST) and direct numerical simulation (DNS). The former is a simplified parallel-flow approximation in which the eigenvalue problem is solved with a complex impedance-type boundary condition. The latter involve a recasting of the Navier-Stokes equations in which the vector of primitive variables is decomposed into a steady quantity accounting for the mean flowfield and an unsteady quantity accounting for the disturbance flowfield. This formulation is modeled using a fully dynamic boundary condition with two-way coupling between the surface and the disturbance. The mean flow solution is subtracted out and we arrive at the governing equations for small-amplitude disturbances. Upon neglecting nonlinear terms, these so-called linear disturbance equations (LDE) can be written as

The transformation matrix from the conservative to primitive state vector is denoted by long with the inviscid and viscous fluxes, respectively, can be written as

The mean flow is solved using a first-order backward-difference scheme, iterating in time until a steady- state solution is reached. For the purposes of this study, PSub dynamics are localized to the fluid-surface interface and are quantified using the amplitude of the complex transfer function and the phase difference between the disturbance pressure and the surface velocity response. In the code PSubs are modeled in two flavors: a mechanical PSub with realizable phases of +/-90 degrees, and a time-delay PSub with realizable phases of 0 to 360 degrees. The latter is inspired by the multi-input multi-output PSub described in [6]. The magnitude of the transfer function is calculated as

An LST-informed optimization procedure was developed to determine the optimal PSub frequency response and spatial design [7]. A single-frequency single-PSub study was performed in which the optimal PSub design was determined on a per-frequency basis. The reduction in disturbance wall pressure amplitude across the spectrum is shown in Figure 1.

From this analysis, a PSub was designed to target a broadband disturbance in the form of a finite-width pulse with a central frequency of 100 kHz. The FFT of the time-varying wall pressure data is shown in Figure 2. The pulse amplitude is shown at three instances in time in Figure 3.

3. Conclusions In this abstract we have successfully demonstrated the application of subsurface phononic materials towards reducing hypersonic boundary layer instabilities (both narrowband and broadband) and promoting transition delay. The phononic materials-disturbance interaction is well-examined and optimized for first and second mode disturbances. In the final paper we will present an analysis of the underlying physical mechanisms behind the stabilization, as well as outline the optimization procedure and provide a full characterization of the design space.

References [1] Leslie M Mack. Boundary-layer stability theory. Jet Propulsion Laboratory, 1969. [2] Alexander V. Fedorov, Norman D. Malmuth, Adam Rasheed, and Hans G. Hornung. Stabilization of hypersonic boundary layers by porous coatings. AIAA Journal, 39(4):605–610, 2001. [3] R. Zhao, T. Liu, C. Y. Wen, J. Zhu, and L. Cheng. Theoretical modeling and optimization of porous coating for hypersonic laminar flow control. AIAA Journal, 56(8):2942–2946, 2018. [4] Zhiyuan Jia, Yuhao Bao, Yangjun Luo, Dazhi Wang, Xiaopeng Zhang, and Zhan Kang. Maxi- mizing acoustic band gap in phononic crystals via topology optimization. International Journal of Mechanical Sciences, 270:109107, 02 2024. [5] Mahmoud Hussein, Sedat Biringen, Osama Bilal, and A. Kucala. Flow stabilization by subsurface phonons. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 471:20140928, 05 2015. [6] Carson L. Willey, Caleb J. Barnes, Vincent W. Chen, Kevin Rosenberg, Albert Medina, and Abigail T. Juhl. Multi-input multi-output phononic subsurfaces for passive boundary layer transition delay. Journal of Fluids and Structures, 121:103936, 2023. [7] Connor W. Klauss, Joshua Batstone, Mahmoud I. Hussien, and Christoph Brehm. Optimization of phononic subsurfaces for hypersonic boundary layer disturbance reduction. In AIAA SCITECH 2026 Forum, 01 2026.

ABSTRACT. Aerodynamic drag reduction on helicopter fuselages remains challenging due to massive three-dimensional flow separation. Although two-dimensional studies demonstrate that exploiting lock-in phenomena can suppress separation, extending these concepts to complex three-dimensional configurations is nontrivial. Previous studies have demonstrated that active flow control is effective for reducing the drag at the back of the ROBIN-mod 7. Motivated by earlier research, this study conducts a detailed parameter study to design a robust feedback controller for pulsed jet operation. The validity of utilizing the time coefficients of Proper orthogonal decomposition (POD) modes—which capture the dominant dynamic features of complex flow fields—as control variables has been extensively demonstrated. To implement this methodology, this study analyzed the relationship between the drag of the fuselage and POD modes. For the feedback control, POD mode’s time coefficient is estimated in real-time using a sparse array of surface sensors. A PID controller is integrated to dynamically adjust actuation parameters based on the predicted time coefficient. The control framework achieves drag reduction through the regulation of separation induced by the fuselage ramp.

ABSTRACT. The pressure distribution significantly influences both the compression performance and low-energy flow displacement capability of the bump. While existing bump design methods already enable inverse design of bumps based on a given three-dimensional pressure distribution, the selection of appropriate pressure distributions as design inputs remains an understudied area. This paper proposes a novel bump design method based on the concept of zoned pressure control. This approach divides the bump surface into distinct regions, where inviscid pressure distributions are specified for each zone. These pressure distributions are then used as inputs for the inverse design of the bump. The resulting bump, under viscous real-world conditions, forms a separation zone near the leading edge with high lateral pressure gradients and a distinct isentropic compression zone in the mid-to-rear section. Consequently, a novel bump is achieved that combines effective low-energy flow displacement with excellent total pressure recovery. Using this method, a novel bump was designed for Mach 2 conditions. First, numerical simulations were employed to analyze and validate the performance of this bump. Subsequently, the variations in aerodynamic performance with angle of attack, ranging from -4° to 4°, and the underlying reasons for these changes were investigated. Finally, the aerodynamic characteristics of the bump were analyzed and studied across a sideslip angle range of 0° to 17°. The results indicate that the bump designed using this novel method achieves a total pressure recovery of 0.936 and a low-energy flow displacement ratio of 86% on the symmetry plane, demonstrating excellent overall performance under design conditions (Mach 2). Moreover, the bump functions effectively across a wide range of operating conditions, maintaining stable performance within angle of attack ranges of -4° to 4° and sideslip angle ranges of 0° to 17°.

ABSTRACT. 1. Introduction The morphology of explosion induced craters in underwater near-bottom explosions significantly affects underwater structures, the intensity of sediment resuspension, and the extent of environmental disturbance.This phenomenon holds substantial engineering relevance for research in marine blasting, submarine pipeline protection, sediment disturbance, and marine ecosystem conservation.[1, 2, 3, 4, 5, 6] Numerical simulations, in contrast, can capture transient processes hard to observe in experiments; they avoid the over idealization of theoretical models; and they offer controllable cost and high repeatability. The numerical simulation of near bottom underwater explosions involves strong shock wave discontinuities, large plastic sediment deformation, and strong fluid-solid coupling. The Lagrangian method can precisely delineate and track multi-material interfaces. Bubble pulsation and the resulting large plastic deformation in sediment during near bottom underwater explosions often lead to severe mesh distortion and numerical divergence. In contrast, the Eulerian method employs a fixed grid, thereby avoiding the mesh distortion drawback inherent in the Lagrangian method. Appropriate interface algorithms are needed for multi material coupling, but they often introduce non conservative errors at material boundaries. Furthermore, the significant differences in properties between detonation products and sediment at the interface often cause severe numerical oscillations. Additionally, the Eulerian method struggles to accurately describe the deformation, failure, and ejection processes of plastic underwater sediment. We proposes a coupled SPH-FVM numerical method for high-fidelity simulation of the complex multi-phase coupled processes in underwater explosions involving saturated clay sediments. The explosives and water medium are discretized using the Finite Volume Method (FVM) grid, which accurately captures fluid dynamic behaviors such as shock wave propagation and bubble pulsation. SPH particles represent saturated clay to model large deformation and ejection, while the immersed boundary method handles fluid solid coupling and ensures interface conservation. 2. Methodology 2.1 Coupled method The schematic of the solution process for the coupled SPH-FVM method is shown in figure1. The immersed boundary method (IBM) is employed to facilitate information exchange between SPH particles and FVM grids. In the IBM, grid cells are first classified based on SPH particle distribution relative to the interface. Type-I cells are FVM grids that contain SPH particles on both sides of the interface; here the interface is treated as a virtual solid boundary. Type-II cells are FVM grids with SPH particles on only one side or on neither side. For Type-II cells, the conventional FVM approach is used to compute the grid fluxes. Type-I cells treat their interfaces as solid-wall boundaries and apply reflective conditions to calculate interfacial fluxes. Coupling SPH particles to FVM grids requires velocity continuity at the interface nodes. The velocity at the interface grid nodes is then obtained based on the momentum conservation condition. The allocation of physical quantities from FVM grid centers to SPH particles is accomplished through shape function interpolation. 2.2 Experiments Experiments were performed under conditions identical to the numerical simulations to verify the accuracy and reliability of the model results. The experimental platform uses an electro explosive device with high repeatability as the detonation source. A specially constructed 2.0×2.0×2.0 m high-strength, high-transparency explosion-proof glass water tank was used. Its excellent optical clarity ensures that high-speed photography can capture high-definition records. This platform provides significant convenience for crater morphology measurement. The crater is fully exposed for measurement by simply pumping the water down to the top of the steel box. A ruler can then be used directly for precise and efficient measurement of crater dimensions, greatly simplifying the post-processing workflow. 3. Conclusions We establishe numerical models for 0.36 g TNT detonated at standoff heights of 0 cm, 3 cm, and 6 cm above the surface of saturated clay sediment. Corresponding underwater near-bottom explosion experiments were conducted under the same conditions as the numerical simulations to validate the accuracy and reliability of the numerical models. Figure 2 presents a comparison of numerical and experimental outcomes for a 6cm standoff height. The results demonstrate that the numerical simulations can fully reproduce the evolution processes of bubbles and craters during near-bottom underwater explosions. They reveal how crater formation and bubble expansion in clay vary with standoff height, and the mechanisms behind these differences. The findings of this research provide strong support for advancements in underwater blasting, underwater pipeline protection, and sediment disturbance studies.

ABSTRACT. 1. Introduction As a hydrodynamic kinetic energy weapon that achieves high-speed, low-drag motion based on cavitation effects, supercavitating projectiles can reach velocities on the order of kilometers per second. When addressing underwater threats, supercavitating projectiles launched across media typically need to enter the water at a small angle and high speed. During this process, disturbances from the muzzle flow field and the external wave environment are unavoidable, subjecting the projectile to asymmetric impact loads. This exacerbates the asymmetry of the cavity shape and may even induce ricochet, affecting strike accuracy and reliability[1-2]. In contrast, underwater launching offers a new approach for the application of supercavitating projectiles. By launching supercavitating projectiles via underwater artillery, the impact of crossing the media interface can be effectively avoided, enabling rapid strikes against underwater targets while circumventing aerial threats to the ammunition. However, the underwater launch environment imposes specific requirements on projectile design: to maintain motion stability, the projectile needs a large length-to-diameter ratio; simultaneously, to ensure damage efficiency, a projectile structure with high specific kinetic energy (i.e., kinetic energy per unit cross-sectional area) must be adopted. Based on ammunition design theory, the discarding sabot structure is the only feasible solution to meet the above performance requirements and has thus become a critical bottleneck in the development of underwater-launched supercavitating projectile technology. APFSDS primarily consists of penetrator and sabot. The sabot separates after the projectile exits the muzzle, influenced by the combined action of the in-bore propellant gas and the fluid resistance ahead of the projectile. Research indicates that the hydrodynamic effects and mechanical interactions between the sabot and the penetrator during the separation process directly influence the smoothness of sabot separation, subsequently affecting the flight stability, accuracy, and armor-piercing performance of the penetrator. Sangeeta et al.[3] quantified the effects of various factors during the in-bore motion and external ballistic phases on the dispersion of 120 mm armor-piercing projectiles. They found that a reduction in residual spin rate by approximately 5% led to a 56% decrease in the first maximum yaw, indicating a significant correlation between the two. The effectiveness of their proposed design was validated through simulations and flight trials. Acharya et al. [4-5] divided the sabot separation process into four stages based on its characteristics and defined three transition points for theoretical analysis and calculation of each stage. They discussed the influences of mechanical collision, shock waves, aerodynamic forces, propellant gas, and gravity on sabot separation. The results indicated that mechanical collision and aerodynamic forces are the major sources of disturbance to the penetrator. Zhang et al.[6] investigated the influence of spin rate on the sabot discard characteristics of armor-piercing fin-stabilized discarding sabot (APFSDS) projectiles launched from a rifled barrel. Their study revealed that an increase in spin rate exacerbates the asymmetry of sabot discard, thereby altering the relative position between the sabot and the penetrator as well as the surface pressure distribution on the penetrator, which subsequently affects its aerodynamic characteristics. They determined that 1000 rad/s is the optimal spin rate for this type of rifled barrel gun. Huang et al.[7] explored the asymmetric discard behavior of sabots under different angles of attack. They discovered that a larger angle of attack intensifies the asymmetry of sabot discard, resulting in a more complex pressure distribution on the penetrator surface and more drastic changes in aerodynamic forces. Furthermore, this aerodynamic interference has a more pronounced impact on firing accuracy and flight stability as the angle of attack increases. Based on the aforementioned research background, and to advance the theoretical framework for underwater launched supercavitating projectiles, this paper establishes a numerical simulation method for the underwater launch process of supercavitating projectiles using the Schnerr-Sauer cavitation model and multiphase flow models. It systematically investigates the sealed underwater launch process for three different sabot structures at a muzzle velocity of v0 = 605.4 m/s. The analysis focuses on cavity evolution characteristics, flow field distribution patterns, projectile motion characteristics, and sabot separation behavior, aiming to provide a reference for trajectory prediction and structural optimization of underwater-launched supercavitating projectiles. 2. Methodology This study employs the VOF model for the numerical simulation of fluid motion, treating fluid elements as a single medium with variable density and viscosity, which effectively captures the interface between immiscible media. During the underwater launch process, the propellant gas inside the barrel is subjected to extreme high-temperature and high-pressure conditions. This study considers the compressibility effects of the gas by introducing the ideal gas equation of state to accurately describe the relationships between gas pressure, temperature, and density. For turbulence numerical simulation, the SST k-ω turbulence model is adopted for closure. Furthermore, the Schnerr-Sauer cavitation model is employed to simulate cavitation effects, treating the fluid medium as a continuous homogeneous mixture, simplifying the handling of phase interfaces, and accurately calculating the water vapor volume fraction through the gas-liquid mass transfer rate. To investigate the sabot separation characteristics and the evolution of the surrounding flow field of APFSDS projectiles launched underwater, this paper selects three typical sabot configurations as research objects: the conical sabot, the flat-nosed sabot, and the chamfered leading-edge sabot(Figure 1). The underwater sealed launching method is adopted. At the initial moment, the gun barrel is filled with high-pressure gas, while the exterior is a stationary water domain. During the projectile motion, the CFD equations and the classical interior ballistic equations are solved jointly. These are used to calculate the propulsive force from the propellant gas behind the projectile and the resistance force ahead of the projectile, respectively. At each time step, the fluid resistance ahead of the projectile and the propulsive force behind it are updated synchronously and substituted into the multi-body (sabot and penetrator) DFBI equations to calculate the projectile's trajectory. At the moment just before the projectile exits the muzzle, the complete in-bore flow field distribution—including pressure, velocity, and temperature—is obtained by solving the interior ballistic equations. This solution is then coupled with the exterior flow field computed via CFD to establish the global flow field state at the instant of muzzle exit. Thereafter, the subsequent evolution of the flow field is entirely governed by the CFD simulation, ultimately providing comprehensive data on the sabot separation process and the associated flow field characteristics during underwater launch.

3. Conclusions This study investigates the sabot configuration of APFSDS projectiles during underwater launch. Three sabot geometries (conical (s2), flat-nosed (s3), and chamfered leading-edge (s1) ) are selected as the research subjects. Their hydrodynamic characteristics, flow field distributions, and sabot separation behaviors are systematically analyzed. The main conclusions are drawn as follows: (1) The underwater launch process of the three types of APFSDS projectiles (s1, s2, s3) exhibits a separation mechanism jointly dominated by mechanical action and fluid action, which differs significantly from the separation behavior of discarding sabot projectiles in traditional air media. The collision between the sabot and the water body marks the initiation of the separation process. At the moment of muzzle exit, the separation behavior has entered the initial development stage. The front end of the sabot and the projectile's cavitator impact the water body at high speed, forming a high-pressure region while simultaneously generating a cavity that envelops the discarding sabot. Concurrently, the rear end of the sabot continues to be propelled by the propellant gas. At the mechanical action separation point, the sabot's motion constraint on the penetrator is released. Notably, at this moment, the tail of the projectile contacts the high-pressure region at the front of the sabot, causing a slight increase in the projectile's velocity (Figure 5). Upon entering the fluid-dominated stage, the cavities and pressure fields of the projectile and sabot gradually separate. The influence of the discarding sabot's flow field on the projectile's motion progressively diminishes until the forces acting on the projectile become completely stable (Figure 4), marking the transition to the external ballistic flight phase.

(2) The conical sabot (s2) exhibits pronounced low-drag characteristics. Due to its streamlined front end, both the mechanical separation point and the hydrodynamic separation point are noticeably delayed, postponing the separation process. However, precisely because of the minimal load borne by its front end, the overall stress on the sabot is uniform and stable, rendering it less susceptible to local buckling or severe deformation. The flat-nosed sabot (s3) achieves an earlier separation response at the cost of an increased frontal area exposed to the flow. Its hydrodynamic separation point is significantly advanced, resulting in a rapid discarding action, which is advantageous for shortening the sabot's disturbance stroke. Nevertheless, the high hydrodynamic pressure acts directly on the frontal plane, readily inducing edge curling and plastic deformation, thereby imposing stricter requirements on the sabot material and structural strength. The chamfered leading-edge sabot (s1) strikes a balance between the two. The chamfered structure mitigates the peak resistance at the front end to a certain degree, causing its separation point to be slightly later than that of the flat-nosed sabot but still significantly earlier than that of the conical sabot. Simultaneously, the chamfered edges effectively alleviate local stress concentration caused by flow separation. Furthermore, its velocity attenuation is the smallest among the three sabot types, making it an effective compromise solution that balances separation rapidity with structural reliability.

Acknowledgements The present work are financially supported by the Fundamental Research Funds for the Central Universities,（No.30925010417，No.30925020220）and National Natural Science Foundation of China Youth Fund（52301374）.

References [1] Zeqing G, Xuepu Y, Shuai S, et al. Experimental study on the high-speed water entry of cylinders at shallow angles[ J]. Physics of Fluids, 2024, 36(10): 104123. [2] Xuan Z ,Yonggang Y ,Xinwei Z .Numerical simulation and analysis of the 3D transient muzzle flow field of underwater artillery[J].Ocean Engineering,2023,284 [3] Panda S S ,Gite L ,Anandaraj A , et al.Dispersion sensitivity analysis & consistency improvement of APFSDS[J].Defence Technology,2017,13(4):316-322. [4] Acharya R, Naik S. Perturbation of Initial Stability of an FSAPDS Projectile[J]. Defence Science Journal, 2006, 56(5): 753-768 [5] Acharya R ,Naik S. Motion analysis during sabot opening process[J]. Defence Science Journal, 2007, 57(2): 229-241 [6] Chengqing Z ,Huiyuan W ,Ting L , et al.Sabot Discard Characteristics under Different Spin Rates of the Rifled Barrel Launching APFSDS[J].Shock and Vibration,2021. [7] Huang Z ,Xia C ,Cao Y , et al.Numerical investigations on the sabots discard process of an APFSDS at different angles of attack[J].The Journal of Engineering,2019,2019(13):373-378.

ABSTRACT. Offshore wind farms are growing fast, but their underwater noise impact is still hard to predict accurately. In this work, we focus on HORSES3D, an open-source high-order solver developed to simulate wind-turbine aeroacoustics and sound transmission into water with higher fidelity than reduced-order models. HORSES3D combines discontinuous Galerkin methods, actuator-line turbine modeling, immersed boundaries, and a diffuse-interface air–water formulation to capture key flow and acoustic mechanisms in a unified framework. It is parallelized with MPI and accelerated with OpenACC for large-scale simulations of offshore turbines and wind farms. Initial studies with a 22 MW reference turbine indicate that aerodynamic noise can contribute to underwater sound levels relevant to marine species, especially when many turbines operate together. These results highlight HORSES3D as a practical tool to quantify offshore wind acoustic footprints and support quieter, environmentally responsible wind farm design.

ABSTRACT. 1. Introduction Recent advances in high-fidelity modelling techniques and computational resources have positioned Computational Fluid Dynamics (CFD) as a transformative tool in marine hydrodynamics, enabling ship performance assessment throughout the design lifecycle. Compared with conventional towing-tank experiments and full-scale sea trials, CFD enables cost-effective full-scale simulations that reduce Reynolds number–related scale effects inherent in model testing [1]. However, despite its growing adoption, the reliability of CFD predictions depends on rigorous verification and validation (V&V), as modelling assumptions and numerical choices remain significant sources of uncertainty. Consistent with ITTC recommendations, robust V&V is therefore essential for the dependable application of CFD in ship hydrodynamic performance prediction. This study introduces a comprehensive, step-by-step CFD framework for predicting the calm-water hydrodynamic performance of a bulk carrier using the SINTEF Ocean Bulk Carrier (SOBC-1) as a reference hull. High-fidelity simulations reproduce the model test results reported by SINTEF Ocean AS [2], including resistance, propeller open-water performance, and self-propulsion characteristics. While many benchmark studies focus primarily on resistance prediction, the present work integrates resistance, propulsion, and wake analyses within a unified full-scale simulation strategy, enabling a more complete evaluation of the propulsion performance chain. The ship scale study aims to eliminate the Reynolds scaling effect. Furthermore, the study seeks to identify the most critical step(s) within the evaluation process that may contribute to significant prediction errors. Numerical predictions of propeller open-water performance were compared with model test data. Resistance and propulsion simulations were carried out at the design waterline (DWL) under free sinkage and trim. The self-propulsion point (SPP) was determined using both a rotating propeller (sliding-mesh) and a body force (virtual disk) method. Their results were compared in terms of accuracy, computational cost, and mesh requirements, with recommendations for the more efficient approach. A three-dimensional wake analysis was also conducted at the same conditions. All simulations, except the propeller open-water analysis, were conducted at full scale. Grid-sensitivity studies confirmed the reliability of the numerical approach. CFD results agreed well with model tests and the predicted outcomes were found to aligned well within the uncertainty limits. 2. Methodology Numerical Approach: In this study, the commercial CFD code STAR-CCM+ (version 2506) was employed. An unsteady RANS solver was used in conjunction with the SST k-ω turbulence model. For both resistance and propulsion simulations, the free surface effects were resolved using the VOF method with a modified High-Resolution Interface Capturing (MHRIC) scheme. Calm-water resistance and self-propulsion simulations were performed using the Dynamic Fluid–Body Interaction (DFBI) module, allowing free heave and pitch motions. The simulations were conducted in accordance with ITTC-recommended procedures for propeller open-water, resistance, and self-propulsion analyses. Computational Domain: Propeller open-water performance prediction analysis was performed in model scale (1:32.0). The computational domain, including the refinement zones and mesh distribution used for the propeller open-water performance prediction, is shown in Figure 1.

Figure 1: Computational domain and mesh distribution for propeller open-water simulations. Calm-water resistance was predicted at full scale using a half-body model to minimize computational cost, while full-body simulations were used for the SPP analysis with a rotating propeller. The computational setup, including domain, refinement zones, boundary conditions, and mesh distribution, are shown in Figure 2. A mesh sensitivity analysis was performed for the propeller open-water and R&P analyses, following ITTC (2008) verification guidelines. A uniform grid refinement ratio of √2 was adopted between successive mesh levels. The resulting convergence ratio satisfied the monotonic convergence criterion (0<Ri<1), indicating consistent grid convergence.

Figure 2: Computational domain and mesh distribution for resistance and SPP analysis. 3. Results and Discussions Propeller Open-water Performance Analysis: Open-water simulations were performed under the same conditions as the physical model tests, with a propeller speed of 14 rps and water temperature of 16 °C. A rigid-body motion (sliding-mesh) approach was applied to model propeller rotation, and a wall-function method with an average y⁺ ≈ 30 was applied near blade surfaces. CFD-derived thrust and torque coefficients showed good agreement with EFD data as shown in Figure 3. However, larger discrepancies appeared at higher advance coefficients (0.7–0.9), where CFD underpredicted thrust and overpredicted torque, reducing efficiency (η₀). Regarding the discrepancies observed in CFD predictions for a model scale propeller, several researchers have noted that standard turbulence models, such as k–ω SST model, exhibit limited accuracy at higher advance coefficients. Transition-sensitive models, such as the γ–Reθ transition model, can capture the laminar-to-turbulent transition on propeller blades in model-scale tests [3]. Thus, better agreement could potentially be achieved by employing the transition model; however, this approach could not be implemented within the scope of the present study. Propeller wake vortices were visualized using the Q-criterion (Figure 3).

Figure 3: Propeller open-water performance curve (left); Velocity distribution Q-Criterion of the model propeller (right). Ship Resistance Analysis: Calm-water resistance simulations were conducted at full scale using the SOBC-1 hull geometry, with hull roughness set to the default value of 150 μm. Vessel speeds ranged from 8 to 17 knots at the design draft (11.0 m forward and aft) under even keel conditions. CFD and EFD results show reasonable agreement across all speeds, with better accuracy at higher speeds as shown in Table 1. Larger discrepancies at speeds below 12 knots suggest that finer mesh resolution near the hull and free surface is needed to capture small-amplitude wave-making resistance, as a uniform mesh was used for all cases. The percentage error ranges from approximately 1.5% to 10%. Table 1: Resistance, trim and sinkage comparison between EFD and CFD results. Vs RES_EFD RES_CFD Trim_EFD Trim_CFD Sinkage AP_EFD Sinkage AP_CFD Sinkage FP_EFD Sinkage FP_CFD [knots] [kN] [kN] [Deg] [Deg] [m] [m] [m] [m] 8 126.24 140.102 -0.036 -0.044 -0.011 0.013 -0.131 -0.131 10 203.77 217.2 -0.063 -0.068 -0.006 0.017 -0.215 -0.209 12.5 345.79 357.2 -0.106 -0.11 0.004 0.03 -0.347 -0.334 15 557.69 566.4 -0.167 -0.168 0.025 0.047 -0.529 -0.511 17 871.96 885.2 -0.228 -0.228 0.049 0.072 -0.709 -0.7

Self-propulsion Point Prediction Analysis: The SPP analyses were performed at 15 knots under DWL loading, allowing free heave and pitch with zero initial trim. Full-body CFD simulations used two propulsion modeling methods: a rotating propeller (sliding-mesh) and a body-force propeller. Since CFD underpredicted thrust, EFD-based KT–KQ curves were applied in the virtual disk method. Both methods showed good agreement with model test results, with about 4% difference in propeller speed and 7% in power. Trim and sinkage differences were negligible. The rotating propeller method required much higher computational cost, while the body-force method offered a practical balance between accuracy and efficiency, though it cannot capture unsteady flow structures or vortex shedding. Therefore, when detailed prediction of unsteady propeller flow characteristics is not required, the body-force propeller method is recommended. Figure 4 (left) presents the free-surface wave profile and flow field for the self-propulsion case with an actual rotating propeller, where the propeller-induced vorticity and wake structures are well resolved. Conversely, when the body-force propeller method is applied, the detailed vortex shedding phenomena are not reproduced, as illustrated in Figure 4 (right).

Figure 4: Flow pattern from SPP analysis with the rotating propeller (left), with a virtual disk (right). In the physical model tests, the 3D wake was measured at the propeller plane for both ballast and design waterline conditions at a ship speed of 15 knots. A comparison between the CFD and EFD wake profiles at DWL will be included in the final draft of the full paper. 4. Conclusions This study used CFD to evaluate the calm-water performance of the SOBC-1, including propeller open-water, resistance, and self-propulsion analyses. CFD results agreed well with model tests, though thrust was slightly underpredicted and torque overpredicted at higher advance coefficients (J = 0.7–0.9). Future work should apply the γ–Reθ transition model for improved accuracy. Self-propulsion simulations at 15 knots using both rotating and body-force propeller models matched experimental data. While the rotating propeller captured detailed unsteady flow at higher computational cost, the body-force method proved more efficient for routine powering predictions. The choice of approach should follow ITTC guidelines based on required fidelity and study objectives.

References [1] V. Krasilnikov, V. S. Skjefstad, K. Koushan and H. J. Rambech, "A Calibration Study with CFD Methodology for Self-Propulsion Simulations at Ship Scale," J. Mar. Sci. Eng., vol. 11, no. 7, pp. 1342, 2023, https://doi.org/10.3390/jmse11071342.

[2] H. J. Rambech, Report: Calm Water performance tests, SOBC-1 (SINTEF Ocean Bulk Carrier no.1), SINTEF Report, Project No. 302006033‐09, August 25, 2021, Enterprise Number: NO 937 357 370 MVA., Norway. [3] K.W. Kim, K. J. Pike, J. H. Lee, S. S. Song, M. Atlar and Y. K. Demirel, " A study on the efficient numerical analysis for the prediction of full-scale propeller performance using CFD," Ocean Engineering, 240, 109931, 2021.

ABSTRACT. In this study, Poisson-free incompressible flow solvers—including the Artificial Compressibility Method (ACM) with its variants (EDAC and BVACM), the KRLNS equations, and the Lattice Boltzmann Method (LBM)—are applied to three-dimensional homogeneous isotropic turbulence on the same GPU platform to compare their computational efficiency. All simulations are conducted at a fixed Reynolds number of Re = 4000 and Mach number of Ma = 0.04. For well-resolved cases, fourth-order ACM achieves higher computational speed than both the fractional step method (FSM) and LBM at the same L2-norm error. Its high-order spatial discretization reduces phase errors, and RK4 time integration allows larger steps. ACM outperforms FSM, partly due to avoiding the global Poisson solve, while LBM may suffer from increased memory traffic due to non-coalesced access patterns arising from the directional streaming process, as well as from a time step limited by the Mach number. For under-resolved simulations, BVACM is more efficient than LBM by damping high-frequency components via bulk viscosity, preserving large-scale structures while improving stability.

ABSTRACT. Graphics Processing Units (GPUs) are becoming widespread in High Performance Computing (HPC) clusters thanks to their ability to highly parallelize computational work. Unlike previous processing units, GPUs do not have a universal programming language that can be used to develop code that runs on different architectures. In this work, we ported an existing Computational Fluid Dynamics (CFD) solver to run on GPUs and able to achieve performance portable results through the Kokkos library. The numerics of the solver and the data structures used to store the physical quantities have been adapted to maximize portability and performance within Kokkos. After porting the code to GPUs, it was tested on the Polaris and Aurora supercomputers at Argonne National Laboratory, which use NVIDIA and Intel GPUs, respectively. The solver achieved nearly perfect weak scaling on both supercomputers and showed consistent performance across GPU architectures and numerical orders. Moreover, we showed that different architectures perform specific functions faster than others, depending on their internal design and optimization. To assess the soundness of the new software, it was validated on a Taylor–Green Vortex simulation against its original version, which had been previously validated against several solvers on the same case. Additional validation is provided by the agreement between the simulation and experimental results for the second-mode wave frequency over a flat plate under hypersonic conditions.

ABSTRACT. The increasing use of GPU-based architectures for large-scale computational fluid dynamics introduces numerical effects that are negligible on traditional CPU platforms but can significantly impact large-eddy simulation results. In particular, the non-associative nature of massively parallel floating-point reductions on GPUs can amplify round-off errors, leading to architecture-dependent behavior. In this work, we investigate the role of the mathematical formulation of the LES spatial filter width in determining numerical robustness and hardware portability. Focusing on the flow past a circular cylinder at Reynolds number 3900, we show that the classical least-squares filter-width formulation, although mathematically sound, exhibits strong sensitivity to floating-point accumulation errors when implemented in a face-centric GPU framework, resulting in non-physical wake behavior. To address this issue, we propose a reformulated filter-width definition that combines a metric-based estimate with a realizability constraint linked to the minimum grid spacing. The new formulation mitigates the amplification of numerical noise and restores consistency between CPU- and GPU-based simulations. Results demonstrate excellent agreement with reference experimental data and CPU computations, highlighting the need to consider numerical robustness and hardware portability as intrinsic requirements of LES modeling on modern heterogeneous HPC systems.

ABSTRACT. Accurate prediction of rotating flows is essential in engineering applications such as aircraft engines and wind turbines, yet it poses significant computational challenges. This study presents a high-performance computational framework for rotational flow simulation by integrating the Lattice Boltzmann Method (LBM) with a level-set moving boundary technique, the Wall-Adapting Local Eddy-viscosity (WALE) Large Eddy Simulation (LES) model, and a Hybrid Recursive-Regularized (HRR) collision model. To enhance computational efficiency on desktop-level CPU/GPU platforms, several CUDA optimization strategies are introduced, including a fast integer-based search method for real-time node updates and a dynamic thread index alignment scheme for improved interpolation at grid boundaries. A multi-GPU parallel computing architecture is developed to enable large-scale simulations. The method is validated through Taylor–Couette flow simulations, showing agreement with analytical solutions. Large-scale numerical simulations of low-altitude rotor flows are conducted for three configurations: a single rotor, a rotor with a stator, and a counter-rotating open rotor. Results demonstrate that the counter-rotating configuration significantly reduces rotational kinetic energy in the wake region while improving axial velocity, highlighting the effectiveness of the proposed framework in capturing complex flow features with high accuracy and computational efficiency.

ABSTRACT. Moment closures based on maximum-entropy principals offer efficient and effective strategies for predicting high-speed and even hypersonic non-equilibrium gaseous flows with shocks ranging from the near-equilibrium continuum regime well into the so-call transition regime, without necessitating the solution of the underlying high-dimensional Boltzmann kinetic equation. A novel non-oscillatory high-order flux reconstruction (FR) scheme is considered for the solution of both the second- and fourth-order maximum-entropy closures associated with a representative one-dimensional kinetic equation having a relaxation-time approximation for the collision operator. In the proposed approach, a solution smoothness detection strategy is combined with standard limiting technique applied to the linear modal content of the high-order FR schemes so as to effectively treat shocks and discontinuities while avoiding unwanted oscillations and still retaining high-order solution accuracy in regions of smooth and resolved solution content. The potential effectiveness of the proposed non-oscillatory high-order FR schemes are demonstrated via application to a range of problems, including those with shocks. The relationship of the proposed non-oscillatory method to other similar strategies for FR and discontinuous Galerkin methods and future application of the approach to the multi-dimensional case are also both discussed.

ABSTRACT. Entropy-stable split-form high-order discretizations have emerged as robust tools for simulating nonlinear hyperbolic conservation laws, yet general convergence results remain limited, even for smooth solutions. This work establishes convergence for entropy-stable split-form continuous summation-by-parts semi-discretizations applied to scalar conservation laws and symmetric hyperbolic systems with homogeneous flux functions on periodic domains, under smoothness assumptions on the exact solution and fluxes. The analysis proceeds directly from the semi-discrete error equation, avoiding projection-based error decompositions and convergence arguments that rely on the linearized (first-variation) stability of the scheme. We derive a nonlinear error evolution inequality in an energy norm and show that an error bound can be found from the solution of a constant-coefficient Riccati initial-value problem. For sufficiently small mesh spacing, this error bound exists on any fixed finite time interval and tends to zero as the mesh is refined, implying convergence despite the possible presence of local linear instabilities. The inviscid Burgers equation is used as a canonical example to illustrate the theory and to clarify the roles of consistency, entropy stability, and nonlinear error growth in determining convergence behavior.

ABSTRACT. Multidimensional summation-by-parts (SBP) operators can be used to construct high-order, conservative, and stable discretizations on complex geometries. While tensor-product SBP operators are widely used, the adoption of multidimensional SBP operators has been slow, in part because they are difficult and tedious to implement. To help address this difficulty, this work presents an explicit formula for first-derivative SBP operators. This formula produces discretizations whose spectrum and accuracy are effectively independent of the quadrature rules used to define the solution nodes. Numerical experiments are provided that verify this behavior.

ABSTRACT. Accurate, robust and cost-efficient algorithms are an ongoing research topic in many fields. For CFD, high-order methods are especially appealing, as they offer the potential for higher accuracy at a reduced cost. Discontinuous Galerkin-type schemes are ideal for large, parallel computations that are needed for CFD, but are still considered high-cost. Cicchino et al. recently proposed Non-linearly Stable Flux Reconstruction (NSFR) schemes within the family of entropy stable schemes, which are provably stable for non-linear problems. This means the scheme is stable without any dependence on mesh resolution, paving the way for mesh adaptation. This could provide significant cost benefits for large eddy simulation (LES), which requires fine grids for accurate simulations. As the solution advances in time, the mesh adaptation algorithm can adapt the mesh as needed, with guaranteed stability from the scheme. In this work, the non-linear stability guarantee of NSFR will be used to perform Implicit LES (ILES), with a focus on mesh adaptation for the LES of turbulent flows.

ABSTRACT. Recently, our group proposed the kinetic energy and entropy preserving flux reconstruction (KEEP-FR) scheme, which is a non-dissipative and robust convective scheme. The achievement of the robust non-dissipative computations for convective terms has enabled the evaluation of diffusion schemes without the interference of numerical dissipation. Based on the robust and accurate convective schemes, the diffusion scheme, whose impact was previously overshadowed by the numerical dissipation, become important in the overall accuracy and robustness. In this study, we propose a new diffusion scheme in high-order FR framework. Existing diffusion schemes, such as the compact discontinuous Galerkin (CDG) and BR2 schemes, evaluate the second-order derivatives by correcting the polynomial distribution of primary variables and the flux. To ensure robustness, the CDG scheme determines common values at cell interfaces using values from only one side, which introduces nonphysical asymmetry into the scheme. Furthermore, we demonstrate that the BR2 scheme suffers from nonphysical numerical oscillations when convective numerical dissipation is absent. To overcome these issues, we introduce the single and symmetric reconstruction (SSR) diffusion-flux scheme, which evaluates second-order derivatives through a single and symmetric reconstruction process. The proposed scheme achieves numerical robustness without introducing nonphysical asymmetries, while simultaneously suppressing the numerical oscillations observed in the existing schemes. Furthermore, in this study, we employ combined mode analysis to analytically clarify the diffusion property of both the proposed and existing schemes. The numerical experiments are conducted, including the double shear layer and the viscous Taylor-Green vortex, to validate the capability of the proposed method.

ABSTRACT. 1. Introduction Non-oscillatory high-order schemes address the central challenge of accurately tracking dynamically evolving interfaces in compressible multiphase flows by mitigating numerical diffusion. To address the computational cost and robustness challenges faced by conventional high-order methods like weighted essentially non-oscillatory (WENO) and discontinuous Galerkin (DG), we have developed the volume integrated average and point value based multi-moment (VPM) family of schemes to better balance accuracy, efficiency, and numerical stability. The fundamental idea of the VPM scheme is to introduce point values (PVs) at cell vertices as independent DOFs, in addition to the volume integrated averages (VIAs) that serve as the sole computational variable in conventional finite volume (FV) methods. Being the conservative quantity, the VIAs are updated via a finite volume scheme, which computes numerical fluxes across cell boundaries. In contrast, the PVs are efficiently computed by solving the governing equations in their differential form, which approximates the pointwise derivative at the cell vertices. Having both VIAs and PVs available within each cell enables the construction of a higher-order polynomial on a compact stencil. However, prior VPM versions are confined to single phase flows, with core limitation lies in the PV update, which requires flux derivatives computed via a derivative-based Riemann solver. Specifically, the Roe solver which linearizes the flux derivatives as the product of the Jacobian matrix and the variable derivatives is employed. This approach is incompatible with more robust, multiphase-capable solvers (e.g., HLLC, AUSM) and is notoriously difficult to extend to multi-fluid systems.

This work overcomes these barriers by introducing the VPM-GF (VPM with Generalized Formulation) scheme. Its primary objectives are: 1) to reformulate the PV update to be compatible with a broad spectrum of standard Riemann solvers, 2) to extend the high-order VPM framework to compressible two-phase flows, and 3) to maintain geometric flexibility, accuracy, and robustness on arbitrary polyhedral unstructured grids. To achieve these, we incorporate the idea of the Generalized Finite Difference (GFD) method. Inspired by Li and Ren's study, we propose a novel VPM scheme featuring a three‑step PV update procedure, which will be presented in the methodology Section. Combined with limiter and positive-/bound-preserving conditions, the proposed scheme show excellent performance in both single- and two-phase flow simulations.

2. Methodology The inviscid compressible simulation in this study are governed by linear advection equation and Euler equations for single phase flows, and Allaire's five-equation model for two phase flows.

The VPM method employs both VIAs and PVs at cell vertices as prognostic variables. The reconstruction, flux evaluation and update procedure for VIAs adopt the methodology proposed in our prior VPM edition, while the PV treatment adopts a three-step procedure. First, a quadratic polynomial is constructed at each vertex using a least-squares method to achieve third-order accuracy. The stencil, hereinafter referred to as the "reconstruction stencil", comprises the directly adjacent vertices and neighboring cell centroids. Second, fluxes are computed at the midpoints linking the target vertex to each member of the reconstruction stencil. This set of midpoints defines the "flux evaluation stencil", which is obtained directly by geometric transformation of the reconstruction stencil, eliminating additional search and storage costs. Third, a least‑squares approximation on the flux evaluation stencil tields the flux derivatives at the target point, which are then used to update the PVs via the differential governing equations.

3. Numerical results A series of benchmark tests has been conducted to evaluate the performance of the VPM‑GF scheme in both compressible single- and multi-phase flow simulations. The sine wave advection and manufactured solution tests confirm the third-order accuracy of the VPM-GF scheme for solving advection and Euler equations on quadrilateral, triangular and polygonal meshes. The canonical Mach 3 forward step and double Mach reflection tests demonstrate the robust shock-capturing capability of the VPM‑GF scheme, which well resolves the shock waves and vortical structures with sharper resolution due to its low numerical dissipation. The performance of the scheme in compressible two-phase flow simulations is validated across one-, two-, and three-dimensional configurations. Two 1D tests are considered. The liquid column advection test presents sharp interfaces and pressure equilibrium over both short and long time scales, confirming the interface-capturing ability and thermo-dynamic consistency. And the water-air shock tube test shows well-captured expansion wave, shock wave and contact discontinuity, revealing the high accuracy, negligible numerical oscillation and strong robustness for problems involving extreme pressure and density ratio. In two dimensions, Complex flow structures with high resolution can be seen in the the canonical shock-R22 gas interaction and triple point problems, further verifying the interface-capturing ability of the VPM-GF scheme with low dissipation. Finally, the three‑dimensional shock-helium bubble interaction test illustrates the high resolution and temporal fidelity of our scheme, with the whole evolution of the flow structures accurately reproduced.

4. Conclusions The proposed VPM-GF scheme successfully extends the VPM method to the multiphase flow simulation through improvements on PV update. Instead of the derivative-based Roe Riemann solver, the improved PV treatment take advantages of the conceptual framework of GFD methods, offering several key advantages, including compatibility with diverse Riemann solvers, a unified Riemann-solver implementation for both VIAs and PVs, applicability to arbitrary polyhedral grids, and, most importantly, the capability to simulate compressible two-phase flows. A series of benchmark tests has illustrated the shock- and interface-capturing ability of the VPM-GF scheme with negligible oscillation and strong robustness, validating the potential of the VPM‑GF scheme for accurate and robust simulations of both compressible single‑ and multi-phase flows.

ABSTRACT. Numerical simulation of compressible multiphase flows based on the four-equation model is known to suffer from two fundamental difficulties with wave propagation and pressure equilibrium preservation. First, the mixture sound speed exhibits non-monotonic dependency with respect to the volume fraction, which leads to robustness issues in the resolution of shocks and acoustic wave propagation across two-phase regions. This difficulty can be mitigated by solving Allaire’s five-equation model augmented with infinitely fast phasic temperature equilibrium, from which solutions of the four-equation model can be recovered. However, when temperature is non-uniform, this augmented five-equation formulation still fails to preserve pressure equilibrium across material interfaces.

In this work, we propose a fully conservative numerical scheme that exactly preserves pressure equilibrium at the discrete level for the augmented five-equation model, for arbitrary initial distributions of temperature and volume fraction. Combined with the monotonic sound speed property of the five-equation formulation, the proposed pressure-equilibrium preserving scheme significantly improves robustness in the presence of strong multiphase interactions, including shock–interface interactions and advection of material interfaces.

ABSTRACT. A one-dimensional slug-capturing framework based on an always-hyperbolic five-equation two-fluid model is presented for gas–liquid flows in inclined and vertical pipes. Unlike conventional conditionally hyperbolic two-fluid models, the proposed formulation preserves hyperbolicity for all admissible flow states, eliminating the appearance of unphysical waves and ensuring that all waves in the domain are physically meaningful. To reduce excessive numerical dissipation in low-Mach-number regimes, a Roe–Turkel preconditioning technique is applied within a fractional-step finite-volume scheme. Although preconditioning typically imposes severe time-step restrictions, a von Neumann stability analysis of the complete hyperbolic–relaxation scheme demonstrates that the overall method retains the standard CFL stability condition for liquid–gas flows. The numerical framework is implemented on graphics processing units (GPUs) to fully exploit parallel hardware capabilities, achieving multiple orders-of-magnitude speedup compared with conventional CPU implementations. The model is validated against experimental data for air–silicone oil flows in inclined pipes, showing accurate prediction of mean slug frequency and void fraction across a wide range of inclination angles.

ABSTRACT. Introduction

Rocket plume impingement into a low ambient environment can entrain and accelerate surface grains to high velocities creating ejecta that may abrade pollute and damage neighboring assets with risk extending beyond the immediate landing zone. The risk is not only affected by whether grains are lifted but also by the distribution of particle exit conditions from the surface such as speed ejection angle and ballistic range in near vacuum where aerodynamic deceleration is reduced. Thus particle trajectories and velocity information are directly related to lander design including component standoff lengths shielding and approach limitations and mitigation strategies including plume deflection barriers.

Plume surface interaction driven dust transport has been predicted using a variety of modeling techniques and previous research indicates that assumptions and closure choices can significantly alter predicted particle velocities and angles. The systematic reliance of particle velocity on particle diameter and initial location in the plume footprint is further demonstrated by trajectory focused studies which motivates the use of quantitative hazard metrics like velocity distribution trajectory angle with respect to the surface and plane crossing particle fluxes. Eulerian Lagrangian approaches are popular because they extract these metrics directly by tracking particles and resolving compressible gas dynamics. However accuracy is highly dependent on drag modeling. It has been demonstrated that compressibility and rarefaction significantly affect projected particle transport.

To this end we present rhoCentralLPTFoam an in house density based compressible solution for particle laden impinging jets that are typical of plume surface interaction at low ambient pressure. With a one inch exit diameter and operating at a high nozzle pressure ratio and high nozzle temperature ratio with ambient pressure and temperature at one Pascal and three hundred Kelvin respectively the basic configuration uses a NASA Acoustic Reference Nozzle that discharges into a low ambient environment and impinges on a particle bath and a flat surface. This contribution establishes a reference dataset by reporting early time particle trajectories and ejection angle distributions particle speed distributions including high velocity tails and sensitivity to particle size and density.

Methodology

Numerical Framework

The continuous phase is solved using a density based compressible formulation derived from the rhoCentralFoam solver in OpenFOAM. The governing equations for mass momentum and total energy are discretized using a central upwind flux scheme providing robust shock capturing capability. Spatial reconstruction is second order accurate and time integration is explicit. Turbulence effects are modeled using the k omega SST model. The solver is fully compressible and suitable for strongly underexpanded jets at high nozzle pressure ratios. The work extends the baseline density based framework by incorporating a Lagrangian particle tracking module tailored for plume surface interaction conditions.

Eulerian Phase

The continuous phase air is governed by conservation of mass momentum and energy. The momentum equation includes pressure viscous stresses gravitational acceleration and drag interaction. In the one way coupled fluid particle interaction the drag feedback to the fluid is neglected. Heat exchange between fluid and particles is also neglected.

Lagrangian Phase

The particle phase is modeled using a multiphase particle in cell approach available in OpenFOAM. Computational parcels represent ensembles of physical particles and particle particle interactions are accounted for through a continuum stress model. This enables representation of locally dense particle regions that may arise during plume surface interaction.

The motion of particles follows Newtons second law and the time advancement is explicit and synchronized with the gas phase solver. In the present configuration fluid particle interaction is treated as one way coupling where the gas phase influences particle motion but particle feedback to the fluid is neglected. Extension to two way coupling is planned in future work to account for momentum and energy exchange for denser particle beds.

Forces Acting on Particles

Gravitational buoyancy pressure gradient and drag forces are considered. The drag coefficient depends on particle Reynolds number and relative velocity between gas and particle. Future developments will incorporate supersonic and rarefied aware drag correlations to account for high relative Mach number and low density effects typical of plume surface interaction in low ambient environments.

Results and Discussion

Flow and Geometry Configuration

The nozzle discharges into a low ambient environment and impinges on a particle bath. The configuration represents a strongly underexpanded jet.

Preliminary Results

Plume surface interaction creates a significantly time dependent dust lifting and transport field with a fast growing near surface high speed jet spreading downstream and the particles forming an expanding ejecta cloud.

Conclusions

A density based compressible Eulerian Lagrangian framework rhoCentralLPTFoam has been presented for plume surface interaction under low ambient conditions. The solver combines a central upwind scheme with SST turbulence modeling and a particle in cell formulation to resolve shock dominated impinging jets and dense particle regions. Preliminary results demonstrate strongly time dependent dust lifting rapid near surface jet spreading and formation of expanding ejecta clouds. Future work will incorporate two way coupling and supersonic rarefied aware drag correlations to improve physical fidelity in high Mach number and low density regimes characteristic of planetary landing scenarios.

ABSTRACT. Accurately and efficiently tracking the interface motion is of great importance in multiphase flow simulations. We propose a novel Front-Tracking method, the Edge-Based Interface Tracking (EBIT) method, to address this challenge. In the EBIT framework, marker points are exclusively placed on the grid edges, and their connectivity is implicitly described through a color vertex field. The interface is advected using a directional split scheme. To maintain the correspondence between markers and cell edges, we carry out a reconstruction step at the end of each advection stage. This step updates the marker positions on cell edges, transforming the Lagrangian markers into an Eulerian field. Such a transformation facilitates the parallelization of the EBIT method, particularly when the method is coupled with quad-/octree meshes, thereby improving scalability compared with traditional Front-Tracking methods.

In this work, we further introduce a simplified strategy to extend the EBIT method to three dimensions (3D). The directional split scheme used in the 2D version is directly generalized to 3D. Specifically, the 3D advection within a cubic cell along a given direction is decomposed into two 2D advection problems on the corresponding cube faces.

The 3D EBIT method on an octree grid is implemented within the free platform Basilisk. To verify the proposed algorithm, we conduct a series of kinematic and dynamic test cases and compare the results with those obtained using the PLIC-VOF method.

ABSTRACT. CFD simulations provide a powerful set of tools for gaining insights into the flow behavior of various fluids and gases. Their main drawback is the large computational effort required to simulate systems beyond laboratory scale.

To address this issue, data-assisted methods can be employed to improve simulation efficiency. The most prominent example is machine learning, where, for instance, deep operator networks have been applied to an expanding range of fluid dynamics problems. BEven though specialized methods can significantly reduce the size of the required training data sets, they often still demand a substantial amount of data. This, in turn, leads back to infeasible computational costs, particularly in cases where custom CFD models are required.

Recurrence CFD (rCFD) is an alternative data-assisted approach that aims to achieve a significant reduction in numerical cost while using as little data as possible. For systems that may exhibit turbulent flow but display recurrent flow patterns over longer time scales, a short CFD simulation can already contain sufficient information about the future system behavior. rCFD uses one or more such simulations to generate a database containing all encountered velocity fields. Assuming that this database covers all relevant system states, the method of analogs can be employed to extrapolate the system forward in time.

The rCFD variant employed in this work stores velocity fields as displacements of cell contents using propagator theory, achieving speed-ups of up to three orders of magnitude compared to conventional CFD simulations.

The method is demonstrated using a submerged entry nozzle (SEN) with an attached mold for the continuous casting of steel. In this setup, argon is injected together with the molten steel to promote favorable flow patterns. As a result, non-metallic inclusion particles (NMIs) present in the steel are more likely to be captured at the slag layer rather than transported downward into the finished product, where they would reduce its quality. The aim of this work is to develop a CFD model capable of simulating the capture of inclusion particles at the slag layer and to calibrate this model against experimental data using rCFD.

ABSTRACT. BYCSprayFoam is an efficient OpenFOAM-based Euler–Lagrange solver developed for gas–liquid two-phase rotating detonation engine (RDE) simulations. Built on OpenFOAM v1812, it extends our gas-phase rotating-detonation solver BYCFoam to incorporate two-way coupled Lagrangian spray parcels, with interphase mass, momentum, energy, and species exchange consistently embedded in the Eulerian multispecies conservation laws. Thermodynamics, transport, and finite-rate chemistry are evaluated through a Cantera coupling, enabling accurate NASA-polynomial thermochemistry and mixture-averaged transport for multi-species reacting flows. To address stiffness and localized time-step restrictions inherent to detonation–spray interactions, the solver employs operator splitting for cell-local chemistry integration and provides dual time stepping to relax physical-time constraints via pseudo-time iterations with consistent Eulerian–Lagrangian coupling. Multiple shock-capturing flux options are supported, and adaptive mesh refinem

ABSTRACT. Particle deposition and detachment in turbulent flows significantly affect industrial processes. In continuous steel casting, deposition of Non-Metallic Inclusions (NMIs) on the inner walls of Submerged Entry Nozzles (SENs) can lead to clogging and flow instability, while detached particles may re-enter the molten steel stream and become trapped in the final slab, reducing steel quality. Deposition is driven by turbulent transport and adhesion mechanisms such as interparticle contact, cavity bridges, and high-temperature sintering. Detachment occurs when hydrodynamic forces exceed the adhesion bonds that hold particles together or to the wall. In this study, a 4-way coupled CFD–DEM approach is used to investigate the formation, morphology, and weak points of near-wall clog structure under turbulent molten steel flow and provide a comprehensive view into both deposition and detachment processes.

ABSTRACT. In-flight aircraft icing is a complex multiphase phenomenon governed by the coupled interaction of airflow, supercooled water droplets, heat and mass transfer, and phase change on the surface. Although numerous numerical icing models have been developed, their predictive capability strongly depends on the physical assumptions used to describe droplet deposition and ice accretion processes. In particular, classical film-based formulations involve internal physical inconsistencies in the treatment of heat removal and phase change during droplet impact and freezing. In this work, a physics-based numerical modeling framework for air–droplet flow and ice accretion in aircraft icing problems is developed and physically justified using dedicated experimental observations. A series of laboratory experiments is performed to investigate the interaction of supercooled water droplets with an icing surface and the characteristic features of ice formation. Based on these observations, a refined physical interpretation of the ice accretion process is introduced into the numerical model. The air–droplet flow is described using the Reynolds-averaged Navier–Stokes equations for the carrier phase and the interpenetrating media approach for the dispersed droplet phase. Ice growth is simulated using a surface control-volume-based mass and energy balance formulation, with convective heat transfer obtained consistently from the flow solution. Numerical simulations are performed for the NACA 0012 airfoil under representative icing conditions. The results illustrate the capabilities of the proposed modeling framework and demonstrate a reasonable agreement between the numerically predicted ice shapes and available experimental data.

ABSTRACT. Numerical simulation provides a way for gaining insight of flow pattern with high-fidelity predictions but are computationally expensive. A hybrid machine learning model integrating a 3D convolutional neural network (CNN), long short-term memory (LSTM) units and multi-attention mechanisms was developed to predict the time-series evolution of two-dimensional flow patterns throughout the reactor domain. A omprehensive data analysis compared predicted and simulated values for solid volume fraction, horizontal and vertical solid-phase velocities, and pressure distribution in bubbling fluidized bed. Results demonstrated that the proposed prediction model preserves essential physical consistency with large reduction of computation cost. It offers a new approach for the enhancement of subsequent CFD simulations for time-series prediction.

ABSTRACT. Large-scale particle-laden flows are characterized by strong scale separation and meso-scale clustering, which render reliable macro-scale predictions with Computational Fluid Dynamics simulations difficult. In gravity-driven systems, particle interactions generate both correlated turbulent motion and collision-induced micro-scale fluctuations, described as granular temperature. Evidence of bidirectional energy transfer between these scales indicates that conventional dissipation concepts are insufficient.

This study analyzes the coupling between particulate-phase turbulent kinetic energy and granular temperature using Particle-Resolved Direct Numerical Simulations of dilute flows with Kolmogorov-scale particles. By applying spatial filtering and tracking individual particle histories, correlated and uncorrelated fluctuating energies are separated. The results show that correlated fluctuations dominate in clustered regions, while uncorrelated fluctuations prevail in dilute regions.

ABSTRACT. Cavitation flow is characterised by two-phase, compressible turbulent flows. Due to strong compressible effects, it requires a compressible flow solver to allow accurate capturing of pressure waves. RANS (Reynolds Averaged Navier-Stokes) simulations have demonstrated an inability to predict flow dynamic in the cavitation context. This is the reason why attention is directed towards ILES (Implicit Large Eddy Simulation) computations, in which numerical dissipation plays a role as a subgrid model. To our knowledge, most of the literature mainly focuses on 2nd-order accurate schemes. This work is consequently devoted to investigate the potential of a 4th-order accuracy scheme

In order to do that we simulate three-dimensional vortex cavitation on the wake of a circular cylinder. Results with second and fourth-order of accuracy will be compared and criticized on the role of the dissipation behavior regarding the cavitation prediction. Comparison and conclusion will be discussed at the conference.