Unstructured Solver Research Articles

Evaluation on different volume of fluid methods in unstructured solver under the optimized condition

We compared the accuracy of volume of fluid (VOF) methods in unstructured solvers using the following five different methods: 1 - the algebraically compressive VOF method, 2 – simple coupled VOF method with Level Set (S-CLSVOF) method, 3 - interface-compressing VOF method incorporated with Laplacian filter (VOFL), 4 - isoAdvector method, and 5 - isoAdvector method incorporated with Laplacian filter (isoAdvectorL) by incorporating them into OpenFOAM®, an open-source software. To evaluate these methods under proper conditions, we compared the calculation accuracy using the optimized parameters, which are explored by Bayesian optimization. The test cases for advection accuracy of volume fraction and for imbalance of surface tension force in static multiphase fluid fields were considered. In this study, we found that the compression parameters and maximum Courant number should be adjusted to obtain high accuracy simulation according to the simulation condition in VOF and S-CLSVOF method. In VOFL and isoAdvectorL methods, the spurious current can be extremely reduced, which means that these methods are suitable for slow flow with higher Laplace number conditions.

Towards Wall-Resolved Large-Eddy Simulation of the High-Lift Common Research Model

In the present study, we performed a sixth-order wall-resolved large-eddy simulation (WRLES) of the high-lift Common Research Model (CRM-HL) on the Leadership Class Computing Cluster Summit. During the 4th High-Lift Prediction Workshop (HLPW-4), wall-modeled LES (WMLES) and hybrid Reynolds-averaged Navier–Stokes (RANS)/LES approaches showed promise in predicting the maximum lift and large flow separations near stall. In those simulations, however, laminar and transition regions were not resolved since the flow was assumed to be fully turbulent. If one wants to resolve these regions, the only viable approach is WRLES. The main purpose of the present study is to demonstrate the feasibility of WRLES for a complex real-world configuration at the wind tunnel Reynolds number. A high-order unstructured mesh LES solver called hpMusic was employed in the study. The simulation at solution polynomial order P=5 has more than 14 billion degrees of freedom per equation. Computational results are compared to experimental data in the form of surface oil flows and pressure coefficient profiles. We highlight lessons learned from attempting a grand-challenge type large-scale simulation: 1) such a simulation is feasible but very expensive, estimated to be at least three orders more expensive than WMLES, 2) oil flows produced with WRLES agree very well with experimental oil flows in fine detail, not observed in those produced with WMLES, and 3) flow visualization is a severe bottleneck, and parallel postprocessing capability is needed.

Large Eddy Simulation of Laminar and Turbulent Shock/Boundary Layer Interactions in a Transonic Passage

Abstract Shock/boundary layer interactions (SBLI) are a fundamental fluid mechanics problem relevant in a wide range of applications including transonic rotors in turbomachinery. This paper uses wall-resolved large eddy simulation (LES) to examine the interaction of normal shocks with laminar and turbulent inflow boundary layers in transonic flow. The calculations were performed using GENESIS, a high-order, unstructured LES solver. The geometry created for this study is a transonic passage with a convergent-divergent nozzle that expands the flow to the desired Mach number upstream of the shock and then introduces constant radius curvature to simulate local airfoil camber. The Mach numbers in the divergent section of the transonic passage simulate single stage commercial fan blades. The results predicted with the LES calculations show significant differences between laminar and turbulent SBLI in terms of shock structure, boundary layer separation and transition, and aerodynamic losses. For laminar flow into the shock, significant flow separation and low-frequency unsteadiness occur, while for turbulent flow into the shock, both the boundary layer loss and the low-frequency unsteadiness are reduced. The time period of the unsteadiness was hypothesized to align with the time it takes for turbulent structures to convect from the shock to the trailing edge and acoustic disturbances to travel back from the trailing edge to the shock, and this is supported by the LES results.

Combination of Advanced Actuator Line/Disk Model and High-Order Unstructured Finite Volume Solver for Helicopter Rotors

