Sensitivity of scale resolving aft-body flow simulations to numerical model parameter variations

The sensitivity of hybrid RANS-LES methods like Improved Delayed Detached Eddy Simulation (IDDES) to numerical model parameter variations related to generic space launch vehicle aft-body flows is investigated. In particular, the changes resulting from the choice of the time-step size, the turbulence model, the fluid modelling, the circumferential grid resolution, the filter length definition, and the data collection period is considered. The results are also compared to experimental and numerical data taken from the available literature. The sensitivity to the time-step size and the turbulence model is minuscule with respect to the obtained mean flow field, wall pressure distributions, azimuthal modes, and wall pressure frequency spectra. However, circumferential resolution, fluid model, and filter length definition affect the solution to a higher extent. Buffeting spectra are very sensitive to the data collection period.


Introduction
The numerical simulation of space launch vehicle aerodynamics features several interesting phenomena. The aft-body region, i.e., the tail section of the launch vehicle, is of particular interest due to the significant mechanical loads that can impact the delicate nozzle structure. Specifically, the so-called buffeting loads can exhibit forces on the nozzle with dominant frequencies [22], which under certain conditions can determine the critical design load for the structural integrity of the nozzle. Since this phenomenon is highly unsteady in nature and is driven by the turbulent shear layer that separates from the base of the launch vehicle main body, it cannot be captured adequately with Reynolds-Averaged Navier-Stokes (RANS) methods. Instead, time and scale resolving simulations like Large Eddy Simulation (LES) are necessary to compute and analyse the expected loads (e.g., Refs. [22,24]).
On the other hand, the high grid resolution requirements connected to LES prohibit to resolve the whole flow around a launch vehicle even at remotely realistic Reynolds numbers. Alternatively, hybrid RANS-LES approaches are available that treat different regions of the flow by RANS and LES(like) approaches (e.g. Ref. [17]). This allows to reduce the computational effort and achieve Reynolds numbers much closer to flight conditions. However, even with more economical approaches, the computational resources required clearly depend on the chosen grid resolution. Furthermore, other parameters affect the required computational resources, as well. For example, doubling or halving the used time step size in implicit time-stepping schemes increases or reduces the required resources by up to factors of 2. Similarly, switching to a different turbulence modelling approach (e.g., from a 1-equation to 2-equation model) for Detached Eddy Simulation (DES) type approaches can increase or decrease computational effort, as well. On the other hand, the used models or parameters can not be chosen purely based on the required computational effort, but the sensitivity of the solution to these changes also has to be taken into account to allow a judgement of the validity, e.g., of the predicted loads.
Similar considerations are required with respect to the choice of the fluid modelling approach. When, e.g., launch vehicle flows with realistic exhaust plumes are considered, the 1 3 fluid properties of the plume differ significantly from those in the free stream. Consequently, a model approach with (at least) two different components (and hence mass conservation equations) is required.
In previous publications, [3,22,24] among others, these sensitivities are not explicitly reported. Instead usually only one, supposedly optimal, set of parameters to investigate the flow phenomena is reported. The present work aims to provide further and systematic insight into the issue which parameters most significantly affect the solution-and hence require closer attention when setting up a simulation-and which parameters affect the solution to a lower extent and thus can be chosen with more liberty. Additionally, the investigation allows to determine the respective deviations that can be expected if simulations with different parameter settings are conducted. For this purpose, the sensitivity of scale resolving launch vehicle aerodynamics simulations to different parameter variations is systematically investigated using the same numerical platform.
After an introduction to the used numerical method in Sect. 2, the considered generic test case is presented in Sect. 3. Finally, the results of the simulations are presented in Sect. 4. These are structured by first presenting the simulation results with our chosen reference settings. Subsequently, the results for changed time-step sizes, turbulence model, and fluid modelling are investigated. Finally, the effect of a changed circumferential resolution and changed filter length are shown before conclusions are drawn in Sect. 5.

