Direct Numerical Simulation of Transitional and Turbulent Flows Over Multi-Scale Surface Roughness—Part I: Methodology and Challenges

Representing realistic engineering geometries within computational fluid dynamics (CFD) simulations is conventionally achieved by solving the equations of fluid motion on computational grids conforming to the shape of the body, that is, using body-conformal grids. More numerous and finer-scaled geometries demand more grid points to resolve, as will higher Reynolds numbers, which tend to demand ever finer grids near the body’s surface for high-fidelity CFD. Additionally, it is also often not clear a priori what the grid resolution requirements may be, such that generating grids for complex geometries tends to require an iterative approach [1]. Thus in the trend toward realistic geometries and flow conditions, solving for the fluid flow to a high accuracy becomes increasingly challenging with increasing geometrical complexity and increasing Reynolds number as encountered in practical engineering applications. Furthermore, given the prevalence and variety of in-service roughness found on turbomachinery components [2], accurate numerical prediction of blade performance demands faithful representation of this roughness. Adequate resolution of micro-scale features on the surface of a body is especially challenging: compared to smooth-wall or regular roughness studies, a much smaller number of direct numerical simulations (DNS) of realistic, irregular rough surfaces has been attempted [3,4]. Especially in the case of generalized, irregular roughness as introduced in the present work, and further discussed in Part II of this paper [5], it becomes impractical to generate new grids for conventional body-conformal CFD for each roughness topology, especially as different roughness parameters, such as height, or the extent of the roughness, must be varied to gain insight on the effect of the roughness on the fluid field.

Immersed boundary (IB) methods, where wall boundary conditions are enforced on grids that do not conform to the geometry of the body or solid boundary, are a valid alternative. The technique was first introduced in the 1970s to simulate cardiac mechanics with elastic boundaries and the associated blood flow as a fluid–structure interaction problem at low Reynolds number [6]. The technique was subsequently extended to model rigid bodies by the addition of artificial body forces to the fluid equations of motion in an effort to drive the velocity field near the boundary to the desired boundary value. The so-called direct forcing method [7] introduces a force distribution that depends on the velocity field, and it can be seen as an artificial response of the body necessary to directly cancel the error in this boundary value at each time-step. For this reason, an additional set of equations (often in the form of a spring-damper system) needs to be solved at every time-step, increasing the computational cost. Moreover, the numerical stability requirements of these equations often pose strong constraints on time-step size. Another drawback of these direct forcing approaches is that the resulting forcing functions are not designed to satisfy the underlying conservation laws for the cells in the vicinity of the immersed boundary, and such methods are not easily extended to general interfacial conditions such as Neumann boundary conditions [8]. In addition, direct forcing methods are in general first-order accurate in the L_∞ norm, and generate spurious pressure fluctuations in the region of the immersed body, resulting from a jump in the velocity across the boundary before the projection step [9]. This is especially problematic for turbomachinery applications where surface pressure on the body is a key quantity of interest. Strict global and local conservation of mass and momentum can only be guaranteed by resorting to a finite-volume approach near the IB, which is the motivation for cut-cell finite-volume IB methods (e.g., for two-dimensional flow [10]). Extending this approach to three dimensions is nontrivial as the cut-cell procedure leads to complex polyhedral cells and discretization of the full Navier–Stokes equations on such polyhedral cells is challenging [1]. Despite the challenges of IB methods, they have been pursued actively in the literature, since, in addition to the complication of grid generation when using body-conformal grids, Cartesian-grid solvers used with IB methods have the potential to generate solutions to complicated problems orders of magnitude faster than conventional body-fitted-grid solvers since, due to grid regularity, solvers can benefit from the use of simple, accurate, and robust numerical methods [8].