In the research field of rotorcraft aerodynamics, there are two fundamental challenges: resolving the complex vortex structures in rotor wakes and representing the moving rotor blades in the ambient airflow. In this paper, we address the first challenge by utilizing a third-order unstructured finite volume solver, which exhibits lower numerical dissipation than its second-order counterpart. This allows for sufficient resolution of small vortex structures on relatively coarse meshes. With this flow solver, the second challenge is addressed by modeling each rotor as an actuator disk (i.e., the actuator disk model (ADM)) or modeling each blade as an actuator line (i.e., the actuator line model (ALM)). Both of the two models are equipped with an improved tip loss correction, which is introduced in detail in the methodology section. In the section of numerical experiments, the numerical convergence properties of the two types of solvers have been compared in the case of two-dimensional infinite wing. In addition, the relationship between the ALM and the lifting line theory is discussed in the cases of fixed-wing calculations. Another goal of these cases is to validate the tip loss correction presented. The validation of the ALM/ADM and comparisons of computational efficiency are also demonstrated in simulations involving both hover and forward flight rotors. It was found that the combination of the third-order finite volume solver and the ALM/ADM with the improved tip loss correction presents an efficient way of performing the aerodynamic analysis of rotor-induced downwash flow.

A comprehensive CFD investigation of tip vortex trajectory in shrouded wind turbines using compressible RANS solver

It is well known that a shroud placed around a wind turbine can increase its power coefficient, but it brings complex mechanisms by which the shroud alters the flow passing through the rotor. Such mechanisms impose numerical challenges, as the shrouded turbines present nonlinear behavior in the wake. This paper deals with a comprehensive analysis of tip vortex trajectory in shrouded wind turbines using Reynolds Averaged Navier–Stokes numerical solutions. The analysis includes aerodynamic performance and vortex characteristics of the whole wind turbine. The Multiple Reference Frame is used on a high-order unstructured compressible solver to study both, isolated and shrouded rotor. The NREL Phase VI Unsteady Aerodynamic Experiment rotor is used as a test case. The accuracy of results for wind speeds between 7 and 25 ms−1 is discussed. Overall, good agreement is achieved between the computed pressure distributions and the experimental reference values. At stalled blade, more efforts are needed to improve numerical solutions, especially for integrated load quantities. The vortex structure is examined, showing that shroud impacts tip vortex trajectory by the increase of the axial induced velocity at the rotor plane. This result, demonstrates that the classical Prandtl tip loss is not accurate for shrouded turbine analysis, and modern finite blade functions are needed. The influence of the flow conditions on the tip vortex trajectory, flow separation and shroud interaction are also discussed.

Parallel optimization and application of unstructured sparse triangular solver on new generation of Sunway architecture

Large-scale sparse linear equation solver plays an important role in both numerical simulation and artificial intelligence, and sparse triangular equation solver is a key step in solving sparse linear equation systems. Its parallel optimization can effectively improve the efficiency of solving sparse linear equation systems. In this paper, we design and implement a parallel algorithm for solving sparse triangular equations in combination with the features of the new generation of Sunway architecture, and optimize the access and communication respectively for 949 real equations and 32 complex equations in the SuiteSparse collection. The solution efficiency of the algorithm presented in this paper outperforms the cuSparse algorithm on NVIDIA V100 GPU platforms in more than 71% of the cases, and the speedup is even better in solving larger cases (matrix size greater than 10,000): our method increases the speedup from 1.29 time of the previous version to an average speedup of 5.54 and the best speedup of 32.18 over the sequential method on the next generation of Sunway architecture when using 64 slave cores.

Calibrating sub-grid scale models for high-order wall-modeled large eddy simulation