Numerical methods
The simulations are performed using the DLR computational fluid dynamics (CFD) solver TAU [6] to solve the compressible Navier-Stokes equations. For multiple species N s , these consist of N s + 4 equations and can be written as where underlined quantities represent vector variables and bold type indicates tensor variables. In this equation are the conservative variables with mixture density , species mass fractions Y s , velocities u , and total energy E. The inviscid fluxes F Eu and the viscous fluxes F NS are defined as respectively, and the source terms are Q = s , 0, 0 T ; however, for the presented test case, no reactions are considered and thus the chemical source term s = 0 . Quantities used in Eqs. (2) and (3) are pressure p, total enthalpy H, temperature T, diffusion coefficient D, viscous stress tensor , Reynolds stress tensor , species enthalpy h s , heat transfer coefficient , turbulent kinetic energy k, and the effective diffusion coefficient k of the turbulent kinetic energy k, with subscript tot denoting the combined contributions from the laminar and modelled turbulence quantities. In this formulation, the equations are already Reynolds-/Favre-averaged to obtain the RANS equations (respectively, Favre-filtered in LES mode) with the overbar denoting averaged (filtered) and the tilde denoting Favre-averaged (Favre-filtered) quantities.
To improve readability, overbar and tilde are dropped in the following. The turbulent diffusion and heat flux are modelled by using constant turbulent Prandtl and Schmidt numbers of Pr t = 0.9 and Sc t = 0.7 . Details on the used quantities and the modelling of the turbulence terms can be found in the literature (e.g., Ref. [8]). The Reynolds stress tensor is modelled using an Improved Delayed Detached Eddy Simulation (IDDES) approach based on a k-SST 2-equation turbulence model [19]. In Subsect. 4.2.2, this model is replaced by an IDDES approach based on the 1-equation Spalart-Allmaras model [20]. The IDDES algorithm automatically switches between RANS mode in regions where attached flow and no flow fluctuations are present and LES-like mode in separated flows and in flows where the grid resolution is sufficient to act in a wall modelled LES-like mode. The filter length definition used in the reference settings is according to Chauvet et al. [2] with the vorticity vector . The sensitivity of the solution to this choice is investigated in Subsect. 4.2.5 by replacing it with an improved definition developed by Mockett et al. [13] where l n = ∕| | × r n and r n are vectors from the cell centroid to every corner of the cell.
TAU uses a second-order accurate finite volume approach to spatially discretize the equations on structured, unstructured, or hybrid grids [6]. For the temporal resolution, a second-order implicit dual time-stepping scheme is used for which the inner iterations are converged using a three-stage Runge-Kutta approach. For the reference settings, the time-step is chosen as the ratio between the free stream velocity and the target grid spacing in the focus region which corresponds to x + 2 in Table 2. This is an approximate best practice approach suggested, e.g., by Spalart [21], but the sensitivity to the exact choice of the time-step is investigated in Subsect. 4.2.1. To allow improved accuracy in resolved regions and improved stability in non-resolved (RANS) regions, a hybrid low dissipation low dispersion central scheme is applied that has been extensively tested and validated for hybrid RANS-LES [14]. The skew-symmetric central scheme used for the convective fluxes [9] for one calorically perfect gas was recently extended for multiple species. It can be written as with the superscripts R, L indicate the values at the left and right side of a face, the face normal vector n , the face normal velocity v n , and the formation enthalpy e 0 . The mass fluxes (per unit area) are computed as which is a trivial extension of the equation for one calorically perfect gas presented by Löwe et al. [10]. The computation of the convectively transported energy, however, differs non-trivially. If multiple species are considered, the internal energy cannot simply be computed as the geometric mean between the left and right states. Instead, this is only possible for the sensible energy, that is The formation enthalpy has to be treated separately, such that it matches the transported species mass flux (see Eq. (6)).
Due to the low dissipative nature of the numerical scheme, a fourth-order artificial matrix dissipation term is added for improved stability [14] which was optimized in LES of channel flows [15]. The artificial dissipation is preconditioned to obtain optimally tuned accuracy also in low Mach number regions [23]. The preconditioning scheme was also recently extended for multi-species flow and has been validated using several two-and three-dimensional test cases [5], including, among others, 2D vortex transport [10] and 3D decaying isotropic turbulence [11].