The boundary data immersion method (BDIM) was first introduced by Weymouth and Yue [8]. Its derivation is based on a general integration kernel formulation allowing the field equations of each domain (i.e., fluid and solid body) and the interfacial conditions to be combined analytically via a meta-equation. The initial formulation of the BDIM is general, therefore it can be applied to a variety of solid/fluid systems including immersed no-slip (such as in the present work) as well as free-slip bodies. By altering the analytic form of the fluid equations, the method ensures exact enforcement of the boundary data and maintains physically consistent behavior near the smoothed interface representing the solid object. Moreover, the BDIM does not incur time-step penalties as other forcing-based methods because it does not require the solution of any additional equations. This is fundamental for high-Re studies, where the time-step is already small and the overall computational cost is high. While the no-slip formulation matches the direct forcing method mentioned above closely, a prominent feature of the BDIM is an additional modification to the pressure term analogous to the discrete operator adjustments of sharp-interface methods (being another class of IB methods where the communication between the moving boundary and the flow solver is accomplished by explicitly modifying the computational stencil near the IB, and both conservation equations are altered). Thus the projection issues for direct forcing methods [9] are avoided. Whereas sharp-interface methods alter the discrete operators, and BDIM alters the analytic equations near the embedded boundary. This allows for simple implementation within existing flow solvers regardless of the geometry being simulated. Building on this initial BDIM formulation, Maertens and Weymouth [11] improved the method’s ability to tackle intermediate Reynolds numbers through the addition of a higher-order term to the integral formulation. The present formulation was extended by Schlanderer et al. [12] to compressible viscous flow.

The present work introduces a framework that exploits an IB method, the BDIM, for high-fidelity simulations of transitional and turbulent flows over multi-scale rough surfaces, with a particular focus on turbomachinery components at engine-relevant conditions. The framework is based on the original mathematical formulation of the BDIM, but includes a new algorithm, here discussed, that extends its applicability to general three-dimensional non-conforming surfaces. This paper also introduces a novel approach to compute and interpret the on-surface data obtained for general non-conforming geometries. In particular, a mathematical approach to derive the wall shear stress and the heat flux on non-smooth walls is introduced. The method is validated in the case of a turbulent channel with egg-carton roughness and on a smooth high-pressure turbine (HPT) blade at engine-relevant conditions. Finally, a new computational setup that exploits the advantages of the BDIM to perform first-of-a-kind high-fidelity simulations of a high-pressure turbine blade with roughness is presented.

The generalized mathematical formulation adopted in this work can be found in Refs. [11,12] and will only be summarized here. The algorithm used to discretize the BDIM equations in the case of general three-dimensional surfaces is introduced.

For static bodies, the weighting functions and the normal unit vectors are computed only once, at the beginning of the simulation. For moving or deformable boundaries, they have to be computed at every time iteration. Equation (1), on the other hand, is enforced at every time iteration, for both static and moving geometries.

On-surface data provide direct quantitative information on the effect of the flow on an immersed body or a rough surface on the flow. For this reason, instantaneous on-surface data are collected at regular intervals during run-time. The data collection procedure is composed of three steps, which can be summarized as follows:

• A series of on-surface data collection points based on the BDIM geometry is defined.

• The location of the collection points, due to the non-conforming nature of the solid boundary, is mapped onto the computational grid. In the case of non-moving solid bodies, as in the present work, these preliminary steps are performed once per simulation, in the pre-processing stage. For moving or deforming geometries, the computation of the collection points and the location mapping have to be repeated at every collection time-step.

• The desired flow variables of interest are then computed at the collection points at each capturing time interval by means of interpolation.

Given the generic flow variable $ψ(x→,t)$ (either a scalar or the scalar component of a vector or a tensor), we indicate with ψ_i the value of ψ corresponding to the ith triangular surface element. ψ_i can be extracted from the flow computational grid by interpolating its value at the centroid of the triangle, or by interpolating at its nodes and then computing its surface average. Both methods are presented.

The Navier–Stokes equations are numerically solved on the flow computational grid. Hence, some variables such as the wall shear stress and the heat flux, because of their dependence on the local orientation of each surface element, require an additional step to be computed.

Since roughness increases the surface area exposed to the flow with respect to an equivalent smooth geometry, it can be useful to introduce a normalization based on a reference smooth area. This is achieved by mapping each triangular element of the rough surface onto an equivalent element of the reference smooth surface by projecting its nodes, as shown in Fig. 4. Indicating with A the area of an element of the rough blade, we can therefore compute its reference area A_s obtained from the projection. Similarly, we can compute a reference tangential direction $t→1s$ defined with respect to the reference element.