AbstractHigh-order methods have demonstrated orders of magnitude reduction in computational cost for large eddy simulation (LES) over low-order methods in the past decade. Most such simulations are wall-resolved implicit LES (ILES) without an explicit sub-grid scale (SGS) model. The use of high-order ILES for severely under-resolved LES such as wall-modeled LES (WMLES) often runs into robustness and accuracy issues due to the low dissipation embedded in these methods. In the present study, we investigate the performance of several popular SGS models, the static Smagorinsky model, the wall-adapting local eddy-viscosity (WALE) model and the Vreman model, to improve the robustness and accuracy of under-resolved LES using high-order methods. The models are implemented in the high-order unstructured grid LES solver called hpMusic based on the discontinuous flux reconstruction method. The length scales in these SGS models are calibrated using the direct numerical simulation (DNS) database for the turbulent channel flow problem. The Vreman model has been found to produce the most accurate and consistent results with a proper choice of the length scale for WMLES.

Development and assessment of an actuator volume method in rotating frame for predicting the flow-field of horizontal-axis wind turbines

The prediction of the unsteady flow-field in horizontal-axis wind turbines using blade-resolved CFD simulations is challenging and computationally demanding. The simulations become more affordable when using a representation of the blades based on the generalized actuator disc concept; however, the projection of the computed forces in the computational domain introduces additional mesh parameters and requires codes that impact the computation time in unstructured solvers and parallel computing. In this work, a new approach to model the effects of the rotating blades using source terms and a sliding mesh technique is proposed. The developed framework also includes the tower and nacelle. The force distributions are obtained using an in-house code, which incorporates Buhl’s correction, three-dimensional effects, and tip and hub loss corrections. The validation was performed against blade-resolved simulations and experimental data at three tip-speed ratios. At design conditions, both models are in good agreement in terms of velocity deficit and angle of attack, being underestimated at off-design conditions. The local effects of the actuator volumes, the induction region, and the near-wake have been analysed in detail. This new approach is flexible, robust, and provides good accuracy, whereas the computational cost is considerably lower than high-fidelity simulations.

ExaWind: Open‐source CFD for hybrid‐RANS/LES geometry‐resolved wind turbine simulations in atmospheric flows

AbstractPredictive high‐fidelity modeling of wind turbines with computational fluid dynamics, wherein turbine geometry is resolved in an atmospheric boundary layer, is important to understanding complex flow accounting for design strategies and operational phenomena such as blade erosion, pitch‐control, stall/vortex‐induced vibrations, and aftermarket add‐ons. The biggest challenge with high‐fidelity modeling is the realization of numerical algorithms that can capture the relevant physics in detail through effective use of high‐performance computing. For modern supercomputers, that means relying on GPUs for acceleration. In this paper, we present ExaWind, a GPU‐enabled open‐source incompressible‐flow hybrid‐computational fluid dynamics framework, comprising the near‐body unstructured grid solver Nalu‐Wind, and the off‐body block‐structured‐grid solver AMR‐Wind, which are coupled using the Topology Independent Overset Grid Assembler. Turbine simulations employ either a pure Reynolds‐averaged Navier–Stokes turbulence model or hybrid turbulence modeling wherein Reynolds‐averaged Navier–Stokes is used for near‐body flow and large eddy simulation is used for off‐body flow. Being two‐way coupled through overset grids, the two solvers enable simulation of flows across a huge range of length scales, for example, 10 orders of magnitude going from O(μm) boundary layers along the blades to O(10 km) across a wind farm. In this paper, we describe the numerical algorithms for geometry‐resolved turbine simulations in atmospheric boundary layers using ExaWind. We present verification studies using canonical flow problems. Validation studies are presented using megawatt‐scale turbines established in literature. Additionally presented are demonstration simulations of a small wind farm under atmospheric inflow with different stability states.

Evaluating performance portability of five shared-memory programming models using a high-order unstructured CFD solver

This paper presents implementing and optimizing a high-order unstructured computational fluid dynamics (CFD) solver using five shared-memory programming models: CUDA, OpenACC, OpenMP, Kokkos, and OP2. The study aims to evaluate the performance of these models on different hardware architectures, including NVIDIA GPUs, x86-based Intel/AMD, and Arm-based systems. The goal is to determine whether these models can provide developers with performance-portable solvers running efficiently on various architectures. The paper forms a more holistic view of a high-order solver across multiple platforms by visualizing performance portability (PP) and measuring productivity. It gives guidelines for translating existing codebases and their data structures to these models.

