Simulation strategies for the Food and Drug Administration nozzle using Nek5000

The use of computational methods for biomedical applications is currently a rapidly emerging area of research. Gaining insight into the complex flow of blood facilitates, e.g., the assessment of the performance of medical devices. This, in turn, allows experts to predict biological responses of interest such as thrombosis and hemolysis.¹ Several studies have shown that large viscous shear stresses can cause detrimental effects on blood constituents during enough exposure time,^2,3 an effect that is enhanced for turbulent flows.^4,5 Under normal blood-flow conditions, shear-sensitive endothelial cells on the vessel wall play an essential role in sustaining vascular integrity and functionality.⁶ These physiological responses can, however, be impaired as the near-wall flow attains more disturbed, turbulent-like characteristics, which increase the susceptibility of vascular diseases.⁷ Here, computational methods have the potential to obtain detailed characterization of these flow regimes, beyond state-of-the-art measurement techniques, which could be used to better understand these currently not well-understood pathological mechanisms. Recently, studies have been conducted in blood filters,⁸ endovascular stents,⁹ prosthetic heart valves,¹⁰ and ventricular assist devices.¹¹ Nevertheless, blood flows possess great complexity, and better characterization is still required to provide accurate and robust predictions. Therefore, further validation and research on modeling techniques is essential to enhance the reliability of computational fluid dynamics (CFD) so that it could be used for safety and regulatory purposes in the biomedical field.

The United States (US) Food and Drug Administration (FDA) initiated the “Critical Path Initiative” program with the aim of developing a robust and standardized CFD methodology for biomedical applications.¹² To that end, a nozzle benchmark case was proposed that resembles an idealized medical device consisting of flow features such as recirculating, expanding, or contracting characteristics, as shown in Fig. 1. In 2011, an inter-laboratory study including 28 independent groups of researchers performed a series of simulations using Reynolds-averaged Navier–Stokes (RANS) simulations to assess the accuracy and limitations of turbulence models in the prediction of flow dynamics for the FDA nozzle case.¹³ Hence, five Reynolds-number conditions in the range of Re_t = 500–6500 at the throat were tested involving laminar, transitional, and turbulent flow regimes. In addition, three groups provided experimental data obtained from particle image velocimetry (PIV) measurements for validation purposes.¹⁴ Nonetheless, evident sources of uncertainties were acknowledged in the experiments related to inaccuracies in fluid properties measurement, PIV algorithm errors, and variability in the inlet boundary condition. The computational results showed considerable discrepancies among the different turbulence models, with no approach presenting a clearly outstanding performance compared to the experimental data. In the low Reynolds-number case of Re_t = 500, most of the laminar inflow-cases yielded satisfactory results, and no flow transition was observed, neither in the experiments nor in the simulations. This observation evidenced the laminar nature of the flow at these low Reynolds conditions, thus reducing the existing uncertainty in both the experimental measurements and simulations. At the moderate Reynolds number (Re_t = 2000), centerline velocities were generally underestimated by the turbulence models in the inlet and nozzle regions, probably due to the laminar and transitional nature of the flow in these zones. At higher Reynolds numbers (Re_t ≥ 3500), the turbulence models failed to accurately predict the velocities and shear stresses, especially the wall-shear stress (WSS), in the recirculation region downstream of the expansion. These observations called for future work to focus on better simulation approaches for transitional effects in the flow, including more refined transition models.

The accuracy of unsteady RANS (URANS) models together with hybrid RANS and large-eddy simulation (LES) methods has also been investigated.¹⁵ A variable performance was observed in the case of Re_t = 3500, which presents laminar inflow conditions and turbulent flow after the expansion. The results showed that URANS models underpredicted axial velocities, especially in the throat and expansion region. Hybrid RANS/LES (HRL) models presented a mixed behavior, and dynamic HRL provided better results than URANS, whereas the commonly used detached-eddy simulations (DES) completely failed to predict turbulence and behaved as a laminar model. With the aim of obtaining a better representation of the physics in the transitional and turbulent regimes, a number of investigations based on LES have been performed.^16–18 Janiga¹⁶ focused on the fully turbulent case at Re_t = 6500 with transitional inlet conditions. Validation of results showed generally good agreement between the temporally averaged flow velocities and the experimental measurements, with a slight underprediction of the values at the throat. Zmijanovic et al.¹⁸ studied the case of Re_t = 3500 and reported the location of the jet breakdown to be very sensitive to variations in the numerical parameters (i.e., grid size, time step, and time advancement scheme). In the same way, high dependency on the numerical scheme was also observed in low-order direct numerical simulations (DNSs) for transitional cases at Re_t = 2000.¹⁹ 