Test case description and setup
The investigated test case [4] represents a generic space launch vehicle that is approximated by a larger cylinder with diameter D = 100 mm for the main body and a smaller cylinder with diameter D second = 0.4D and length L second = 1.2D for the nozzle structure. The coordinate system origin is located on the symmetry axis and x = 0 is located at the base of the larger cylinder and all boundaries are farfield boundary conditions that are non-reflecting inflow condition for acoustic modes and an open-end condition for subsonic outflow [7]. The turbulent profile prescribed as an inflow condition was determined in a precursor RANS simulation, such that the growing turbulent boundary layer reaches the experimentally given thickness of 0.2D at x = −2.45D. The Reynolds number is Re D = 1.2 ⋅ 10 6 and other free stream (velocity, pressure, and Mach number) and plume (stagnation pressure, stagnation temperature, and expansion ratio) conditions used are summarized in Table 1.
The two considered grids consist of hexahedral elements in the region of interest and feature the same in-plane resolution, which is shown in Fig. 1 and a circumferential resolution of 1.88 • and 0.94 • , respectively; other details on the grids can be found in Table 2. Exemplary non-dimensional grid spacings + = √ w are given at a location in the early part of the shear layer ( x∕D = 0.1 , r∕D = 0.5 , location 1) and in the region above the mean reattachment location ( x∕D = 1 , r∕D = 0.3 , location 2), respectively. A first non-dimensional wall distance y + w < 1 is maintained at all walls. The grids feature hexahedral cells in the focus region and use prismatic cells to reduce the resolution towards the boundaries. To ensure sufficient resolution, a grid sensor [16] is employed. A detailed discussion of the overall grid design as well as an intensive analysis of the grid parameters and validation of the used grid sensors can be found in Schumann et al. [18].
Since both external flow and plume are comprised of air, the fluid can be treated as a calorically perfect gas with constant Prandtl number Pr = 0.72 and ratio of specific heats = 1.4 . However, for configurations with realistic plumes, this is no longer the case, and thus, in Subsect. 4.2.3, the present fluid is treated as a mixture of thermally perfect gases with mass fractions of 0.74 N 2 , 0.26 O 2 and trace amounts of N, O, and NO. With the species properties taken from Ref. [1], the resulting gas mixture possesses approximately the same properties of Pr = 0.72 and = 1.4 at free stream conditions. With this approach, the influence of using a more complex fluid model suitable for the treatment of reacting flows can be analysed.
All computations are initialized from an unsteady precursor simulation on a similar grid. Subsequently, before data acquisition starts, the simulations are continued with the respective parameter settings until no statistically transient behaviour can be observed. Unsteady data of the computations are recorded for 38 CTUs, which was determined to be sufficient to obtain a time averaged solution of the flow. Here, 1 CTU is defined as the ratio between free stream velocity and main body diameter. Using Welch's method [25] with three segments and 50% overlap, the collected data allow a spectral analysis with

Results
The mean flow field downstream of the separation location can be found in Fig. 2 where color contours of the mean axial velocity and streamlines based on the mean axial and radial velocity as well as a snapshot of the circumferential vorticity are displayed. From the mean flow field, the recirculation region enclosed by the base, the second cylinder, and the shear layer is clearly visible. In the snapshot, the break up of the turbulent shear layer can be observed leading to the creation of turbulent structures. These impinge on the external nozzle surface near the mean reattachment location at x r ∕D = 1.172 . A part of the structures is transported further downstream and partially interact with the supersonic plume, while others are entrapped in the recirculation region and are transported upstream towards the base. This general behaviour and the same qualitative mean flow field are obtained for all considered investigations and hence are not shown again in the respective subsections.
For a quantitative evaluation, the results obtained with the reference settings are first compared to data from the literature to show a general agreement with the previous investigations. Subsequently, the sensitivity of the obtained solution to parameter changes is investigated. In particular, the sensitivity is investigated in terms of the mean reattachment location, mean and rms wall pressure distribution, wall pressure spectral content, circumferential coherence modes, and buffeting loads.