Verification of the γ-Reθt Transition Model in OVERFLOW and FUN3D

The findings of the recent AIAA and NATO workshops on transition modeling have underscored the importance of verifying the transition models that are based on Reynolds-averaged Navier–Stokes. To that end, the present work aims to verify the Langtry–Menter (LM2009) transition model coupled with Menter’s shear-stress transport (SST) turbulence model. Specifically, the SST-2003-LM2009 model, as implemented in the OVERFLOW structured overset solver and the FUN3D unstructured solver, is applied to both natural and bypass transition scenarios on a flat plate and to natural and separation-bubble-induced transition on the NLF-0416 airfoil. The Richardson extrapolation method was used to analyze and report the grid convergence of integrated and local load coefficients upstream, within, and downstream of the transition region, based on solutions obtained using the two solvers on identical grids. Reasonable agreement is obtained between the two codes for both global and local surface load coefficients when further refining the workshop mesh series, although with a grid-convergence index for the local friction coefficient that may exceed 1% in certain cases. The study also highlights the solution sensitivity to the choice of turbulence inflow boundary conditions and model options, along with the necessity to unambiguously specify this information for benchmark exercises.

Analysis of uncertainty factors in the prediction of missile pitch damping derivatives by unstructured CFD solver

Based on the transient planar pitching method, the dynamic derivative calculation function is implemented on an unstructured computational fluid (CFD) software NNW-FlowStar. Then the pitch-damping dynamic derivatives of the standard model of the ANF missile at zero angles of attack and different Mach numbers (0.5∼4.5) are calculated. The influence of various numerical and physical parameters on the calculation results of dynamic derivatives is studied in detail. It is found that the dynamic derivatives of pitch damping are independent of the values of amplitude (0.125° ≤ A ≤ 1°) and reduction frequency (0.025 ≤ k ≤ 0.2). When the number of calculation steps of the sine oscillation period is 200 and the number of inner iterations is 40, the prediction accuracy and computing efficiency can meet the engineering requirements simultaneously in this case. The influence of different grid resolutions on the dynamic derivative calculation results is studied. Results show that the grid requirements of dynamic derivative calculation are consistent with those of steady-state calculations. Finally, the calculation results are compared with the free flight test results, which verifies the reliability of the methods and parameters suggested in this paper.

A massively parallel implicit 3D unstructured grid solver for computing turbulent flows on latest distributed memory computational architectures

An implicit unstructured grid density-based solver–PRAVAHA based on a parallel variant of the lower upper symmetric Gauss-Seidel (LUSGS) method is developed to compute large-scale engineering problems. A four-layered parallel algorithm is designed to efficiently compute three-dimensional turbulent flows on massively parallel modern multiple instruction multiple data-stream (MIMD) computational hardware. This data parallel approach achieves multiple layers of parallelism including continuity of flow solution, transfer of solution gradients, and calculation of drag/lift/solution residuals, right up to the innermost implicit LUSGS solver sub-routine, which is relatively less explored in the literature. Domain decomposition is performed using the METIS software based on multi-level graph partitioning algorithms. Non-blocking message-passing interface library functions are used to manage inter-processor communication through explicit message passing, efficiently. Super-linear scalability of the parallel solver is established on the current state-of-the-art supercomputing facility, the 838 teraflops PARAM seva on up to 6144 cores. Linear or even super-linear speedup on problems of significant size is observed even on ad-hoc parallel computing platforms like workstations and multi-node clusters, for turbulent flow simulations.

Cost vs Accuracy: DNS of turbulent flow over a sphere using structured immersed-boundary, unstructured finite-volume, and spectral-element methods