The BDIM has been extensively validated in the context of two-dimensional simulations over a range of test cases of different complexity including the flow around a cylinder, turbulent boundary layers, and airfoil trailing edge noise [12]. The purpose of this work is to provide a numerical framework to investigate transitional and turbulent flows over multi-scale roughness. Hence, the validations discussed here serve to demonstrate that the numerical method is able to accurately predict the various flow mechanisms that are relevant to the topic. Schlanderer et al. [12] already demonstrated the capabilities of the BDIM to correctly simulate the transient growth of oblique instabilities in a supersonic flat plat boundary layer. The framework not only accurately captured the amplitude of the velocity, temperature, and density perturbations but also the growth rate, which is a strong endorsement for the method’s suitability for compressible transitional flow. The accuracy has been demonstrated also in the case of a plate not aligned with (and hence not conforming to) the underlying fluid computational grid. It is worth repeating that the numerical formulation of the BDIM has not changed in this work, we only propose a new point-mapping algorithm that can extend the current capabilities to simulate general complex three-dimensional geometries.

In the context of transitional and turbulent flows over complex surfaces, it is worth mentioning the recent work of Ananth et al. [14], in which the BDIM is successfully applied to assess the performance of riblets with different shapes. In this case, the BDIM has also been shown to be able to simulate the effect of a roughness element to trigger the boundary layer transition.

Here, we present a validation of the BDIM applied to two different cases: a three-dimensional channel with egg-carton roughness and a smooth HPT blade at engine-relevant conditions. In the case of the channel, we will see that the BDIM is able to correctly represent the effect of surface roughness on the boundary layer velocity profile and higher-order statistics, provided the grid in the vicinity of the wall is sufficiently refined to accurately resolve the interface between flow and solid body. This test case is also employed to validate the data-capturing framework. The smooth HPT case will serve to assess the capability of the framework to capture complex aero-thermal features such as strong pressure gradients, boundary layer transition, shock waves, and strong turbulent mixing processes. All the simulations are carried out using HiPSTAR, a well-validated in-house numerical solver for compressible flow [15].

DNS of a turbulent channel flow with roughness are performed to validate the performance of the present three-dimensional BDIM. Since the main focus of the validation is on the near-wall region, the minimal-span channel with roughness used by Macdonald et al. [16] is adopted. It has in fact been shown that the minimal-span channel can accurately capture the turbulent near-wall flow using a reduced spanwise domain size with respect to a full-span simulation, hence significantly reducing the computational cost. A schematic of the computational setup for a rough channel and a reference smooth channel is shown in Fig. 5. The dimensions of the reference smooth channel, normalized with respect to the channel half-height h, are L_x = 2258/Re_τ, L_y = 360/Re_τ, and L_z = 113/Re_τ, where x, y, and z are the streamwise, wall-normal, and spanwise directions, respectively. Re_τ is the friction Reynolds number based on the reference channel half-height h. Simulations are carried out at Re_τ = 180. The domain is periodic in the x and z directions and it is discretized using a body-fitted Cartesian grid, with a number of grid points N_x = 1130, N_y = 169, and N_z = 56. The spacing in the streamwise and spanwise directions is constant and it is equal to Δx⁺ = 2.0 and Δz⁺ = 2.0, respectively. The superscript ⁺ is employed to indicate lengths expressed in viscous units. The grid in the wall-normal direction is refined in the near-wall region, with a minimum spacing of $Δyw+=0.40$ at the wall and a maximum spacing of $Δyc+=6.1$ at the channel centerline. The flow is driven by a pressure gradient P_x = −1 acting in the streamwise direction, constant both in time and space, and it is an input to the simulation. A uniform heat sink H_s = −20 is applied to the flow, while both smooth and rough walls are iso-thermal, with a uniform temperature T = 1. The values of P_x, H_s, and T are normalized using reference primitive variables or a combination of thereof. The reference Mach number is set to M = 0.018, small enough to ensure that any compressibility effect is negligible.