A major part of the previous studies relied on laminar (parabolic) inflow conditions,^{15,17,19–21} even in certain cases where the flow was known to be at least transitional, as in the case of Re_t = 6500 studied by Janiga.¹⁶ Nevertheless, studies performed by Zmijanovic et al.¹⁸ revealed that turbulent injection had a strong impact on the fidelity of turbulent and transitional flow features. In this way, greater robustness was achieved by introducing small perturbation in the inflow, leading to a more accurate and stable prediction of the jet breakdown location. The effect of injecting noise perturbations also led to better agreement with the PIV results in the study from Bergersen et al. ²¹ Fehn et al.²² employed a precursor simulation to prescribe the inflow field, which yielded overall satisfactory results in the turbulent cases of Re_t = 2000–6500. However, the breakdown location was slightly later than the PIV data, which was attributed to the presence of larger flow fluctuations in the experiment. This variation increased when imposing a laminar inflow profile. In the transitional case of Re_t = 2000, the jet was observed to remain stable in the simulations, in contrast to the PIV measurements, which presented high uncertainties in the mass flow rate.

In terms of spatial discretization techniques, finite-element methods (FEMs),^17,19,21 finite-volume methods,^15,16,18 finite-difference methods with immersed boundary conditions,²⁰ and discontinuous Galerkin spectral element²² methods have been employed. In this regard, there is room for investigating the performance of high-order methods such as the spectral-element discretization. The spectral-element method (SEM) is based on the weighted residual principles of finite elements, relying on high-order basis functions employed in spectral methods.²³ Consequently, it offers geometrical flexibility as traditional FEM and the high accuracy provided by the high-degree spectral functions. The exponential convergence of the latter allows using a fewer number of elements in SEM than the traditional low-order FEM approach.²⁴ In the same way, a comparative study of the finite-volume solver OpenFOAM and the SEM solver Nek5000 showed the latter to be more computationally efficient for high-fidelity simulations.²⁵ 

The present work consists of an analysis of the performance of high-order SEM capabilities in the flow characterization within the FDA nozzle. High-order DNSs have the potential to provide complete and comprehensive benchmark data with a low degree of uncertainty as, e.g., compared to existing PIV measurements. Four different flow regimes including laminar, transitional, and turbulent conditions were studied (i.e., Re_t = 500, 2000, 3500, and 5000) using a DNS approach in order to achieve maximum accuracy and thus fidelity in the results. A parabolic velocity profile and the synthetic-eddy method (SyEM) were used in each case to study the impact of different inflow conditions on the flow. A three-dimensional mesh with a hemisphere configuration was developed, and the accuracy of the results was evaluated against available PIV experimental data.

The paper is organized as follows: in Sec. II, the numerical setup, including the meshing strategy, is introduced; in Sec. III, we show the results of the study, including instantaneous flow fields and statistics; the results are thoroughly discussed in Sec. IV; finally, the conclusions of the work and the outlook are presented in Sec. V.