We report a comparative study of three numerical solvers for the direct numerical simulation of the flow over a sphere at Re = 3700. A high-order spectral-element code (Nek5000), a general purpose, unstructured finite-volume solver (OpenFOAM) and an in-house Cartesian solver using the immersed-boundary method (IBM) are employed for the analysis; results are compared against previous numerical and experimental data. Numerical results show that Nek5000 and the IBM code operate within a similar computational performance range, in terms of cost-vs-accuracy analysis, for both global parameters as well as local flow features. On the other hand, OpenFOAM needed a significantly higher number of degrees of freedom (and, overall, a higher cost) to match some of the basic features of the flow, such as the length of the recirculation bubble forming downstream the sphere. For the finest grid resolutions, the three codes are in good agreement for most of the analyzed flow metrics. Overall, our results suggest that high-order methods and second-order, energy-conserving approaches based on the IBM may be both viable options for high-fidelity scale-resolving simulations of turbulent flows with separation.

A parallel unstructured multi-color SOR solver for 3D Navier–Stokes equations on graphics processing units

In this paper, we develop a GPU parallelization technique to accelerate the pressure Poisson solver. We apply it in three-dimensional Navier–Stokes simulations with arbitrary geometries. The numerical method is based on a projection method, and the discretization is obtained by combining a second-order unstructured finite-volume method and a σ transformation. The pressure solution is responsible for this formulation’s most time-consuming computation tasks. This is because the Poisson equation must be approximated. This work proposes a GPU-based Multi-Color Successive Over-Relaxation (MC-SOR) method to solve the resulting sparse linear system. The proposed parallel finite-volume code is verified and validated by solving the Poisson problem. We employ analytical solutions and well-known benchmark problems with complex flows. These include the Taylor–Green vortex problem and the problem of turbulent flow around a circular cylinder.

Efficient spectral implementation of ODE wall model and the extension of integral wall model to unstructured LES solvers

While wall modeling has become an indispensable tool for scale-resolving simulation of high Reynolds numbers wall-bounded turbulent flows, discussion of efficient implementation of wall models in parallel unstructured-grid solvers is lacking. In this paper, we discuss physical and numerical improvements to two zonal wall models based on integral form of the boundary layer differential equations within an unstructured-grid cell-centered finite-volume LES solver. The first model is a novel implementation of the ODE equilibrium wall model, where the velocity profile is expressed in the integral form using the constant shear-stress layer assumption and the integral is evaluated using a spectral quadrature method, resulting in a local and algebraic (grid-free) formulation. The second model is the modified form of the integral wall model of Yang et al. (2015) [14], which is based on the vertically-integrated thin-boundary-layer PDE along with a prescribed composite velocity profile in the wall-modeled region. A modification to the inner layer velocity profile which fixes the coordinate invariance issue of the original model is proposed. Additionally, several numerical challenges unique to the implementation of these integral models in unstructured mesh environments, such as the exchange of wall quantities between wall faces and LES cells, and the computation of surface gradients, are identified and possible remedies are proposed. The computational performance of the wall models is assessed both in a priori and a posteriori settings against the traditional finite-volume based ODE equilibrium wall model, showing a comparable computational cost for the integral wall model, and superior cost scaling for the spectral implementation over the finite-volume based approach. Load imbalance among the processors in parallel simulations seems to severely degrade the parallel efficiency of finite-volume based ODE wall model, whereas the spectral implementation is remarkably agnostic to these effects. The integral wall model, however, is even more prone to parallel efficiency degradation than the finite-volume based ODE wall model due to communication overhead among processors in highly parallelized settings.

Acceleration of the data-parallel lower-upper relaxation time-integration method on GPU for an unstructured CFD solver

The acceleration of CFD solvers on GPUs allows for large speedups, which enable more complex and detailed simulations. However, not all numerical methods can be efficiently executed in a massively parallel fashion. The Runge–Kutta (RK) and Lower-Upper Symmetric Gauss–Seidel (LU-SGS) time-integration methods are widely used on GPUs and multicore CPUs, respectively; however, the RK method suffers from low convergence speed, while the LU-SGS method is not adapted to a many-core environment. In this paper, we propose to accelerate the matrix-free version of the Data-Parallel Lower-Upper Relaxation (DP-LUR) time-integration method on the unstructured solver FaSTAR, developed by JAXA. The DP-LUR method outperformed the RK and LU-SGS methods on an NVIDIA Tesla V100 GPU when tested on the ONERA M6 and NASA CRM geometries, and reached up to 3.83% of the performance peak of one device.