Comparison to literature data
Experimental pressure measurements on the cylinder and base walls are available by Depres et al. [4] and Meliga and Reijasse [12] and numerical comparison data can be found, e.g., in Weiss et al. [24] and Statnikov et al. [22]. In Fig. 3, the wall pressure distribution obtained with the reference settings is plotted together with one numerical and two experimental data sets (taken from Weiss et al. [24]). Additionally, the distribution obtained when increasing the data recording period by a factor of 2 (i.e., to 76 CTUs, labeled "long") is plotted. The figure shows that the obtained data agree well with the experimental data and is independent of data collection period. A slightly further downstream pressure minimum than in the numerical data of Weiss et al. is visible, which leads to a better agreement with the experimental data for x∕D > 0.6 . However, one should note that certain deviations (e.g. c p ≈ 0.023 at x∕D ≈ 0.72 ) exist between both experimental data sets, indicating experimental uncertainties when comparing to this data. Additionally, differences in the exact free stream conditions (e.g., free stream turbulence or velocity profiles which are not provided) and the wind tunnel geometry might affect the numerical solution, as well. Hence, the sensitivity to parameter changes is evaluated primarily in comparison to our own reference solution and only secondarily to the literature data.
The rms pressure distribution is well captured by independent of data collection period. Particularly, in the center of the recirculation region, experimental and numerical data points match extremely well. For x∕D < 0.2 , the present computation shows a slightly higher amount of pressure fluctuations that match the experimental data better, whereas for x∕D > 0.8 , the rms pressure predicted by Weiss et al. [24] is slightly larger and matches the experimental data points better. Especially in this region, the longer data recording period leads to an increase of the pressure rms values by about 10%, improving the agreement with the experimental and numerical data from literature.
A key feature in the wall pressure spectra that is present in all publications is a dominant peak at Sr D ≈ 0.15 … 0.2 , particularly in the region 0.4 < x∕D < 0.8 . Additionally, all publications show an increase of spectral content for 0.3 < Sr D < 1 towards the end of the nozzle ( x∕D > 1 ). However, the quantitative spectral content at different frequencies differs, sometimes drastically, between publications. For example, in Meliga et al. [12, Fig. 4], a clearly dominant peak at Sr D ≈ 0.2 appears at x∕D = 0.55 , whereas in Depres et al. [4, Fig. 16]), it is not visible at all. This indicates an uncertainty in the amplitudes as well as axial extend and location of the spectral features, and hence, in the following, the spectra are compared mainly between different parameter variations to estimate the sensitivity of the solution, but not directly with the literature data. Exemplary, the scaled premultiplied power spectral density (PSD) [3] of the wall pressure is shown for one experimental [12] and one numerical [22] data set together with our reference solution with shorter and longer data recording period at a position near the center in Fig. 4 and the end of the recirculation region in Fig. 5. The peak at Sr ≈ 0.2 in the spectral content at the center is clearly visible, whereas at the end of the nozzle, a broadband content around Sr ≈ 0.6 can be observed. The longer data collection period does not impact these features significantly.
Weiss et al. [24] analyse the coherence between wall pressure sensors at the same axial location and different circumferential positions to determine dominant pressure modes. For this, the coherence function is computed from the crossspectral density S 12 and auto-spectral densities S 11 and S 22 at two circumferential locations 1 and 2 with = 2 − 1 and same axial position x using [24]   with frequency f. Subsequently, the real part of the coherence function is transformed into modes applying where the index m indicates the mode number. The shape of these pressure modes is shown schematically in Fig. 6. Both numerical [24] and experimental [4,12] investigations indicate that the symmetric mode C r,0 is dominant at Sr D ≈ 0.1 and the antisymmetric mode C r,1 at Sr D ≈ 0.2 , while all other modes exhibit only small amplitudes. This is also the case for the current investigation, and hence, only these first two modes are plotted together with the literature data in Fig. 7 for the axial position x∕D = 0.72 . The dominant peak for C r,1 around Sr D ≈ 0.2 is clearly visible for all evaluated spectra, but the numerical data sets show reduced amplitudes. The peak for C r,0 at Sr D ≈ 0.1 is not well captured with the reference settings computation, but clearly appears with the extended data collection period.
Finally, the buffet loads on the nozzle structure can be analysed. The scaled premultiplied PSD of the forces in the y-and z-direction are displayed in Fig. 8 for the reference solution with shorter and longer data recording period compared to the data obtained by Statnikov et al. [22]. Additionally, the average of the two directions is displayed. The figure shows that for a shorter data collection period, the spectra in the two directions show significantly different peak values, e.g., 0.95 in the y-direction and 0.4 in the z-direction for Sr ≈ 0.2 and 0.3 in the y-direction and 1.15 in the z-direction for Sr D ≈ 0.3 . With a longer data collection period, the resemblance between the spectra in both directions is increased, e.g., peak values are around 0.75(y) and 0.5(z) at Sr ≈ 0.2 and 0.4(y) and 0.8(z) at Sr ≈ 0.3 . This is most likely attributed to the fact that the side loads were found to be dominant in one sporadically changing azimuthal orientation for short periods of time as was shown by Statnikov et al. [22, Fig. 13]. Hence, even though both directions should ultimately yield the same spectra due to the axisymmetry of the model, with a limited data collection time, less changes of the side load orientation can occur, and consequently, the side loads may appear more dominant in one direction. Interestingly, the average PSD over both directions appears much less affected by the data collection period which also supports the above argument. It also agrees with the literature data, with the main difference being a slight merging of the two distinct peaks that can be partially attributed to the lower  Fig. 7 Amplitudes of circumferential coherence modes C r,0 and C r,1 at x∕D = 0.72 . Comparison data taken from Weiss et al. [24] frequency resolution and total data collection period in the present evaluation. Note that the data by Statnikov et al. [22] were collected over an approximately 5 times longer period (approx. 410 CTUs) than in the current computation with extended data collection period. Even in their data, there is still some difference in peak location (e.g., Sr D ≈ 0.02 for the major peak) and amplitude (e.g., spPSD = 0.15 for the secondary peak). For the following investigation of the sensitivity of the obtained solution to parameter changes, the shorter data collection period of 38 CTUs is used, since for most of the considered quantities of interest, the impact is small. Additionally, this is a pragmatic choice, since a total of 8 different parameters are investigated, and thus, increasing the data collection period by a factor of two (or possibly more) is merely not economic. However, when evaluating differences in the more sensitive quantities like the rms pressure near the end of the nozzle, the behaviour of C r,0 at low frequencies, and the details of the buffeting loads, the possible impact of the data collection period will be shown.