The numerical simulations in this work have been performed using the Nek5000 CFD solver, an open-source code developed by Fischer et al.,²⁶ which is based on the high-order spectral-element method. It features relevant advantages compared to low-order methods (such as minimal dissipation and dispersion errors), as well as excellent parallel performance. Furthermore, the code is characterized by its high computational efficiency as it employs a tensor–product formulation, thus reducing the computational cost and memory usage. The code is mainly written in Fortran 77, with specific sections written in C, such as the communication library. The spatial discretization is based on the spectral-element method developed by Patera,²³ which combines the geometrical flexibility of the finite element approach and the high accuracy of spectral methods. In this way, it offers a high level of accuracy with fewer grid points compared to traditional low-order FEM methods. Hence, the SEM exhibits an advantageous computational efficiency for high-fidelity simulations. In the same way as in the FEM, the SEM formulation is based on a variational principle. Consequently, the global domain Ω is divided into a finite number of non-overlapping elements Ω^e on which the solution is approximated by a polynomial expansion following the Galerkin approach, and global continuity at the element boundaries is enforced. In Nek5000, the basis functions are determined by the Lagrangian interpolation based on Legendre polynomials P_N of order N. Furthermore, the P_N–P_N−2 formulation proposed by Maday and Patera²⁷ is employed, where the basis functions associated with the pressure field are two polynomial orders below N, which is the order used for the velocity field, effectively eliminating spurious pressure modes. It should be noted that one of the greatest advantages of the SEM is the exponential decay of the error as the order of the polynomial increases. These Lagrangian basis functions exhibit local support on each element and are defined in the interval [−1, 1]³ for a three-dimensional domain; thus coordinate transformation from the actual to a reference element $Ω^≔[−1,1]3$ is needed. The selected interpolation points within each element are the Gauss–Lobatto–Legendre (GLL) points (⁠$ξj=0N$⁠), which are the roots of the polynomial (1 − ξ²$PN′$(ξ)), where ξ ∈ [−1, 1]. The GLL points are characterized by a non-uniform distribution within the element, exhibiting a higher density of grid points near the element boundaries. One remarkable characteristic of the Lagrange interpolants ϕ is that they satisfy ϕ_j(x_i) = δ_i,j (where δ_ij is the Kronecker delta function), and therefore, each basis function differs from zero only at one specific GLL location. The orthogonal nature of the Lagrange basis functions leads to a (near-) diagonal mass matrix, which greatly benefits the computational performance.

For the simulations, the DNS approach has been chosen, which means that all relevant flow scales are resolved on the numerical grid, and, consequently, no turbulence model is needed. The four Reynolds numbers Re_t = U_bD_t/ν shown in Table I are considered, where D_t is the throat diameter, ν is the kinematic viscosity, and U_b is the bulk velocity at the throat. The obtained results were non-dimensionalized, considering these reference quantities, i.e., U_b and D_t at the throat. Furthermore, 1% of the energy of the highest frequency mode is chosen to be explicitly removed by a high-pass filter based on previous experience with the numerical stability of Nek5000.^28,29 The fundamental principle of DNSs is based on resolving the unsteady Navier–Stokes equations for a wide range of existing spatial and temporal scales. Thus, it represents the most powerful numerical approach that provides the high accuracy and precise characterization of the flow field. The main drawback of this method is its computational cost. Thus, the grid resolution was evaluated in terms of the viscous length ℓ^* = ν/u_τ, where the friction velocity is $uτ=τw/ρ$⁠, which is expressed in terms of the wall-shear stress $τw$ and the fluid density ρ. An inner-scaled wall distance of $rmax+≤1$ (where + denotes scaling with ℓ^*) to the first grid point was maintained for all the cases. In addition, $Δzmax+≤15$ in the streamwise direction, $RΔθmax+≤12$ in the azimuthal direction, and $Δrmax+≤8$ in the radial direction were achieved over the whole domain. In the same way, small time steps are needed in order to accurately retrieve the flow statistics and thus to properly compute the fluctuating terms. In the present work, a second-order backward difference scheme (BDF2), coupled to a second-order extrapolation scheme (EXT2), was used to evaluate the temporal discretization of the solution. The choice of the time step Δt was made according to a Courant–Friedrichs–Lewy (CFL) condition so that the Courant number remains below 0.5. Furthermore, the natural stress-free condition was used at the outlet in all the cases under consideration. Previous DNSs with Nek5000 were taken as the reference in order to set the time resolution criteria.³⁰ 