Computation of drag crisis of a circular cylinder using Hybrid RANS-LES and URANS models

Numerical simulation of flow past a circular cylinder across the “drag-crisis” region is extremely challenging for turbulence models because the boundary layer undergoes laminar–turbulent transition and variable-locus separation. We investigate the SA-DDES hybrid model along with two variants, namely, SA-kLES and SA-ILES, based on Spalart–Allmaras (SA) model, and include for comparison the SA-BCM transition and the SA-URANS models, for Re ranging from 50,000 to 5 million, using an in-house unstructured grid solver. All hybrid RANS-LES models produced clearly turbulent-like behavior, as evident from the Q-criterion, while the URANS models did not. A decline in the drag coefficient is noticed in all the turbulence models, but not the sharp decrease observed experimentally, with one exception: the SA-BCM transition model, which predicted the drag coefficients much closer to the experiments. The hybrid RANS-LES models outperformed the URANS SA-BCM model only in the fully turbulent trans-critical region and better represent the physics in the wake region for all Reynolds numbers studied. All the hybrid RANS-LES models produced similar results, suggesting comparatively equal performance in predicting separated flows. We believe that the performance of a hybrid model for mid-range Reynolds numbers will be greatly enhanced if the model is equipped to handle the laminar–turbulent transition.

Investigation of low-dissipation low-dispersion schemes for incompressible and compressible flows in scale-resolving simulations

A comprehensive study is conducted on a second-order low-dissipation low-dispersion (LD2) scheme in scale-resolving simulations of both incompressible and compressible flows, using a node-based unstructured CFD solver. The scheme deploys a higher order central reconstruction of the face values (up to fourth-order on structured meshes) and a matrix dissipation formulation to reduce the dispersive and dissipative numerical errors. The LD2 scheme is examined for compressible flow cases involving shock discontinuities, LD2-Compressible (LD2C), and is verified in a classical shock-tube problem. The scheme is then further verified in Large-Eddy Simulations (LES) of decaying isotropic turbulence (DIT) in comparison with available experimental data. It is shown that in scale-resolving simulations, the LD2C scheme is able to significantly improve the prediction as compared to a conventional second-order central scheme. The scheme is then further assessed and verified in hybrid Reynolds-Averaged Navier–Stokes (RANS)-LES computations for the subsonic and supersonic turbulent channel flow, where excellent agreement with reference DNS and correlations are observed. Moreover, a supersonic base flow is simulated using hybrid RANS-LES, where improved predictions are observed. The LD2C scheme exploits a shock sensor incorporating vorticity and is shown to improve the prediction of the resolved shear stress in the shear layer of compression.

Influence of isolator section on the shock augmented mixing in SCRAMJET engine

The paper aims to numerically study the effect of isolator design on shock-augmented mixing in a generic supersonic Ramjet (SCRAMJET) engine and reduce what would be a very complicated thermodynamics+aerodynamic optimization in the design to a simpler aerodynamic one. Changing the length of the isolator section effectively varies the thickness of the boundary layer that is formed in the isolator, which then enters into the SCRAMJET combustor. The study is conducted on a three-dimensional geometry using the Menter's SST model for turbulence on an in-house unstructured grid RANS solver. In this paper, the authors quantify the effect of the incoming boundary layer thickness on the shock-augment air-fuel mixing in the SCRAMJET combustor. It is observed that for the maximum variation in the isolator length, the normalized mixing volume only improves by 7%, which is only a marginally better mixing of air and fuel obtained as with an increase in the isolator length. The authors conclude that the Prandtl-Meyer expansion fan (PMEF) formed at the flame-holder effectively isolates the combustor from the isolator section to a degree that the isolator dimensions can therefore be chosen entirely on considerations of flow homogenization of inlet air before it enters the combustor section of the SCRAMJET engine. Although, the paper analyses a generic SCRAMJET with specific dimensions, the results are more generally applicable than for the specific design actually considered in the paper. Thus this study is also relevant to SCRAMJET engines that may differ radically from our considered design in form and dimensions.