Sensitivity of the mean flow quantities
In the following, the effect of the parameter changes on mean reattachment location, mean pressure distribution, as well as rms pressure distribution is analysed.

Sensitivity to time-step size
As described in Sect. 3, the choice of the time-step is guided by approximate best practice procedures based on the convective CFL number. However, the obtained time-step can differ depending on the selected target grid spacing location. Furthermore, the free stream velocity is not always a good indicator of the local velocity in the region of interest, adding uncertainty in this regard as well. Hence, in the following, the time-step is varied to investigate possible effects on the solution.
Two additional time-step sizes are tested; one reduced time-step size t 1 = 0.5 t ref and one increased time-step size t 2 = 2 t ref . The remaining settings are unchanged reference settings.
The resulting mean flow fields for both the reduced and the increased time-step differ in general only slightly from the reference solution. The mean reattachment location agrees to within 0.4% with the reference solution ( x∕D t 1 = 1.170 , x∕D t 2 = 1.177 ), and the obtained mean and rms wall pressure distributions shown in Fig. 9 also indicate only minor quantitative changes in the flow field. The mean pressure distributions are essentially indistinguishable for all three considered time-step sizes. Similarly, the pressure fluctuations illustrated by the rms value of the pressure differ only slightly, with the largest differences appearing in  The differences can partially be attributed to small changes in the fluid properties due to the different fluid modelling resulting, e.g., in slightly different transport coefficients. The viscosity, e.g., differs by about 0.3%. Other sources for the deviations are differences in the numerical treatments and algorithms. One such example is the computation of the internal energy at the face in Eq. (8) where numerical errors in the summation of the formation enthalpy can occur that are not present if just one species is used or the formation enthalpy is neglected. Similarly, the preconditioning algorithm also contains additional contributions for a more complex fluid model that could explain deviations. The deviations also match results from previous simple validation cases where small differences were also observed for different fluid modelling approaches that could be traced down to differences in the numerical operations [5]. Hence, even though both approaches should yield identical solutions at the considered conditions, numerical changes in the algorithms due to the expanded modelling capabilities are the most likely reason for the observed quantitative differences.
The pressure distributions for both the reference case and that with changed fluid model still agree satisfactorily with the experimental data. Furthermore, the average and maximum deviation in mean pressure coefficient between the reference settings and the more complex fluid model are only around 0.01 and 0.03 (relative pressure difference 0.4% and 1%), respectively. In terms of pressure fluctuations, the difference in rms pressure coefficient from the reference solution is on average 0.0017 (8%) with a maximum deviation of 0.005 (21.5%).