The flow was simulated for 400 time units in the cases of a turbulent inflow, which as described below was obtained by means of the synthetic-eddy method (SyEM). The breakdown location quickly stabilized, and no major variations were observed within this period. In the cases of a parabolic inflow, the breakdown location shifted continuously toward the outlet, as also reported in previous studies.²¹ Therefore, results were extracted for 300 time units for comparison with the SyEM cases, which was sufficient to obtain converged mean profiles and fluctuations.

In the present study, the effect of the inflow conditions is evaluated through a comparison with the PIV experimental database generated by Hariharan et al.¹⁴ Their error estimation was performed considering a confidence interval of 95%, an approach previously employed by Stewart et al.¹³ Furthermore, the parabolic velocity profile and turbulent inflows are tested in the present study, where the latter was applied in the transitional and turbulent cases of Re_t = 2000–5000. The synthetic-eddy method (SyEM)³¹ used in this work to generate turbulent inflow conditions is based on the divergence-free approach of Poletto et al.,³² implemented in Nek5000.³³ In this method, a series of synthetic eddies are generated within a box around the inflow plane and convected by the mean velocity. As the eddies exit the box, new coherent structures are continuously regenerated. To apply this method, the input of the mean velocity U(r), the desired turbulent kinetic energy k(r) profile, and the dissipation rate ϵ(r) of the eddies as functions of the radial position are required. The turbulence statistics used for the SyEM are based on the results obtained from a previous DNS study,³⁴ considering a straight pipe with Re_τ = Ru_τ/ν = 180. The method was verified for Nek5000 by Hufnagel et al.,³³ and the minimum inlet length criteria determined in that study are followed.

The geometrical model used for the present simulations replicates the configuration proposed by the FDA, consisting of a pipe with a nozzle and an adjoined sudden expansion. The geometry dimensions are based on the experimental model and are normalized with respect to the throat diameter D_t, as shown in Fig. 2. Note that all corners are sharp. In the laminar case at Re_t = 500 and the transitional case at Re_t = 2000 with a parabolic inlet, a short configuration consisting of 84 256 elements was used with 25D_t after the expansion, as seen in Fig. 2. In the turbulent cases at Re_t = 3500 and Re_t = 5000, as well as in the case of Re_t = 2000, with SyEM, the outlet channel was extended to twice its original length, up to 50D_t, due to the presence of flow breakdown into the turbulent regime. The latter configuration consisted of 162 976 elements. The mesh configuration and polynomial order used in each simulation are summarized in Table I.

The simulations were performed using a grid consisting of hexahedral elements with conformal nodes. The meshing strategy was based on the hemisphere mesh previously tested in studies for jets in cross-flow.³⁵ Hence, a three-dimensional mesh, as shown in Fig. 3, was created by means of the open-source Gmsh meshing tool.³⁶ In this way, high resolution is achieved locally downstream of the jet expansion due to the non-uniformity of this configuration. Furthermore, the grid resolution has been studied by testing two polynomial orders of N = 5 and 7 in the cases of Re_t = 2000 and Re_t = 5000 with a parabolic inflow. The results obtained in the case of Re_t = 2000 matched for both polynomial orders, and therefore, only the results of N = 5 have been included. In the case of Re_t = 5000, higher refinement was not observed to provide a better prediction of results. Due to this reason and the high computational cost required for simulating a high polynomial order of N = 7, the rest of the simulations were performed considering N = 5.

Instantaneous streamwise velocity contours are displayed in Figs. 4–7 to provide a first overview visualization of the flow development, the effect of inflow conditions, and grid resolution. Note that the color scales have been adjusted for visualization purposes. As observed in Fig. 4, steady laminar flow is predicted for Re_t = 500 throughout the whole domain. This outcome is not surprising as the bulk Reynolds numbers in the throat as well as at the inlet are much lower than the Reynolds number established for sustained pipe turbulence (Re = 2040, see Avila et al.³⁷). In the case of Re_t = 2000, the flow was observed to be transitional in the experiments.¹⁴ Nevertheless, the cases carried out with a parabolic inflow did not show a jet breakdown in the short domain, as observed in Fig. 5(a), where the velocity remains stable and steady after the expansion. The reason for not capturing the transitional behavior of the flow can be related to the subcritical nature of the flow in the pipe sections, for which the background noise level in the simulation is presumably too low [Fig. 10(a), parabolic inflow] to trigger any nonlinear interactions. However, by employing a longer domain and increasing the fluctuation levels via the SyEM [Fig. 10(a), SyEM], the transitional nature of the flow was captured, as shown in Fig. 5(b). Note, however, that the bulk Reynolds number in the wider pipe is only 667, meaning that the turbulent patches further downstream are eventually dying out again.

