Particle-resolved simulations are essential for studying dense-suspension rheology because they resolve the microstructure that controls macroscopic stress, but their computational cost remains substantial due to the tight two-way coupling between the fluid and the suspended particles. In the OpenFOAM-based fictitious-domain solver considered here, the dominant cost is the iterative determination of the forcing field that couples the two phases by enforcing rigid-body motion inside the particles. We address this bottleneck by learning a solver-internal quantity rather than replacing the solver itself: a hybrid graph neural network (GNN) and U-Net model predicts how the forcing field changes each time step and initializes the solver’s existing iteration. The governing equations, solver loop, and convergence criterion remain unchanged.On simulation intervals not used during training, the learned initializer cuts the mean forcing error – an accuracy measure – by about 95% relative to naive persistence (repeating the previous value). Embedded in the full solver, it yields an iteration-reduction factor of 2.54 and an overall speedup of 2.13× after accounting for inference cost. The time-averaged shear stress changes by only −1.65% relative to the unmodified solver. Trained on a single simulation, the model retains iteration speedups between 1.48× and 2.78× across six further test cases, not used for training, that vary the initial configuration, domain size, volume fraction, friction coefficient, particle–particle interactions, and imposed flow type. These results indicate that an accurate learned initializer can accelerate the original solver without altering its governing equations or convergence criterion.
Supercritical carbon dioxide (sCO2) has garnered significant interest in advanced thermal power cycles due to its high efficiency and compact system design. Nonetheless, under heating conditions, sCO2 displays highly intricate thermo-fluid behavior, especially in the proximate wall region. As the fluid traverses the Widom line, abrupt variations in thermophysical properties induce pseudo-boiling and other transport pseudo-phenomena, thereby markedly modifying turbulence structures and heat transfer mechanisms. Precise modeling of these effects necessitates fine near-wall resolution and sophisticated turbulence models, which consequently elevate computational costs. For large-scale systems, high-fidelity simulations become time-consuming and often impractical. To address these limitations, the present study develops a supervised machine learning (ML) framework aimed at predicting near-wall characteristics of heated sCO2 flow within horizontal tubes. A comprehensive numerical dataset was generated employing an Eulerian approach coupled with the Shear Stress Transport (SST)k−ω turbulence model, encompassing heat fluxes of 25-100 kW/m2, mass fluxes of 200-1000 kg/(m2.s), pressures of 7.5-20 MPa, and tube inner diameters of 4-12 mm. Various learning algorithms, loss functions, and optimizers were systematically evaluated to identify the optimal model architecture. The results demonstrate that the proposed modified deep residual ANN, incorporating residual connections, SiLU activation, a combined Huber–MAE loss function, and AdamW optimization, achieves significantly improved statistical and physical performance compared to the other ML models. The model attains a coefficient of determination (R2) exceeding 0.97 for both upper and lower wall temperatures, while accurately capturing localized thermal peaks associated with heat transfer deterioration in horizontal tubes with sCO2.
This paper presents a design optimization for zero-thickness corrugated airfoils realizable only in computational domains. To isolate the aerodynamic effects of surface corrugation from those of finite thickness – particularly the influence of leading-edge geometry on flow separation and vortex formation – a zero-thickness configuration was adopted. A multi-objective optimization approach based on the non-dominated sorting genetic algorithm II was employed. Aerodynamic performance was evaluated via a computational fluid dynamics based on the immersed-boundary method implemented on a multi-layer Cartesian grid that supports numerical calculations of zero-thickness objects. The objectives were to minimize drag and maximize lift. The airfoil geometry was parameterized using nine evenly spaced control points along the chord. Two optimization cases were investigated with two angles of attack (0°, 4°) and Re = 104. The drag–lift trade-off was examined using flow visualization and parallel coordinate plots, which revealed contrasting trends in leading-edge control parameters between drag minimization and lift maximization. Flow field visualizations revealed that lower-drag solutions typically exhibited pronounced corrugation patterns, whereas higher-lift solutions resembled cambered geometries. Low-drag designs show distinct corrugation, reducing drag below that of flat plates, as recirculation in concavities forms airfoil-like streamlines and reduces viscous drag. High-lift designs have smooth cambered shapes and trailing-edge dips to trap vortices. Further, trailing-edge dips can increase lift by trapping vortices. At a high angle of attack, the highest-lift solution exhibited flow separation and vortex formation near the leading edge. These findings suggest that corrugated airfoils are primarily advantageous for drag reduction, rendering them particularly effective for applications prioritizing low-drag performance.
The atomization of a liquid jet in crossflow is revisited for the density ratio corresponding to conventional Jet A-1 aviation fuel by numerical simulations. To correctly capture the large range of spatial scales involved, the coupled methods of momentum conserving Volume-of-Fluid and Lagrangian Particle Tracking methods are used. Both solvers are combined via momentum coupling and two-way droplet conversion. Additionally, we take advantage of Adaptive Mesh Refinement to save computation time. The Lagrangian method is free from coalescence and breakup models, therefore a global criterion is applied to avoid droplet transfer from the VOF solver to the Lagrangian solver in the region with a high number of coalescence events. To achieve this, a new post-processing method for identifying coalescence and disintegration of droplets is described. The Lagrangian method is validated regarding single droplet trajectories, and a qualitative validation of the Euler–Lagrange conversion is presented. Two different sets of inflow velocities for the gaseous and liquid phases are investigated. The results obtained with the coupled Euler–Lagrange method show a good agreement with those from pure Eulerian VOF simulations, regarding droplet size and velocity distributions and column breakup point. Jet penetration is validated, and the Sauter Mean Diameter is close to experimental correlations. A significant reduction of computational cost is achieved by the coupling to the Lagrangian method.
Reduced-order simulation of turbulent particle-laden flows is required in numerous configurations, where the resolution of the whole spectrum of turbulent scales through DNS is out of reach. Whereas structural or stochastic models have been derived in order to provide a synthetic turbulent model for the non-resolved scales of the fluid flow field, reproducing particle dynamics is challenging because it requires capturing both spatial and temporal correlations. We present a reduced-order framework that combines wavelet-based structural modelling with stochastic evolution. Using compactly supported divergence-free wavelets within a multiresolution analysis, the method provides direct control over spatial structures and correlations of synthetic multiscale incompressible velocity fields. In contrast to Fourier modes, the wavelet basis functions are localized in space and spectrally non-sharp in Fourier space, and spread over a range of wave numbers, which requires a dedicated procedure to enforce a prescribed turbulent energy spectrum. The stochastic evolution of wavelet coefficients further ensures consistent temporal correlations. The proposed framework is evaluated in homogeneous isotropic turbulence under a fully reduced setting, where all turbulent scales must be provided by the model. When coupled to a disperse phase in the one-way coupled framework, results show that it reproduces particle preferential concentration across a wide range of Stokes numbers as well as the pair-dispersion regimes, achieving similar agreement with DNS data as for classical Fourier-based Kinematic Simulation. This establishes a physically consistent turbulence model, which combines structural fidelity with stochastic dynamics, providing an alternative framework for synthetic turbulence modelling and investigating particle–turbulence interactions resolution.
This work introduces R13-ML, a machine-learning–enhanced regularized 13-moment framework for shock-dominated rarefied gas flows in the transition regime. The analytically derived higher-order closures and collisional production terms of classical R13 are replaced by fully connected neural networks trained on high-fidelity Direct Simulation Monte Carlo (DSMC) data. The networks learn mappings for the third- and fourth-order moments (mijk,Rij) and the production terms (Qij,Qi) from a physics-based, ten-dimensional feature set, combining normalization by local mean free path and sound speed with Gaussian filtering, polynomial interpolation, and a weighted Huber loss. The trained surrogates are embedded in a high-order discontinuous Galerkin spectral element solver for dynamic online evaluation of the closures. Within the one-dimensional training range Ma=1.2–8.0, R13-ML reproduces higher-order moments and collision integrals in close agreement with DSMC, while errors of the analytical R13 closures exceed an order of magnitude for Ma≳5. For one-dimensional shocks, the model remains accurate at Ma=9 and stable at Ma=12, where classical R13 develops severe oscillations. Robust generalization is further demonstrated for unsteady shock-interaction problems and for two-dimensional rarefied flow over a cylinder at Ma=7, where R13-ML captures the spatial topology and peak values of higher-order moments with stress and heat-flux deviations below about 6% relative to DSMC.
Studying the near-field jet/wake is essential for contrail formation because it sets the ice crystal number and persistence. In contrail simulations, two large eddy simulation (LES) formulations are common: a spatial approach, which explicitly resolves spatial development but is computationally intensive, and a temporal approach, which assumes the frozen turbulence hypothesis and is computationally efficient. Both approaches have been utilized in prior contrail studies; however, comparisons between them are limited. In this paper, temporal and spatial formulations of LES were compared to analyze near-field jet contrails under cruise conditions representative of an Airbus A320neowith a LEAP-1A engine. The wake vortex was then initialized from both spatial and temporals jets at two plume ages (tj = 0.12 s and tj = 0.5 s) corresponding to the pre-development and developed microphysical stages, respectively. Microphysics were treated with an online-coupled Lagrangian scheme to evaluate ice activation and growth for (i) soot-only and (ii) soot+ambient nuclei. In the jet phase, the temporal jet produced ice microphysical properties that overall were equivalent to those produced by the spatial jet under soot-rich plume and supersaturated conditions; differences observed were within the same order of magnitude for both the mean ice radius and the ice number-based emission index. In the vortex phase, the ice number concentration was inherited from the jet phase, and differences between the temporal and spatial jet initializations diminished with the plume age. Delaying the vortex onset led to relatively smaller particle radii and a slightly higher ice number concentration. Overall, ambient aerosols increased the ice number but slowed crystal growth via water vapor competition, thereby reducing any sensitivity to vortex initialization.
This paper generalizes the two-dimensional (2D) subcell volume-of-fluid (SVOF) method proposed by Sha and Jia (Comput. Math. Appl. 160, 86-107, 2024) to three dimensional (3D) cases. On this basis, we develop a staggered multimaterial arbitrary Lagrangian–Eulerian (MMALE) method for hexahedral grids. The present 3D SVOF method operates by dividing each interfacial cell into eight subcells. It then determines the normal vector of the linear interface through a least-squares minimization of the discrepancy between the reconstructed and given volume fractions of the reference fluid across these subcells. The plane constant is computed by enforcing the given volume fraction of the reference fluid in the entire interfacial cell. This 3D SVOF method maintains a key characteristic inherited from the 2D SVOF method that it does not rely on information from adjacent cells. However, in contrast to the 2D SVOF method that uses non-overlapping triangular subcells, the subcells in our 3D method are overlapping. The developed MMALE method follows a Lagrangian-plus-remap formulation. Its Lagrangian stage employs a compatible finite volume discretization, enhanced with edge-based artificial viscosity for shock capturing and Tipton’s pressure relaxation model for interfacial cell closure. The remap stage is performed via a polyhedron subdivision and intersection algorithm. A comprehensive numerical assessment is conducted to evaluate the proposed methods, encompassing static interface reconstruction, dynamic advection, and MMALE tests. Results validate that the 3D SVOF method attains second-order accuracy in static tests, outperforming some conventional volume of fluid (VOF) methods and matching the moment of fluid (MOF) method in precision. Dynamic tests further illustrate its capability to accurately capture interfaces undergoing large deformation. Finally, MMALE tests confirm its robustness, with the SVOF method producing smoother interfaces and fewer fragments than the MOF method.
The multiple-relaxation-time (MRT) lattice Boltzmann model is generally limited to high viscosity, which restricts its application to simulate high Reynolds flow on relative coarse mesh. This study investigates the numerical stability and accuracy of the two-dimensional MRT lattice Boltzmann model in simulating high-Reynolds-number flows with very low viscosity. The D2Q9 lattice with the lattice sound speed cs=1/3 is used in this paper. Particular attention is given to the construction of transition matrix, which is usually obtained by a Gram-Schmidt orthogonalization of several chosen moments. A linear stability analysis for two-dimensional MRT lattice Boltzmann model was conducted to examine the effects of orthogonalization strategies and relaxation parameters. Three orthogonalization strategies are compared: no orthogonalization, classical unweighted orthogonalization and weighted orthogonalization. The results reveal that orthogonalizing the transformation matrix is essential for stable simulations at low viscosities. The linear stability analysis suggests that acoustic modes are mainly affected by se, which is associated with the bulk viscosity, and by the ghost relaxation rate sϵ, whereas shear modes are governed by sν, which determines the kinematic shear viscosity, and by the ghost relaxation rate sq. Weighted orthogonalization significantly enlarges the feasible range of sϵ, whereas the dependence of sq on sν is intrinsic and cannot be removed. Based on these findings, feasible relaxation parameter regions were identified, providing practical guidelines for stable D2Q9 MRT simulations. The results were further validated through a series of benchmark problems. The results demonstrate that MRT reproduces the main flow features with the convergence rate very close to 2. Furthermore, weighted orthogonalization enabled stable simulations at viscosities as low as 10−5 and Reynolds numbers up to 104. These findings highlight the potential of MRT with a weighted orthogonalized transition matrix as a robust and computationally efficient approach for modeling low-viscosity, high-Reynolds-number flows in practical applications.
An immersed boundary method (IBM) with displaced forcing for the simulation of particle-laden flows is presented. The method is based on the direct forcing IBM in which the regularized Dirac delta function is used as a kernel when transferring information between the Eulerian frame for a fluid and the Lagrangian frame for a particle through velocity interpolation and force spreading. Improvement in accuracy is made by spreading the fluid–solid interaction force not around Lagrangian points marking the immersed boundary but around displaced forcing points. The forcing points are displaced inwardly from the immersed boundary to reduce the fraction of the force applied to a fluid around a particle, with the majority of the spread force being placed inside a fictitious particle domain, while the force term is formulated to impose the no-slip and no-penetration condition at the immersed boundary. A semi-implicit forcing iteration scheme is used to impose the no-slip and no-penetration condition effectively with the displaced forcing. The method is tested for several particle-laden flows. Results show that the method leads to improved accuracy in the prediction of drag and of a velocity field near the immersed boundary. It is also found that the displaced forcing IBM captures subtle fluid–particle interactions associated with the instabilities of a settling spherical particle in the oblique regimes, retaining the performance of the direct forcing IBM for moving particles. A new equation for torque acting on a particle is presented and, by enhancing the angular momentum conservation, shown to improve the prediction of angular velocities of the particle in the oblique regimes.
To achieve high-fidelity direct numerical simulation of forced isotropic turbulence that shows consistency with the Kolmogorov theory, two issues associated with forcing and time-integration of the Navier–Stokes equations are addressed in the present work. The proposed forcing extends a well-known approach in which a force proportional to fluid velocity is applied to low-wavenumber modes. It is challenging to achieve appropriate resolution at both ends of the resolved wavenumber range. We constrain the kinematic viscosity of the fluid so that the simulated turbulence adjusts itself until the required small-scale resolution is obtained. Results show that conventional restriction of forcing to the lowest wavenumbers amplifies the influence of periodic boundary conditions, causing significant spatial correlations in the velocity field. We systematically study this issue and provide guidelines concerning the selection of forcing parameters, ensuring that both longitudinal and transverse velocity autocorrelation functions take qualitatively correct forms. In addition to that, we indicate that the round-off errors following from the discrete Fourier transform can make specific modes cease to satisfy the Hermitian symmetry, causing the development of an undesired velocity field. Assuming appropriate spatial resolution, it turns out that careful selection of forcing parameters is important for obtaining results consistent with statistical properties expected for homogeneous and isotropic turbulence, whereas effective mitigation of round-off errors preserves the Hermitian symmetry, ensuring the solution remains stable over long periods of time.
A novel integrated numerical framework is developed for the direct numerical simulation of compressible two-phase flows with phase change. The framework combines three core components: (i) a low-Mach compressible solver that maintains accuracy in the low-Mach limit, (ii) a sharp-interface phase-change model utilizing the Level Set (LS) method and Ghost Fluid Method (GFM), where the mass transfer rate is determined by the interfacial heat flux, and (iii) an efficient mesh adaptation strategy based on Multiresolution (MR) analysis. The solver is validated against a series of benchmarks of increasing complexity, including an oscillating water column, Rayleigh–Taylor instability, bubble expansion under time-varying pressure, and classical phase-change problems (Stefan problem, sucking problem, and vapor bubble growth). The results demonstrate the solver’s robustness in handling high density ratios and large Jakob numbers. The performance of the MR-based adaptation is evaluated by varying the threshold parameter ɛ, which governs the trade-off between solution accuracy and mesh compression. The MR approach yields smaller L1 errors than a narrow-band AMR strategy at comparable compression rates for the Rayleigh–Taylor instability test case, demonstrating the advantage of wavelet-based refinement over purely geometric criteria. For the evaporating bubble case, substantial mesh compression is achieved while preserving accuracy comparable to that of the uniform-grid reference. These findings highlight the potential of the proposed framework for the accurate and efficient simulation of two-phase flows with phase change in the low-Mach regime.
A prediction of the reentry trajectory of a spacecraft end-of-life confers a growing importance on the accurate prediction of aerodynamic characteristics of spacecraft covering all flow regimes. In the process of spacecraft re-entering the atmosphere, the flow regimes around them cover the molecule flow regime, the rarefied regime, the slip flow regime and the near continuum regime. Due to the great differences between the geometric scales of different parts, there will be plural phenomena with multiple flow regimes in the flow field. This phenomenon will cause many new and more complex multiscale hypervelocity aerodynamic problems. But traditional numerical simulation methods, such as Navier-Stokes equation solvers and the DSMC method, are difficult to simulate multiscale and multiple flow regimes integratively. On the other hand, Boltzmann equation in kinetic theory can describe the transportation process in molecule flow regime, rarefied regime, slip flow regime and near continuum regime, solving Boltzmann equation is a felicitous method for all flow regimes in a sense. Therefore, based on the Gas Kinetic Unified Algorithm (GKUA) by solving the computable model of Boltzmann equation, the aerodynamic characteristics of simplified debris models from disintegrated reentry Tiangong-type spacecraft are simulated and compared with experimental results. The results show that the max error between GKUA and the experimental data is less than 6.1%, and GKUA can well simulate the aerodynamic force when Kn is extremely small. In the experimental state, when the axial spacing of simplified multi-body is more than twice the characteristic length, Type III shock wave interference formed by two detached bow shock waves will not impact the surface and have little influence on the normal force; It is further verified that Type IV shock wave interference form high pressure and heat flux on the surface. This study can provide an effective way to simulate the aerodynamic characteristics of spacecraft reentry and its disintegrated debris around multiple flow regimes, and provide support for predicting disintegrated debris footprint.
A liquid jet in crossflow penetration model is adapted from the literature to a fully Euler–Lagrange dispersed-phase framework for the two-phase simulation of aeronautical burners. To this end, specific Lagrangian elements are used to mimic the liquid jet body presence taking into account the flattening of the column and change of trajectory under aerodynamic loading. In addition, conventional Lagrangian droplets are emitted from the jet to reproduce surface and column breakup following the boundary layer stripping model and experimental correlations. An adequate statistical strategy is proposed to homogenize the volume associated to these emitted droplets so as to smooth the injection process while limiting the potential coupling issues linked to the Eulerian discretization. In that respect, a point source correction is implemented to overcome the large size of the liquid jet elements compared to the Eulerian cells. Although semi-experimental since the droplet diameter distribution of the injected spray is to be known a priori, the model is implemented in the simulation code AVBP and tested on a high pressure experimental bench. The model is seen to recover better fuel distribution than direct spray injection approaches usually used in such a context while keeping computational cost low.
Mid-fidelity quick-turnaround CFD methods are attractive for aerodynamics characterizations. Reduced-order techniques such as immersed boundary models of bluff bodies and actuator-line models for propeller or rotor blades are of high interest, and these models can be run highly efficiently using structured-block Cartesian GPU-accelerated solvers. In particular, the Lattice-Boltzmann Method (LBM), an alternative fluid simulation method on Cartesian grids, has shown strong performance on GPUs. Yet, comparisons of the LBM with traditional Cartesian finite-volume Navier–Stokes methods on GPUs are limited. In the present work, the accuracy and computational performance of LBM and finite-volume Navier–Stokes methods was assessed for mid-fidelity aerodynamic applications, including an isolated actuator-line rotor and a ship airwake characterization. Both solvers were run using identical computational grids and GPU hardware, allowing for a controlled comparison between the two methods. The results showed that LBM was generally much less dissipative than a 2nd-order accurate finite-volume method and provided speedups of around 4–6X in time to solution.
In this study, we numerically evaluate a characteristic explicit pressure calculation method for the improved lattice kinetic scheme (LKS) for two-phase flows. The explicit scheme utilizes a pressure evolution equation and incorporates inner iterations to increase computational stability. For the comparative analysis, we introduce an implicit approach in which the pressure is calculated by solving the Poisson equation using the simplified marker and cell (SMAC) method. The primary objective is to elucidate the relationship between the number of inner iterations and the acoustic error in the improved LKS, which is verified through a comparative analysis with the implicit method. We perform benchmark simulations under identical conditions for three canonical two-phase flow problems: (i) a stationary droplet, (ii) a single rising bubble, and (iii) Rayleigh–Taylor instability. Our results reveal that the number of inner iterations in the improved LKS affects not only numerical stability but also pressure oscillations at the acoustic scale. Specifically, the frequency of the acoustic waves depends on the number of inner iterations. Particularly under external forces such as gravity, the use of many iterations leads to a lower frequency and faster decay of acoustic waves. The implicit scheme successfully suppresses such acoustic wave scale effects, resulting in smoother pressure evolution; however, in terms of computational efficiency, the explicit method significantly outperforms the implicit method. These findings elucidate the trade-offs among numerical stability, acoustic wave errors, and computational cost in two-phase simulations.