The extra resolution with respect to the smooth channel is required to adequately represent the effect of the rough walls using the BDIM. In fact, the grid should offer enough resolution in the streamwise, spanwise, and wall-normal directions to resolve the BDIM smoothing region. The grid requirements therefore highly depend on the solid body geometry (in particular the local orientation of its surface with respect to the Cartesian axis) and on the characteristics of the flow. The ideal case is obviously a grid with equal uniform spacing in all directions, such that Δx = Δy = Δz, to guarantee the best representation of the soothing region. However, this is generally impractical from a computational perspective, because the resulting overall point count is often prohibitive. It is therefore necessary to find a trade-off between accuracy and number of points. Some general guidelines can be given for the specific case of wall-bounded flows, for which the direction that requires the largest resolution is the wall-normal, due to the large velocity gradients. For wall-bounded flows, $ε$ should be prescribed as a multiple of the wall-normal spacing Δy_min (uniform within the BDIM smoothing region), such that $ε=αΔymin$ and 2.0 ≤ α ≤ 2.5. This ensures that the smoothing region contains a minimum of four grid points in the wall-normal direction. For the present case, we have chosen α = 2.0. For the streamwise and spanwise directions, a convergence study should be performed to ensure that the smoothing region is represented by at least two grid points. This preliminary study can be performed just by observing the resolution of the smoothing region on the grid, without performing the simulation. This approach performed in the x–y plane for the rough channel is presented in Fig. 6 by showing the values of the first weighting function μ₀ for different streamwise resolutions. The grid resolution is progressively increased, as summarized in Table 1. By visual inspection, the coarsest grid has been discarded due to the inadequate aspect ratio of the cells at the wall. A canonical grid convergence study has then been performed on grids (c) and (d) from Fig. 6 to assess their accuracy in capturing the relevant fluid dynamics. Grid (c) with N_x = 1130 has been preferred because it ensures adequate accuracy at a lower computational cost.

Due to the higher resolution required in the near-wall region, the grid used for the rough channel simulation is 1.5 times larger than the equivalent smooth case. In terms of computational requirements, owing to the additional operations required to enforce the boundary conditions using the BDIM, the cost for the rough channel simulation was 1.64 times higher than the smooth simulation.

The results for the wall-normal mean streamwise velocity profile $u¯+$ and for the root mean square of the velocity fluctuations $urms′+$⁠, $vrms′+$ and $wrms′+$ are shown in Fig. 7 as a function of the wall-normal viscous coordinate y⁺. The velocity profiles for both the smooth and the rough channels have been averaged in time, as well as in the streamwise direction. The averaging is performed on an interval of 100 flow-through units, corresponding to the time it takes for flow to convect through the channel 100 times. The results from the current simulations agree well with the reference data obtained from body-fitted simulations by Macdonald et al. [16]. The minor differences observed fall within the accuracy range that is expected when comparing BDIM against body-fitted data and they can also be attributed to differences in the numerical solver, the grid resolution, and the simulation parameters. Similar trends have been observed in Ref. [14] for the same validation. It is also important to remember that the minimal-span channel can adequately capture only the near-wall turbulent structures up to a distance of 0.4L_z from the wall. Hence the comparison of the present method with respect to the reference is restricted to that region.

The challenges of simulating surface roughness and its effects on the performance of a turbine blade are investigated and discussed in the case of a VKI LS-89 HPT blade [17] in a linear cascade using DNS. This configuration without surface roughness has already been explored in previous studies in the context of high-fidelity simulations, thus offering a well-established reference [18–20]. Simulations are carried out at engine-relevant conditions, with a Reynolds number of Re = 590, 000 based on the axial chord C_ax, the outlet flow velocity, and an exit Mach number M = 0.9. The spanwise extent of the HPT blade is set to 0.4C_ax, which is wide enough to ensure the correct development of the largest inflow turbulence structures [20]. All quantities are presented in non-dimensional form, with lengths normalized by C_ax and flow quantities such as velocities, density, and temperature normalized by reference inlet conditions. At the inlet, synthetic turbulence is introduced in the domain employing a compressible version of the digital filter [21] method. The turbulence intensity T_k, obtained from the three components of the velocity fluctuations as u′, v′, and w′, is defined as $Tk=1/3(urms′2+vrms′2+wrms′2)$ and is equal to $8%$ of the axial inlet velocity U_in. The integral length scale is $Lx/Cax=8%$⁠. At the outlet, a non-reflective characteristics boundary condition [22] is enforced in combination with a stretching of the grid in the streamwise direction to introduce numerical dissipation and to avoid non-physical reflections of the acoustic waves.