Finally, turbulent flow after the expansion was observed in the experiments for both cases of Re_t = 3500 and Re_t = 5000, as also seen in Figs. 6 and 7. It should be noted that the simulations performed with the parabolic inflow predicted the transition to turbulence to take place further downstream compared to SyEM cases. In the same way, increasing the polynomial order in the case of Re_t = 5000 shifted the breakdown location slightly toward the outlet. The fluctuating behavior imposed by the SyEM at the inlet can be observed in Fig. 8. The local bulk Reynolds number, based on the pipe diameter and bulk velocity, is Re = 1667, which is below the threshold Reynolds number for sustained pipe turbulence. However, since new fluctuations are constantly fed into the domain through the SyEM and the distance to the accelerating section of the nozzle is not large, quasi-stationary turbulence can be expected to reach the nozzle.

To summarize these first visualizations, one can say that the experimental observations are essentially recovered using the present simulations, i.e., the relevance of the inflow turbulence is clearly established. In addition, a slight sensitivity to resolution could also be observed for the highest-Re case.

The results obtained for the mean axial velocity along the centerline of the FDA nozzle are shown in Fig. 9 together with the experimental data by Hariharan et al.¹⁴ As one can observe, for all Reynolds numbers, the flow rapidly accelerates as it enters the nozzle and it reaches the maximum velocity at the throat. In general, good agreement with PIV data is achieved at the inlet, nozzle (z = −15.67 to −10), and throat regions (z = −10 to z = 0) with the simulations based on a parabolic inflow, whereas the SyEM-cases underpredicted the centerline velocity near the inlet. Note that a fully laminar flow reaches a non-dimensional velocity of U_z = 2/3² = 0.22 at the inlet, whereas the turbulent inflow is about 0.17 due to the fuller profile. The centerline velocities for the turbulent cases slightly increase, which is most probably due to the low Reynolds numbers as discussed.

Beyond the sudden expansion located at z = 0, the flow starts to decelerate due to the sudden increase in the cross-flow area. The results obtained for Re_t = 500 match those of the experimental data very well, as seen in Fig. 9(a), throughout the whole geometry. Differences between experiments and simulations appear after the expansion when predicting this velocity drop in the transitional and turbulent cases. In general, the experiments show an abrupt decrease in velocity, with less intermittency for higher Reynolds conditions, due to the earlier turbulence breakdown. In the simulated case of Re_t = 2000, the centerline velocity was observed to be quite stable, starting to decrease only after 20D_t in the case of the SyEM as opposed to the sudden drop which took place after 12.5 D_t in the experiments. The turbulent cases at Re_t = 3500 and 5000 also showed this abrupt reduction further downstream compared to the experiments, especially with the parabolic inflow (i.e., the ones with a lower fluctuation level). The predicted breakdown locations with a turbulent inflow (SyEM) shift toward the data from PIV with increasing Reynolds number, as can be seen when comparing the results from Figs. 9(c) and 9(d). Also shown in panel (d) is the comparison of the simulation with a higher polynomial order N = 7. For Re_t = 2000, the two lines were practically indistinguishable throughout the domain.