Sensitivity to circumferential resolution
The sensitivity of the solution to in-plane and circumferential resolution changes for a configuration without plume was previously investigated [18]. The investigation showed that the circumferential resolution drastically changed the flow field, but it was hypothesized that this behaviour might change when a configuration with an active plume is investigated. This hypothesis stated that the strong qualitative changes might be predominantly due to an insufficient resolution in the secondary recirculation region behind the second cylinder that is not present with an active plume, whereas the resolution in the main recirculation region is sufficient to capture the qualitative flow field.
Like in the previous investigation, the present fine grid features the same in-plane resolution as the grid used for the reference solution, but twice the amount of circumferential planes, i.e., the coarse grid has a resolution of 1.88 • and the fine grid of 0.94 • . The obtained mean and instantaneous the reattachment region where the increased time-step size predicts an increased level of pressure fluctuations by about 11% which is similar to the deviation observed for the longer data collection period.

Sensitivity to turbulence model
To test the sensitivity of the results with respect to the turbulence model, an additional solution is obtained with an IDDES approach based on the Spalart-Allmaras (SA) turbulence model [20]. Again, all other settings remain unchanged from the reference settings. The computational time with the SA-based approached is reduced by about 20% compared to the SST-based approach.
The mean and rms pressure distributions are also displayed in Fig. 9. Similarly to the changed time-step size, the mean flow field and pressure distribution remain essentially unchanged from the reference solution and the mean reattachment location is x∕D = 1.176 . The pressure fluctuations show small differences, especially in the reattachment region, with the one equation model showing a slightly lower overall level of pressure fluctuations that also deviates further from the experimental data. However, these differences are approximately of the same magnitude as the differences observed with a longer data collection period. The small differences between the different turbulence models might be intuitively surprising, especially considering that for RANS computations, the two models show significant quantitative differences in the computed flow field. For example, the reattachment location in axisymmetric RANS computations for the SA model is at x∕D = 1.028 , whereas the k-model predicts a reattachment at x∕D = 1.188 , leading to corresponding differences in the pressure distributions. However, one has to consider that in the IDDES approach, the majority of turbulent fluctuations are resolved and only the smallest scales and regions with attached wall bounded flow are dominated by the turbulence model. These parts are modelled well by both one-equation and two-equation models. Hence, if the resolution is sufficient, it is reasonable for both turbulence model approaches to yield similar results.

Sensitivity to fluid model
The resulting flow field for the computation with a more complex multi-species fluid modelling approach agrees qualitatively well with the reference solution. However, upon closer analysis, the mean reattachment location of the shear layer occurs slightly further downstream at x∕D ≈ 1.185 . This is also visible in the mean pressure distribution displayed in Fig. 9, where a downstream shift in the pressure distribution is observed. The distribution of the pressure Fig. 10 Mean (left) and rms (right) wall pressure coefficient distribution for the reference solution and the cases on the fine grid flow fields are qualitatively similar to that obtained with the reference settings in Fig. 2, but with the finer circumferential resolution, the reattachment occurs further downstream at x∕D∕approx 1.187 . Consequently, the mean pressure distribution in Fig. 10 is also slightly shifted downstream. Additionally, the level of pressure fluctuations is reduced for all axial locations.
These results differ significantly from those obtained previously for a configuration without plume [18]. With an active plume, the flow field for the coarse and fine resolution remain qualitatively similar with a clear mean reattachment on the surface of the second cylinder. Hence, only relatively small changes in the size of the recirculation region and in the pressure distribution are visible. In contrast, in the investigation without plume, the main recirculation region merged with the secondary recirculation region behind the second cylinder and thus changed the flow field qualitatively by completely removing a mean reattachment of the shear layer on the cylinder surface. Since the only difference between the configurations is the presence of a second recirculation region downstream of the second cylinder, it seems reasonable to deduce that the drastic flow field change observed for the configuration without plume is due to a strong sensitivity of this second recirculation region to the circumferential resolution. For the configuration with plume, these flow features do not occur and the overall flow field sensitivity to the circumferential resolution is lower.
However, the quantitative differences in the pressure distributions indicate that the solution is still relatively sensitive to the circumferential resolution. Hence, even though the agreement with the experimental data is slightly reduced, a circumferential resolution of at least 1 • seems necessary to capture the pressure distribution and reattachment length appropriately, since grid refinement objectively reduces the errors of the numerical scheme. A further improvement with continued refinement is possible, but seems uneconomical considering the non-dimensional circumferential resolution in the reattachment region would drop to values of z + < 100 . Additionally, circumferential resolutions reported in the literature are also not significantly finer (e.g. Refs. [22,24]).