A schematic of the computational domain in the x–y plane is shown in Fig. 8. The computational grid is composed of three blocks with different topologies interfaced using an overset method [23], which allows the exchange of information between blocks by means of interpolation in the overlapping regions. This allows to locally optimize the grid topology and resolution to suit the needs of regions of the domain subjected to different physical dynamics and different computational requirements. As shown in Fig. 8, the blade passage is discretized using an H-type Cartesian background grid, indicated as block 1. This grid contains the inlet and the outlet of the simulation and it is subjected to periodic boundary conditions in the pitchwise direction. Block 1 is discretized using N_x × N_y × N_z = 1470 × 716 × 576 grid points in the axial, pitchwise, and spanwise directions, respectively. Two O-type grids are employed to discretize the flow region around the blade. These grids are arranged in two layers, an outer (block 2) and an inner one (block 3), and their topology follows the mean blade geometry. Block 2 has a resolution of N_s × N_n × N_z = 8165 × 239 × 576 points, where the subscript s denotes the blade tangential direction and n the wall-normal direction. Block 3 has been introduced to provide additional resolution in the near-wall region, necessary to accurately resolve the interactions between roughness and turbulent eddies using the BDIM. The grid count for block 3 is N_s × N_n × N_z = 35,595 × 110 × 1200 points for the rough blade and N_s × N_n × N_z = 35,595 × 100 × 576 points for the smooth blade. Introducing a normalization based on the local viscous length scale obtained from the reference smooth blade, the grid spacing is such that Δs⁺ ≤ 3.0, Δn⁺ ≤ 0.9 and Δz⁺ ≤ 6.5, resulting in a much higher resolution with respect to previous high-fidelity simulation studies [18,24].

In the present work, the roughness that characterizes the surface of the blade is irregular and three-dimensional and it has a spatial statistical distribution that is near Gaussian. This type of roughness has been shown to well represent the surface characteristics of real-life blades subjected to operational metal erosion [27]. The key topographical parameters of the current distribution are shown in Table 3. The root mean square roughness height k_rms of the selected distribution is k_rms = 0.0004 C_ax, corresponding to an equivalent value of Nikuradse sand grain roughness of k_s = 0.002 C_ax (k_s = 5.0k_rms). For a blade with an axial chord of 50 mm, this is equivalent to a surface roughness distribution with k_s = 100 μm (and k_rms = 20 μm). The value of k_s normalized with respect to the local averaged viscous length scale on the suction side of the reference smooth blade, indicated as $ks+$⁠, is shown in Fig. 9. Both the skewness and the excess kurtosis of the distribution are approximately 0. The streamwise effective slope, $ESs=As−1∬As|∂hr(s,z)/∂s|dAs$⁠, and the spanwise effective slope, $ESz=As−1∬As|∂h(s,z)/∂z|dAs$⁠, are 0.16 and 0.17, respectively.

A comparison of the blade surface at the trailing edge between the smooth and the rough blade is shown in Fig. 10.

Validation of the reference smooth blade in terms of surface pressure coefficient, skin-friction coefficient, and wall-normal heat flux are carried out by comparing the BDIM results against the numerical simulations of Pichler et al. [20] and Zhao and Sandberg [28]. Validation of the reference results against experiments is included in the original papers.

Figure 11(a) shows the spanwise and time-averaged pressure coefficient $C¯p$ on the blade surface, exhibiting excellent agreement with the reference. The skin-friction coefficient $C¯f$ is shown in Fig. 11(b), demonstrating that the present method is able to correctly simulate the boundary layer transition location and evolution in the vicinity of the blade trailing edge. The BDIM data show a slight under-prediction of $C¯f$ at x = 0.25 on the suction side (⁠$≈10%$⁠). The agreement of the overall $C¯f$ trend with the reference is, however, excellent. A more detailed discussion of the computation of the skin-friction coefficient $C¯f$ is presented in the next section. Finally, the heat flux $Q¯$ is compared in Fig. 11, again showing good agreement of the overall trend, despite some discrepancies that are slightly larger, especially on the suction side, with respect to the previously discussed quantities.