The variance of the fluctuating axial velocity $uz′$_,rms along the centerline is displayed in Fig. 10 for all the cases. Re_t = 500 remained steady, and the fluctuations are identically zero in the simulations. In the experiment, there are small fluctuations at the throat, slightly lower than those at Re_t = 2000, although the jet breaks down to turbulence in the latter as discussed. For all other cases, for the upstream of the nozzle z < −15.67, the simulations with a laminar inflow show, as expected, no fluctuations (see the inserts in Fig. 10) and some minor fluctuations in the converging and throat sections (up to z = 0), respectively. It is interesting to note that for the SyEM inflow, the fluctuation levels upstream of the nozzle z < −15.67 decrease for Re_t = 2000, whereas they are more or less constant for higher Re_t; this is again in agreement with the low local Reynolds number in the upstream pipe discussed. Once the converging section reached z = −15.67, the PIV measurements show an unexpected increase in fluctuations, which are most probably due to experimental imperfections in the nozzle section. Conversely, the simulations show a rapid decrease in fluctuations between z = −15.67 and −10, consistent with flow acceleration. Thus, one can observe that the flow state reaching the sudden expansion at z = 0 is noticeably different for various cases, which is mostly independent of the Reynolds number. Beyond the sudden expansion, a peak is observed in all the cases due to the flow breakdown. As discussed in conjunction with the centerline velocities, the breakdown position is different for the various cases.

Figures 11–13 show the mean axial velocity profile along the radial axis at six different downstream locations for Re_t = 500, 2000, and 5000. The data obtained from the simulations performed at the lowest Reynolds number of Re_t = 500 are in very good agreement with the experimental data at all stations (Fig. 11). For higher Re_t, velocity profiles similar to the experiments were obtained in the cases of a parabolic inlet, whereas clear differences were observed in the SyEM cases at the inlet due to the modified mean profile and the imposed fluctuations. Nonetheless, both inflow conditions yielded good agreement that generally matched the experiments in the throat region (at z = −12), where the profile became fuller as a result of acceleration. In the same way, good agreement is observed right after the sudden expansion where the weak recirculation zone was properly captured and the negative velocity profile was predicted between the wall and the jet areas. Further downstream, great discrepancies in the profiles were observed at 20 D_t in the case of Re_t = 2000 between experiments and simulations (Fig. 12) because the simulations remained laminar with a more pronounced backflow region, whereas the experimental flow transitioned to turbulence. Similar findings were observed at Re_t = 3500 and the parabolic-inflow case.

In this section, both the wall-shear stress along the axial coordinate and the viscous stresses in various cross sections are evaluated. It should be noted that the absolute value of these terms was computed for comparison with the available PIV experimental measurements. Figure 14 shows the evolution of the WSS along the streamwise direction. First, the comparably wide error margins in the experimental results must be remarked. Due to issues with the measuring procedure (i.e., PIV for the analyzed data), near-wall mean statistics yield considerable uncertainty, particularly in the higher-Reynolds-number pipe section. Nevertheless, for all Re_t, the best agreement can be observed in the region containing the (upstream) nozzle throat. As mentioned previously, the simulations with a SyEM inflow tend to predict a higher shear stress due to their fuller non-parabolic mean velocity profiles. Along the throat region (before z = 0), all simulations capture the plateau in the smaller pipe section, agreeing well with the mean experimental results despite their considerable error bars. The maximum shear stress is obtained exactly at the beginning of the throat, and it reaches values approximately four times higher than those in the following straight section. This peak was not measured in the experiments, but given the fine mesh in this region, it is likely to be captured with high fidelity in the simulations. Overall, it can be concluded that the WSS is accurately predicted for z < 0 in simulations and experiments alike.

The insets in the various panels of Fig. 14 show a close-up of the evolution of the wall-shear stress. At the lowest Reynolds number, which was fully steady in both experiments and simulations, even the low-shear-stress values were predicted with good accuracy, including the recovery just after the expansion. In the proximity of the expansion, a peak in WSS magnitude is captured in both the experiments and the simulations. This increase in WSS magnitude is observed to take place slightly after the flow breakdown, and thus, the simulations predicted this location further downstream compared to PIV measurements.

As previously mentioned, the turbulent inflow with the SyEM provided better results in this region than the parabolic-inflow cases since the flow was less stable after the nozzle. In addition, for the highest Reynolds number, a slight shift in the reattachment point could be observed toward the outlet when increasing the polynomial order. This shift, however, was much smaller than the difference between the laminar and turbulent inflow.