Sensitivity to filter length definition
Finally, the effect of the filter length definition on the solution is evaluated. Since the general superiority of the ̃ definition over the definition is easily shown theoretically [13], this investigation is conducted using the fine grid. This allows not only an estimation of the sensitivity of the solution to the changed filter length definition, but also evaluates the general achievable accuracy with optimal settings (i.e., combination of best filter length definition and finest investigated grid).
The reattachment of the shear layer for this case occurs at x∕D = 1.169 . The wall pressure distributions can be found in Fig. 10. Interestingly, the reattachment location and pressure distributions are very similar to those obtained for the reference settings, i.e., on the coarser grid with the suboptimal length scale definition .
The main effect of the finer circumferential resolution is the reduced numerical dissipation due to the improved resolution. Since this reduced dissipation leads to a downstream shift of the pressure distribution, it seems logical that the upstream shift observed with the optimized filter length can be attributed to an increase in dissipation. The source of this increased dissipation can only be the eddy viscosity, since neither grid nor numerical scheme is changed. Comparing the solutions on the fine grid shows that the optimized

Spectral analysis
The wall pressure spectra for the different cases at the location x∕D = 0.6 are displayed in Fig. 11. The peak at the non-dimensional frequency Sr D ≈ 0.2 is clearly visible for all cases on the coarse grid, but the amplitude is slightly reduced. For the cases on the finer grid, the peak is even less distinct, and instead, the amplitude at higher frequencies is increased, with a clearer peak at Sr D ≈ 0.35 . The reduced peak for the fine grid is due to a redistribution of the energy to higher frequencies 0.3 < Sr D < 1 that make up about 31% of the total energy for this case, whereas for the reference case, it is only 22% at this axial position. For the case with the changed filter length, there is also a higher content in this higher frequency range (29%).
However, the detailed behaviour significantly depends on the axial location that is being evaluated as is evident from the additional spectrum at x∕D = 0.4 for the changed filter length case shown. At this location, the peak at Sr D ∕approx 0.2 is a clearly dominant peak. This behaviour is further visualized in Fig. 12 where the distribution of the scaled premultiplied PSD at Sr D = 0.212 is displayed for different axial locations. It is clearly visible that for the case with optimized filter length, the peak at this frequency appears with a similar intensity as in the reference case at a further upstream ( x∕D ≈ 0.45 ) position. Hence, the peak does not disappear completely, but is merely disrupted at the evaluated position in Fig. 11 by other frequency content. On the other hand, for the case using the fine grid without optimized filter length, the peak is only very faintly visible around x∕D ≈ 0.6 , and instead, contributions at higher frequencies are increased over nearly all axial positions. This behaviour might be explained similarly to the shift in reattachment location by the excessive reduction of eddy viscosity that allows high-frequency disturbances to persist, whereas these are more damped on the coarser grid or with an improved definition of the filter length.
Towards the end of the nozzle at x∕D = 1.15 , the spectra for all cases look very similar, with no distinct peaks but a broadband increase of spectral content between 0.3 < Sr D < 1 (see Fig. 13).
The evaluation of the wall pressure in terms of circumferential modes of the coherence function in Fig. 14 shows that for all cases on the coarse grid (left), a peak of C r,1 at Sr D ≈ 0.2 and a smaller one around Sr D ≈ 0.35 are visible. However, the amplitude of the main peak is smaller than in the reference solution for all cases. On the fine grid (Fig. 14,  right), the peak is again less defined, even completely disappearing on the fine grid with the reference filter length definition. A similar behaviour as for the pressure spectra can be observed for the case with improved filter length definition, filter length indeed increases the average eddy viscosity in the recirculation region by about 10%, which explains the upstream shift of the pressure distribution. Consequently, it seems that the better solution on the coarse grid is most likely due to a correct total dissipation for the wrong reasons: the excess numerical dissipation makes up the lack in turbulent dissipation, leading to an approximation of the total dissipation on the coarse grid which otherwise can only be obtained with a more sophisticated combination of grid resolution and filter length definition. This finding is a good example showing how important systematic investigations of grid and modelling parameter sensitivities are.
with only a minor peak appearing at the evaluated position x∕D = 0.72 , but a more pronounced one further upstream (cf. Fig. 12).
The spectra of the buffeting forces displayed in Fig. 15 show that the peak around Sr D ≈ 0.2 is visible for all cases with the exception of the case on the fine grid with improved filter length. For that case, the force spectra show a small peak between Sr D = 0.1 … 0.2 and a strong peak around Sr D ≈ 0.4 instead. It can also be observed that the detailed spectra differ quite significantly for all cases with some showing two very distinct peaks at Sr D ≈ 0.18 and Sr D ≈ 0.35 (time-step change, fine), whereas for others, only  14 First two modes of the circumferential coherence for the cases on the coarse (left) and fine (right) grid one peak around Sr D ≈ 0.23 (turbulence model, fluid model) appears. Additionally, the peak height for all cases is less than that of the reference solution, indicating less energy of the buffeting force and instead more energy at both small and high frequencies. The large deviations in the buffet loads even for those cases that showed very good agreement in the observed wall pressure distributions and spectra indicate that the buffeting loads are one of the most sensitive quantities in the analysis and likely also the quantity that most benefits from a longer data collection period. Hence, the obtained sensitivities for this quantity might be not completely representative and have to be interpreted with caution.
Increasing the data collection period to 76 CTUs for the case with optimized filter length (Fine + filter length (long)) confirms this. The spectra change significantly, shifting the main peak towards Sr D ≈ 0.21 and the secondary peak towards Sr D ≈ 0.37 , also improving agreement with the literature data; that this is indeed a converging behaviour is confirmed by an additional case with a further increased data collection period of 114 CTUs, i.e., three times the initial data collection period. It is visible that the spectra for the cases with 76 and 114 CTUs are much more alike, indicating that the initial data collection period of 38 CTUs is not necessarily sufficient for a qualitative evaluation of the buffeting force spectra.
However, all other investigated quantities are much less sensitive and a further extension of the data collection period does not change the results discussed in Table 3 (reattachment

Conclusions
The sensitivity of the transonic flow around a generic space launch vehicle geometry to different numerical parameter variations is investigated. None of the investigated parameter variations changed the qualitative flow field, and for several of the investigated parameter, even quantitative changes were minuscule.
The exact choice of the time-step size to within a factor of 2 larger or smaller and the choice of the underlying turbulence model essentially does not affect first-and secondorder statistics. The spectral content of the wall pressure also agrees very well with the reference solution for these cases by displaying the same dominant peak with very close agreement in terms of frequency and amplitude.
A changed fluid modelling as well as a finer circumferential resolution shift the reattachment location about 0.01D further downstream, whereas an optimized filter length shifts it approximately that same distance further upstream. Similarly, the mean wall pressure distributions are shifted downstream or upstream, respectively. However, with the changed fluid model, the spectral content of the wall pressure is only slightly affected, whereas the grid refinement and filter length definition change the spectra significantly.  A spectral analysis of the buffeting loads shows large differences between all cases, with the case on the fine grid with an optimized filter length showing the largest deviations. However, it is also shown that the buffet loads appear to be the quantity most sensitive to the data collection period, making quantitative changes in peak frequencies and amplitudes difficult to interpret.
As a consequence, it is recommended to always ensure a sufficient circumferential grid resolution and a state-ofthe-art filter length definition is used to obtain accurate results. On the other hand, the exact choice of time-step size and turbulence model can be chosen more freely and other considerations such as computation effort should be taken into account when making these choices. One also has to take into account that the change of the fluid model can affect the mean results. For the analysis of mean flow features and pressure spectra, a relatively short data collection period (here: 38 CTU) is sufficient, but for the analysis of force spectra, a significantly longer period (2-3 times, here: 76-114 CTUs) is necessary for accurate results and has to be planned for accordingly.
Equipped with and employing these findings and recommendations, it is possible to investigate space launch vehicle aft-body flows with high confidence in the numerical results. This allows to examine the impact of different parameter variations on the solution and observed loads. Interesting effects are, for example, expected from changes in the geometry such as nozzle length, from variations in exhaust plume conditions in terms of fluid properties and nozzle pressure ratio and from increased wall temperatures.