The instantaneous snapshots of the flow field around the smooth and the rough blades shown in Fig. 12 provide a visual representation of some of the complex physics that characterize HPT flows. Iso-volumes obtained from the Q-criterion offer qualitative information on the turbulent structures present in the boundary layers and in the wake region, highlighting the differences between the two blade configurations. For the smooth blade, on-surface turbulent eddies can only be observed in the trailing edge region, suggesting that the boundary layer behavior is predominantly laminar. This is not the case for the rough blade. The surface roughness acts as a trip and has the effect of promoting a much earlier transition. As a result, the boundary layer on the suction side is predominantly in a fully turbulent state. In Part II, where additional roughness amplitudes and distributions are also presented, we will see that this has direct implications on the loss and on the overall performance of the blade. The instantaneous velocity divergence field collected at mid-span and displayed in Fig. 12 highlights a series of normal shocks on the blade suction side, suggesting that the flow is subjected to a strong acceleration in the HPT blade, reaching transonic conditions. The two flow snapshots suggest that the surface roughness has an impact on the location and the strength of the shock waves. For the rough blade, the shocks are shifted upstream and their intensity seems to reduce with respect to the smooth blade. Further analysis is necessary to fully quantify and characterize this behavior. It is interesting to observe that the normal shocks interact not only with the boundary layers but also with the wakes of the other blades present in the linear cascade, once again an indication of the complex dynamics observed in HPT flows.

A comparison of the skin-friction coefficient averaged in time and in the spanwise direction can be found in Fig. 13. For the spanwise averaging (see Sec. 2.2.5), in the calculation of $Cf*$⁠, the pressure term has been replaced by $pi−p¯bin$⁠, where $p¯bin$ is the pressure averaged in each bin. This is necessary to reduce the noise present in the curve of $Cf*$ due to the fact that p is several orders of magnitude larger than τ_w. This is justified by the fact that what contributes to the pressure drag is not the absolute value of p, but the pressure variations in the tangential direction. From Fig. 13, we can see that in the case of the rough blade, $C¯f$ only accounts for about 50% of the overall skin-friction coefficient $C¯f*$⁠. Roughness has the effect of increasing the overall skin-friction coefficient on both the suction and the pressure side, promoting the transition of the boundary layer from laminar to turbulent, as already inferred from Fig. 12. We refer to Part II for a more detailed discussion on the effect of roughness on skin friction and heat flux.

A framework for high-fidelity simulation of transitional and turbulent flows over multi-scale roughness has been presented. The backbone of the framework is a three-dimensional, second-order accurate BDIM that allows us to simulate complex three-dimensional geometries on non-conforming structured grids, hence retaining high computational efficiency. The general formulation of the BDIM has been presented, with a particular focus on the procedure used to map the surface of the solid body, discretized with triangular elements, onto the computational grid. The effect of a solid boundary on the flow is obtained by directly enforcing local velocity and temperature values to those grid points located in the vicinity of the boundary. A BDIM data-capturing method used for on-surface data collection and analysis has also been introduced. It has been shown that collecting on-surface data at the element nodes generally yields more accurate results, due to the presence of larger gradients in these regions.

The BDIM and the surface data-capturing have been validated on fully developed turbulent channel flow over sinusoidal egg-carton roughness in a minimal-span domain simulated by means of DNS. The boundary layer profiles have been compared to reference data present in the literature, showing good agreement. The integral streamwise force and the total heat flux, computed numerically using the BDIM data-capturing method, have been validated using momentum and energy conservation principles. The rough channel case also served to discuss the BDIM grid resolution requirements.

The BDIM has then been employed to investigate the effect of irregular homogeneous surface roughness on the performance of a high-pressure turbine blade at engine-relevant conditions using DNS. The resolution required to accurately resolve the interaction between the flow and the surface roughness has been achieved by an overset three-block grid, in which two overset o-grids with different resolutions—a refined inner layer for the near-wall region and a coarser outer layer for the external boundary layer—have been employed to resolve the blade, while a Cartesian background grid has been used to resolve the incident turbulence and the far wake. Flow visualizations have shown that the surface roughness promotes the boundary layer transitions on the suction side, shifting the transition region upstream with respect to a reference smooth blade. This has also been confirmed by comparing the time-averaged and the spanwise averaged skin-friction coefficient. We have seen that when comparing on-surface data between rough and smooth blades, the results should refer to the geometry of the smooth blade: in particular, surface quantities should be normalized by the same reference area and the same tangential direction should be considered when comparing the skin-friction coefficient. As a consequence, in the case of a rough blade, the skin-friction coefficient does not only depend on the wall shear stress but also on the tangential component of the pressure force, which accounts for about $50%$ of the total coefficient.

Support from the ARC is acknowledged. This research used data generated on the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR2272.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.