The interior, viscous stresses were computed based on the strain rate tensor, properly averaged in the azimuthal direction and over time. The resulting profiles along radial cuts are displayed in Figs. 15–17 for Re_t = 500, 2000, and 5000, respectively. At the lowest Reynolds number, as shown in Fig. 15, a completely linear profile is recovered close to the inflow, confirming a parabolic velocity profile for both experiments and simulations. The accelerated flow through the nozzle increases the stresses, particularly at the center region, which gradually starts to show a linear profile again when approaching z = 0, albeit with a value 27 times larger. After the expansion, the highest shear rates are obtained at r = 0.5, along the strong shear layer exiting the small pipe, while in the outer region, the gradients initially are much smaller. The emergence of the region with a reversed flow can clearly be identified by the kink in the shear-stress profiles at r ≈ 1 and the subsequent negative stress at the wall (at r = 1.5). Overall, there is excellent agreement with the experimental data even though the provided error bars are very wide in certain regions.

Similar to the results obtained for the velocity profiles, a notable mismatch with respect to the experiments was observed at the inlet in the cases with the SyEM (Figs. 15–17). The linear behavior of the experiments is correctly predicted by the laminar simulations due to the parabolic velocity profile at this region. For locations downstream, the profiles become increasingly complex as the flow is accelerated but are captured correctly taking into account the large errors of the mean experimental profiles close to the centerline. After the sudden expansion, the experimental shear-stress profiles show a near-Gaussian curve with a distinctive peak corresponding to the interface between the jet and the recirculation area, which are accurately replicated by the laminar simulation of Re_t = 500. For simulations with the SyEM, due to the injection of random spots for turbulence generation, the statistical profile was not parabolic at the inlet, and hence, the resulting shear-stress profile was not linear. After the expansion, differences between the simulation and the experiment arose due to the late prediction of the flow breakdown, especially in the cases with parabolic-inlet conditions.

The most relevant challenge when simulating the FDA configuration is to obtain an accurate prediction of the transition location in the jet, as reported in previous studies documenting significant difficulties when reproducing the experimental results.^13,18,19,21 The main observation from these studies is the high sensitivity of the jet-transition point to the inflow conditions and consequently the amplitude and nature of fluctuations present in the upstream flow. Therefore, injecting noise perturbations into the inflow yielded a better prediction of the flow breakdown in previous studies,^18,21 in a similar way to the use of a precursor simulation.²² The PIV measurements were carried out in the presence of experimental noise generated at the throat, which affected the downstream evolution of the flow. This issue can be clearly seen in the experimental velocity fluctuations shown in Fig. 10. Moreover, significant uncertainties in the measurements were observed in the case of Re_t = 2000, which amount to nearly 10% of the flow rate.¹³ 

Regarding the results from the present study, in Fig. 18, we show the jet-transition location z_j,tr for all the simulations compared with the reported PIV results. Here, we define the transition location as a streamwise position where the axial velocity is starting to decrease to the turbulent level. We first note that z_j,tr decreases both in the experiments and the simulations for increasing Re_t, i.e., the transition in the jet takes place earlier. This is, of course, expected as the lower viscosity makes the flow prone to instabilities. On the other hand, the jet-transition location is observed further downstream in the simulations than in the experiments, at the three Reynolds numbers under study. In particular, the values of z_j,tr are higher for the parabolic-inflow cases than for the SyEM simulations. The reason for these discrepancies is the presence of higher fluctuations at the throat in the experiments, as previously shown in Fig. 10. The differences in z_j,tr are also noticeable in the streamwise evolution of centerline velocity (Fig. 9) and the shapes of the mean velocity profiles (Figs. 12 and 13) at various streamwise locations. It should also be noted that the jet-transition location was predicted further downstream when increasing the number of grid points in the case at Re_t = 5000 when using the parabolic inflow. Zmijanovic et al.¹⁸ obtained an inconsistent trend in the results when refining the mesh, whereas Bergersen et al.²¹ reported a shift in the jet-transition point toward the outlet. The latter result is consistent with our observations and probably also what one would expect from DNSs: higher resolution will prevent the flow from breaking down prematurely (see, e.g., the LES studies of the transition in a channel flow³⁸). It is interesting to note that the simulations performed with the SyEM provide values of z_j,tr closer to those in the experiment than in the cases where a parabolic inflow was considered. In particular, the jet-transition locations obtained in the simulations are progressively closer to the ones in the experiment as Re_t increases, an effect that is more pronounced in the SyEM cases. Note that the present numerical results show the high sensitivity of the jet to the conditions at the throat, which may significantly affect the transition location.

Regarding the upstream behavior of the contraction, the parabolic-inflow cases exhibit better agreement with the experiments than the SyEM cases near the inlet. The reason for this is that all the experimental cases have laminar conditions in this region, whereas the SyEM introduces fluctuations on top of a turbulent mean flow at the inlet. In the same way, previous studies revealed that turbulence models fail to predict the flow behavior at the inlet, whereas laminar simulations yielded satisfactory results.¹³ Fully laminar conditions were observed in the experimental results for the case of Re_t = 500. Simulation results at this condition also exhibited the same laminar behavior throughout the FDA pipe, and good agreement in the results was obtained. As clearly observed in Fig. 9(a), the flow progressively decelerated after the expansion with no breakdown. Previous studies with the DNSs,¹⁹ LES,¹⁷ and laminar¹³ simulations with the same inflow condition yielded similar results, a fact that supports the laminar observations. Therefore, no synthetic-eddy method was needed in this case for capturing any turbulent behavior. This was expected at such low Reynolds number since it lies far below the critical Reynolds-number condition for a transitional flow.

In this study, we conduct DNSs of the FDA nozzle configuration at Reynolds numbers of Re_t = 500, 2000, 3500, and 5000 using the high-order spectral-element code Nek5000. The effect of laminar (parabolic) and turbulent (generated via the SyEM algorithm) inflow conditions on flow development was analyzed and compared to PIV experimental data taken as the reference. At Re_t = 500, the flow was fully laminar along the pipe, and good agreement with the experimental data was achieved. For higher Reynolds numbers, our results clearly show that the jet-transition location is very sensitive to the inflow conditions. Although the Re_t = 2000 case was observed to be transitional in the experiments, this behavior was only captured in the SyEM case with a long domain. In the turbulent cases at Re_t = 3500 and Re_t = 5000, the jet transition was captured in the simulations performed with both inflow conditions. Here, the parabolic inflow provided predictions closer to those in the experiment near the inlet due to the fact that a laminar inflow was also employed in the latter. On the other hand, the SyEM predicted the transition of the jet downstream from that in the experiment, although in better agreement than the ones obtained in the laminar-inflow cases. This was due to the fluctuations present at the throat in the experiment, which was obviously not considered with the smooth pipe assumption in CFD predictions. Furthermore, it should be noted that excellent agreement between simulations and experiments was observed in the nozzle region for all the cases and both inflow conditions. In the same way, all the simulations properly captured the recirculation zone and the low negative velocities after the expansion. This essentially means that the accelerating parts of the flow are insensitive to the inflow. Therefore, it is noteworthy that the low-amplitude fluctuations that pass through the nozzle then lead to the large sensitivity downstream of the nozzle.

Overall, the relevance of the present study lies in being the first one to tackle the FDA nozzle problem using high-order spectral-element methods and the SyEM to generate inflow conditions. Our results also highlight the extreme sensitivity of this test case when it comes to producing an experimental realization of the flow, particularly regarding the disturbances produced at the throat. This supports the use of high-fidelity numerical data as a benchmark for the FDA nozzle case and sensitive flow cases alike. With a DNS-based reference model, boundary condition uncertainties can be substantially relaxed while also allowing for much more in-depth flow analysis than state-of-the-art experimental techniques. Future studies should be targeted at analyzing the influence of higher Reynolds numbers further, which, of course, also increases the overall computational cost, in addition to quantifying the influence of secondary effects of, e.g., swirl in the inflow and potentially insufficiently developed flow conditions and rheology. Once a robust simulation setup is available, more refined analysis tools, such as modal decompositions and reduced-order modeling of the flow at hand, will be possible.

The simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at HPC2N (Umeå) and NSC (Linköping).