Turbulent impinging jets on rough surfaces

This work presents direct numerical simulations (DNS) of a circular turbulent jet impinging on rough plates. The roughness is once resolved through an immersed boundary method (IBM) and once modeled through a parametric forcing approach (PFA) which accounts for surface roughness effects by applying a forcing term into the Navier–Stokes equations within a thin layer in the near‐wall region. The DNS with the IBM setup is validated using optical flow field measurements over a smooth and a rough plate with similar statistical surface properties. In the study, IBM‐resolved cases are used to show that the PFA is capable of reproducing mean flow features well at large wall‐normal distances, while less accurate predictions are observed in the near‐wall region. The demarcation between these two regions is approximately identified with the mean wall height km of the surface roughness distribution. Based on the observed differences in the results between IBM‐ and PFA‐resolved cases, plausible future improvements of the PFA are suggested.


INTRODUCTION
Technical fluid flows occur over surfaces that are rough as a result, for instance, of the manufacturing process or wear. From the fluid-dynamics perspective, a wall is considered either rough or smooth according to the flow characteristics. For wall-bounded turbulent flows, the interaction of the smallest scales of motion with roughness elements usually results in enhanced wall-friction and heat exchange coefficients with respect to a comparable flow over a smooth wall. While engineering tools for the prediction of drag and heat transfer on rough surfaces exist, they are typically calibrated in canonical wall-bounded turbulent flows, like pipe or channel flows, where a large amount of reference data dating back to the pioneering work of Nikuradse, Moody and Colebrook is available [16]. Rough wall reference data for more complex flow conditions is much more sparse. A thorough review of turbulent flows over rough surfaces is presented in Chung et al. [2].
The present study considers the impingement of a turbulent jet onto a surface characterized by a random roughness type distribution (in contrast to regularly arranged surface structures). While the drag increasing effect of roughness in internal flows, like pipes or channels, is often assessed indirectly through pressure drop measurements, this is not possible for the impinging jet configuration. The effect of roughness is thus experimentally assessed through optical flow field measurements of such a stagnation flow. For the corresponding numerical analysis, investigating the effects of roughness via fully resolved numerical simulations of turbulent flows over rough walls remains a challenging task, despite the constantly increasing computational power. In this context, immersed boundary methods (IBM) have become standard practice for simulating complex boundaries in direct numerical simulation (DNS) (see for instance [8], [17]). Despite their versatility, IBM require heavy mesh refinement when complex geometries are involved. In consequence, resolving roughness geometries within Cartesian grids may exceed the grid resolution requirements demanded to capture the smallest eddies in the flow. One possible alternative is the parametric forcing approach (PFA) introduced by Busse and Sandham [1]. Similarly to IBM, no body-conform mesh is required to describe the roughness. The PFA adds a forcing term into the Navier-Stokes equations, which is active only within a thin layer where surface roughness is present. The added unit volume force serves to account for the form drag induced by roughness elements. The modified PFA presented by Forooghi et al. [6] includes an additional viscous drag contribution to the original formulation. The same work shows the successful application of the method for predicting the flow characteristics of fully developed turbulent channel flows with rough walls. Recently, the PFA roughness model was also successfully applied to turbulent flows over localized (heterogeneous) roughness [15]. In the present work, IBM and PFA approaches are both compared among each other and to experimental data, providing novel reference data for spatially developing flows over rough surfaces. Figure 1 depicts a sketch of the flow configuration. The jet originates from a fully developed turbulent pipe flow and impinges on a plate placed 2D away from the jet exit section (D represents the pipe diameter). The Reynolds number is Re = 10 000 (Re is based on the pipe diameter and the bulk mean velocity in the pipe). The flow is incompressible and free from volume forces, except in the near-wall region where a distributed forcing is applied to introduce the effects of surface roughness into the flow. Fully developed inflow boundary conditions are enforced at the inlet section of the computational domain using a periodic precursor simulation of a fully developed turbulent pipe flow with a length of 12D. We note that inflow conditions could alternatively be generated with a synthetic eddy method. However, these methods require a spatial development length that differs for different statistical properties of the turbulent flow. Since impinging jets are known to be sensitive toward inflow conditions [5] the precursor simulation was chosen instead to ensure clearly defined inflow conditions. In the numerical simulations, a confinement plate is placed at the same height as the nozzle exit section (see Figure 1 on the right). The adoption of a confinement plate greatly simplifies the modeling of the flow configuration since it prevents entrainment of fluid from outside of the computational domain. Both the impingement and confinement plates have a diameter of 20D. The large size of the computational domain is adopted in order to minimize the influence of the outflow boundary condition on the flow region of interest. The outflow boundary condition is applied on the lateral surface of the computational domain following the stabilizing outflow boundary condition presented by Dong [3]. No-slip and nonpenetration boundary conditions are applied for the velocity field on both the confinement and the impingement plates. The latter corresponds to the bottom of the computational domain and surface roughness is deposited homogeneously over it. The flow solver employed is Nek5000 [4]. It is an open-source code based on the spectral element method proposed by Maday and Patera [12]. Nek5000 is well known for its efficient parallel scaling characteristics and for retaining high-order spectral accuracy in the solution.

Flow configuration
Time advancement is performed using a third-order accurate implicit-explicit scheme in which linear terms are treated implicitly, while nonlinear terms are treated explicitly. Flow field statistics are accumulated in time on-the-fly during each simulation in order to avoid the storage requirements of saving a large number of instantaneous samples for later post-processing. Averaging in the circumferential direction is successively applied and the resulting mean flow field is two dimensional, as any mean flow property depends only on the radial distance from the jet axis r and the wall-normal distance x 3 . In what follows, an over-bar (⋅) is used to indicate the process of averaging in both time and circumferential

Computational mesh
The computational mesh adopted for the jet flow domain consists of 2 034 995 spectral elements, while the precursor pipe-flow mesh counts a total of 608 175 elements. Details of the mesh are reported in Figure 2 and in Table 1 for the jet flow domain. The pipe mesh parameters in the cross-sectional area of the pipe are the same as those reported in Figure 2B and, in addition, a uniform elements distribution, consisting of 256 elements, is used in the streamwise direction for the pipe.
For a fixed number of elements, the spatial resolution is set by the polynomial order of the local element basis onto which the solution is discretized. In particular, a seventh-order solution is adopted for all the IBM-resolved cases, while a fifth-order solution is chosen for the PFA-resolved cases. In the former case, the higher resolution is exploited in order to better resolve the geometrical features of the surface roughness in the simulation. For all the cases reported in this study, a P n − P n formulation has been used. To check whether the adopted resolution is sufficient to resolve all relevant scales of motion, the local grid size is compared to the Kolmogorov length scale = ( 3 ∕ ) 1∕4 , where indicates the kinematic TA B L E 1 Mesh parameters of the jet flow domain. All the parameters reported in the table are depicted in Figure 2 for clarity. The first column in the table indicates the inner radius of the mesh, while the other columns report the number of elements in the respective direction  Figure 3 reports the distribution of the ratio (Δx 1 Δx 2 Δx 3 ) 1∕3 ∕ in the x 1 − x 3 plane of the jet flow domain for a representative case of both the IBM-and PFA-resolved cases (here Δx i , with i = 1, 2, 3, represents the local grid size in the coordinate direction i). From the inspection of the figure, the grid resolution can be deemed appropriate for both cases, as the local grid size to Kolmogorov length scale ratio is always smaller than ≈ 3 in all flow regions of interest. Differences observed in Figure 3A and B are mainly due to the different resolution adopted for the IBM-and PFA-resolved cases, respectively. A grid independence check suggests that even the coarser computational grid employed for the PFA cases is sufficient for a consistent reproduction of basic statistics of the velocity field in the region of interest of the flow. For instance, Figure 4 represents the comparison between third-and fifth-order solutions for one representative PFA-resolved case. The figure compares various mean radial velocity profiles at different radial locations. No major discrepancies between the two different order solutions can be appreciated at any location, suggesting that the adopted resolution suffices in resolving accurately the mean flow features.
A similar grid independence study was not performed for the IBM-resolved cases due to prohibitive computational costs. Because of the resolution requirements, it is estimated that the computational cost associated with the IBM-resolved cases is approximately three times greater than that associated with the PFA-resolved cases. Indeed, within the framework described so far, a complete DNS using the IBM requires ≈ 920 ⋅ 10 3 core-hours, while a DNS with the PFA uses only ≈ 290 ⋅ 10 3 core-hours. Therefore, the good agreement of IBM-resolved results with experimental data (shown in Sections 4.1 and 4.2) is assumed to be sufficient proof that the adopted grid resolution is appropriate.

Parametric forcing approach
The approach followed in the present work is equivalent to the "modified PFA" introduced by Forooghi et al. [6]. A forcing term f i is introduced into the momentum equation of the incompressible Navier-Stokes equations: The above equations are valid for (i, j = 1, 2, 3) and summation over repeated indices is implied. Equations (1) and (2) are made dimensionless using the mean bulk velocity U b in the pipe and the pipe diameter D. Therefore, the Reynolds number appearing in Equation (2) is defined as Re = U b D∕ , where indicates the kinematic viscosity of the fluid. Within the present formulation, forcing is applied only in the wall-parallel directions x 1 − x 2 , therefore f 3 = 0. Furthermore, it is assumed that the volume force f i is made up of a linear (with the local fluid velocity) and a quadratic term. The former represents the viscous drag contribution of roughness elements, while the latter accounts for their form drag. By introducing two shape functions A and B, the volume force can thus be expressed as: The functional dependency of the functions A and B on the wall-normal direction x 3 only is specific for the present case and will be clarified in the following sections. The analogy with the Darcy-Brinckmann-Forchheimer equation in the context of flow through porous media is evident in Equation (3) (see, for instance, [18]). Furthermore, the similarity with modeling of flow through porous media is additionally exploited to determine the function A appearing in Equation (3). As derived in Forooghi et al. [6], A is expressed as: where * is the viscosity of the fluid, while K * is the permeability of the porous medium (asterisks in Equation (4) denote dimensional quantities). Different approaches can be followed in order to relate the permeability K to the geometrical features of the porous medium or, in this case by analogy, of the rough wall geometry. One possibility, as suggested by Forooghi et al. [6], is to use the Kozney-Carman theory and, hence, the similarity between flow through porous media and flow through bundles of capillary conduits (the Kozney-Carman theory is explained, for example, in Kaviany [11]). By following this approach, the proper nondimensional form of the function A is found to be: In Equation (5), K k is the Kozney constant arising from the adoption of the Kozney-Carman theory, s(x 3 ) is the total wetted surface area per unit total volume, and (x 3 ) is the porosity. In particular, the porosity is defined as the ratio of the fluid volume to the total volume. Both the wetted area and the porosity are here considered functions of the wall-normal direction x 3 because the surface roughness is modelled homogeneously in the x 1 and x 2 directions.
By recalling the physical meaning of the quadratic term in Equation (3), it is straightforward to identify the shape function B using the analogy with the form drag induced by a bluff body. This results in: where C d is interpreted as the drag coefficient, and s f (x 3 ) is the total frontal projected area (averaged in the x 1 and x 2 directions) of the roughness per unit total volume.

Method
Parametric forcing approach Immersed boundary method = 2 ⋅ 10 5 = 200 From Equations (5) and (6), it is observed that the shape of the two functions A and B is completely determined by the geometric characteristics of the roughness. The constants K k and C d appearing in the definition of the volume force f i can be used to tune the model. In this respect, it is remarked that the model constants could be also related to the features of the rough surface in order to have a fully predetermined model which does not necessitate any ad hoc tuning. However, a clear connection between these constants and geometrical features of the surface roughness is not yet available and, for this reason, the constants are considered as free parameters in this work.

Immersed boundary method
Reference DNS are obtained by using the IBM of Goldstein et al. [8]. Similarly to the PFA, this particular IBM introduces a forcing term to the momentum equation. The forcing is applied only for grid points lying within the solid phase of the computational domain, thus a complete description of the solid body geometry can be achieved with the IBM, provided that the grid resolution is sufficiently fine. The forcing term appearing in Equation (2) is in this case: where and are two constants that can be adjusted, and U body,i is the ith component of the prescribed velocity of the solid body. In the present study, only stationary walls are considered, therefore U body,i = 0 for i = 1, 2, 3. The volume force introduced in Equation (7) acts as the input of a proportional integral controller which aims at keeping the local fluid velocity to the prescribed value U body,i (thus, in the considered case, U body,i = 0).

Models constants
As it is clear from the preceding discussion, two constant parameters must be set for both the PFA and the adopted IBM.
Regarding the IBM, this topic is well covered in the original work of [8] where it is shown that the order of magnitude of and needs to be properly adjusted for each considered flow configuration, but the actual value of the two constants is irrelevant on the final result. On the assumption that the response of the system to an input force as that in Equation (7) can be approximated by the response of a linear system, the parameters and assume a clear significance. In particular, the integral constant determines the natural frequency of the system, while the proportional constant is related to the damping ratio. Hence, large values of allow to track fast fluctuations of the local fluid velocity, while large values of help to damp out oscillations. Clearly, the major limitation encountered when setting the two parameters is the numerical stability of the time marching technique adopted. Similar reasoning can be performed also for the parameters K k and C d of the PFA, but no clear meaning can be attached to the two constants in this case. Since no prior knowledge about the PFA applied to impinging jet configurations is available, model constants used in this study are picked without any modification from turbulent channel flow simulations in which the PFA was successfully applied to account for surface roughness effects [6]. The set of constants adopted for all the presented cases is reported in Table 2.

Surface roughness
Rough surfaces consist of wall-height distributions in the x 1 − x 2 plane having prescribed statistical properties. The algorithm employed to design roughness geometries is completely based on the procedure presented in Pèrez-Ràfols and Almqvist [14]. This is an iterative approach which allows to generate roughness topologies having prescribed power spectrum (PS) and probability density function (PDF). The algorithm starts from a surface roughness distribution characterized by a prescribed PDF and, in general, unknown PS. The initial wall-height distribution is transformed in the . k 99 ∕D = 0.12; Fourier space where a correction is applied in order to match the desired PS. Successive transformation back into the physical domain leads to a modified wall-height distribution which has, possibly, a PDF that differs from the one specified at the step zero of the algorithm. A correction in the physical space is hence performed to match again the prescribed PDF. The procedure is repeated iteratively until convergence is reached for both the PDF and the PS. Further details about the algorithm can be found in Pèrez-Ràfols and Almqvist [14]. By scaling properly the surface height distribution obtained with the algorithm just presented, it is possible to characterize different rough surfaces. In this respect, the main quantity adopted in the present work to distinguish between different rough surfaces is the k 99 value (i.e., the 99% confidence interval of the PDF distribution of roughness height). Three different rough surfaces are investigated numerically in this study; these are characterized by k 99 = 0.05D, k 99 = 0.12D and k 99 = 0.15D. A Gaussian PDF is adopted for all the cases, while the PS is prescribed to match the PS of a realistic rough surface as in Forooghi et al. [7]. The effective slope parameter (ES) (as defined in Napoli et al. [13]) is kept constant ES = 0.41 for the three cases. Samples of the obtained roughness geometries are displayed in Figure 5. The x 3 ∕D = 0 plane is located at the lowest valley of the roughness height distribution. This choice of origin is consistent with surface roughness caused by deposition in real applications. The three rough surfaces are modeled within the PFA via the two roughness shape functions A and B defined in Equations (5) and (6), respectively. The normalized shape functions for the three investigated surfaces are reported in Figure 6. From the figure, it is evident that the function A, which tunes the effect of the viscous drag term of the forcing (i.e., the linear term), is relevant only over a thin layer very close to x 3 ∕D = 0. On the other hand, the function B has a wider support which spans the entire thickness of the roughness height distribution. It is remarked that the overall contribution of the form drag term of the PFA forcing (i.e., the quadratic term in Equation 3) becomes weaker close to the wall where the local fluid velocity tends to zero.

REFERENCE EXPERIMENT
The experiment is designed to produce highly similar flow conditions to the DNS setup. The experimental realization of the axial symmetric impingement jet configuration is shown in Figure 7. A circular, seamless stainless steel tube with an The impingement plate with a diameter of 12D (= 300 mm) is made of aluminum and is kept at the same temperature as the air jet. In addition to an impingement plate with a smooth surface (milled, roughness R a < 3 μm), a wall with a precisely defined rough surface was manufactured using a 3D printing process (Selective Laser Melting, FIT Prototyping GmbH). The desired surface topology was generated according to the process described in Section 2.5, transferred to the 3D CAD model of the wall and printed in layers of 50 μm thickness. The process was previously tested on smaller model geometries and the surface structure was checked using white light interferometry. Deviations from the specified surface structures were smaller than 75 μm on average. The statistical properties of the rough surface topology used in the experiments (k 99 ∕D = 0.04) are comparable to that in Figure 5A.
Velocity fields are measured using particle imaging velocimetry (PIV). The light source is a frequency doubled Nd:YAG double cavity laser (Spectra Physics PIV 100) firing two pulses with a variable delay Δt with a repetition rate of 10 Hz. The laser light sheet is spanned by a cylindrical lens system and has a height of ≈ 40 mm and a thickness of ≈ 65 μm (1∕e 2 radius) at the waist. The air stream is seeded with di-ethyl-hexyl-sebacat oil droplets using a seeding generator (Palas AFG 10) with an average droplet size of 0.5 μm. The scattered light is detected by a sCMOS camera (LaVision sCMOS Imager) with 2560 × 2160 pixels and a 16 bits analog-to-digital converter. A Nikon 105 mm f ∕1.8 objective lens is used for imaging. Close-up rings of 12, 36, and 48 mm are employed for different spatial resolutions and fields of view. Camera and laser sheet are placed in a Scheimpflug arrangement (viewing angle approximately 2 • -3 • ), in order to minimize spatial clipping of the scattered light by the wall and to facilitate measurements very close to the wall. Commercial software (LaVision Davis 10) is used for image acquisition, parameter setting, perspective correction, and image correlation. The time delay between the double images has been set according to the correlation window size and the range of velocities present in the experiment. Table 3    Typically each measurement series consists of about 1000 vector fields. The PIV data presented in this paper consists of multiple independent measurements. They differ in the field-of-view, Δt, the mass flow controllers used, and multiple re-alignments of the optical setup to account for systematic errors. The bulk velocity U b is determined for each measurement individually by integrating the velocity profile at the outlet. It varies by less than 3.5 % between measurements. All velocities reported hereafter are normalized by the bulk velocity of the respective measurement. Finally, based on a statistical analysis close to the nozzle exit, the statistical error (confidence interval) of the major velocity component is < 1.5 % for the mean (depending on the turbulence level) and between −7% and +15% for the root-mean-square values at a confidence level of 99 %.

Validation of the DNS setup: smooth wall case
The case of a turbulent jet impinging on a smooth flat plate is considered at first in order to check the adequacy of the numerical setup and of the experimental framework adopted in the study. The flat wall is resolved in the DNS using the IBM. This unconventional approach of modelling the flat wall through IBM allows to confirm the IBM implementation directly within the present DNS set-up. In addition, in this case the IBM forcing can be directly linked to the wall shear stress, that is, the momentum transport from the fluid to the wall, which is exploited for a direct comparison with the PFA model in Section 4.3. Flow field statistics are computed by averaging in time (on the fly during the simulations) and in the circumferential direction. Hence, the resulting averaged flow field is two dimensional, with any statistical quantity being a function of the radial distance r from the jet axis and the wall-normal distance x 3 . In what follows x 3 and z will be used interchangeably to indicate the wall-normal distance.
In the comparison between the DNS results and the experimental measurements, discrepancies between mean flow variables could be expected due to the slightly different configurations adopted in the two cases. Indeed, a confinement plate is adopted in the DNS, but not in the experiments. Nevertheless, for the examined nozzle-to-plate distance of 2D, the influence of the confinement plate on the near-wall mean flow is assumed to be negligible with respect to the uncertainties associated with the loss of accuracy of the measurements in this area. This assumption is confirmed by the very good agreement observed everywhere away from the near-wall region where experimental measurements provide much more accurate data.
Experimental and DNS inflow conditions are reported in Figure 8. The figure shows the comparison of basic flow statistics sampled in the symmetry plane at 1 mm distance from the jet exit section. The mean axial velocity profile u z ∕U b is displayed in Figure 8A. Very good agreement is observed between the measured and computed mean axial velocity profiles at all radial locations. Reynolds stress tensor components u z u z ∕U 2 b and u r u z ∕U 2 b are depicted in Figure 8B and C respectively. They are also in very reasonable agreement between measurements and simulations with slight indications of nonperfect symmetry conditions in the experiment.
Mean radial velocity profiles close to the smooth surface and at different radial locations are presented in Figure 9A.   Figure 9B shows the boundary layer growth along the impingement plate. As it is customarily done in the analysis of wall-jets, the inner layer height is identified as the wall-normal location at which the radial velocity reaches its maximum value (see for instance Wu et al. [19]). Similar to the inner layer, an outer layer thickness is defined as the greatest wall-normal distance at which the radial velocity reaches one half of its maximum value at the same radial location. Very good agreement between the numerical and experimental results is observed for, both, the inner and outer layer development. The wall-stress distribution along the plate is reported in Figure 9C. For the DNS, the wall-stress at a fixed radial location is computed from the mean radial volume force introduced by the IBM and by applying a control volume analysis to a small volume that encloses the considered radial location. On the other hand, for the PIV measurements, the wall-stress is computed from the slope of the mean radial velocity profile at the wall. The latter is estimated from a polynomial fit of the radial velocity profile close to the wall and that is forced to be zero at the position of the wall. The DNS-computed and PIV-measured wall-stress distributions agree reasonably well at all radial locations. In Figure 9C, the PIV-measured wall-stress values are reported along with 3 error bars (here indicates the standard deviation of the slope as obtained from the polynomial fit).

Validation of the DNS setup: rough wall case
The present numerical framework has also been validated by further considering an impinging jet over a homogeneous rough surface. To this aim, the impingement plate with prescribed surface roughness was 3D-printed, as described in Section 3 of this work. The roughness distribution chosen for this case has the same statistical properties as for the surfaces presented in Section 2.5. However, the wall height distribution was scaled in order to have k 99 ∕D = 0.04. The same surface roughness is resolved in the simulation using the IBM. The comparison of the DNS results and PIV measurements for this case is reported in Figure 10. Mean radial velocity profiles at different radial locations are depicted in Figure 10A Here k m indicates the mean wall height of the roughness distribution. Experimental and computed profiles match very well at all the considered radial locations. Nonetheless, it is noted that assessing experimentally the near-wall region of the rough surface is challenging and, for this reason, it is difficult to obtain accurate measurements in this area. Similar to the smooth wall case, Figure 10B reports the development in the radial direction of the inner and outer layers for the rough wall case. Consistent with the radial velocity profiles in Figure 10A, discrepancies are visible for the inner layer height, especially at small radial distances from the jet axis and directly reflecting the larger uncertainties close to the rough surface. Some systematic deviations are also visible for the outer layer at r∕D > 2.5. The Reynolds shear-stress u r u z ∕U 2 b distribution at x 3 ∕D = 0.06 and x 3 ∕D = 0.1 is shown in Figure 10C. The distribution for the smaller wall-normal distance shows very nice agreement between the DNS results and the PIV measurements. Visible discrepancies are observed only near r∕D = 1. At greater wall-normal distance, that is, at x 3 ∕D = 0.1, deviations are more pronounced for 1 < r∕D < 2, and, especially, a larger discrepancy is observed in the peak value of the distribution. These differences might be ascribed to experimental uncertainties. For example, the exact determination of the location of the rough surface in the experimental images is associated with a larger uncertainty (Δx 3,wall ≈ ±0.008D). In addition, the PIV measurements currently lack radial averaging of the velocity profiles, making the results sensitive to the influence of local roughness elements that can affect the flow field at larger radial distances (compare the discussion in Section 4.3 related to Figure 12). Overall the agreement between experiment and DNS is reasonable and differences are not considered crucial to invalidate the present numerical setup.

PFA assessment
The PFA is compared here to fully resolved surface roughness cases that were simulated with the IBM. In particular, we present results obtained by applying the two methods to the three roughness geometries described in Section 2.5. The mean flow field is investigated in order to reveal possible areas where the PFA might have shortcomings. A first comparison is done for the mean radial velocity profiles. After the jet impingement, a radial wall-jet develops along the rough plate where a boundary layer develops and thickens rapidly due to the entrainment of the outer quiescent fluid.
Mean radial velocity profiles at different radial locations are depicted in Figure 11 for the three rough wall cases. The velocity profiles are plotted using a logarithmic scale for the wall-normal coordinate x 3 to better visualize the near-wall region. In the figure, a black dashed line represents the mean wall height of the roughness distribution, while blue and red lines represent IBM and PFA-resolved results, respectively. The velocity profiles shown in the figure correspond to the radial locations r∕D = 0.6, 0.8, 1.0, 1.2, 1.4, 1.5, 2.0, but they are plotted shifted away from each other by 1.0 units for better visualization. It is observed that PFA-resolved profiles display a non-zero velocity at all the locations below the mean roughness height of the considered rough surface. This is in accordance with the PFA, and contrary to the IBM adopted, as the applied forcing does not aim at enforcing a zero velocity within the forcing layer. At the same time, the large discrepancies observed especially for profiles close to the jet axis can be ascribed to the lack of wall-normal forcing for the PFA. As discussed in Section 2.2, in the present study f 3 = 0, consistent with the original version of the PFA presented in Busse and Sandham [1] and Forooghi et al. [6].
The near impingement region is characterized by a strong wall-normal velocity component. In this area, the wall-normal forcing applied by the IBM acts effectively to damp the wall-normal velocity in order to guarantee zero flux through the rough wall. Hence, on average, the deflection of the wall-normal jet into a radial wall-jet takes place at wall-normal distances greater than the roughness mean wall height.
The behavior just described can be clarified by considering how the IBM-resolved flow develops radially at different wall-normal heights. Figure 12 depicts the time averaged distribution of the radial velocity component on planes at different wall-normal heights for the IBM-resolved case with the largest roughness size (i.e., the roughness distribution having k 99 ∕D = 0.15). In particular, slices at x 3 ∕D = 0.2, 0.18, 0.16, 0.11, 0.10, 0.05 are shown in Figure 12A,B,C,D,E,F, respectively. Contours of roughness elements at the same wall-normal height are visualized with the aid of white solid lines. The mean height of the roughness distribution for the case reported in the figure is k m ∕D ≈ 0.13; therefore, the three panels at the top of the figure (Figure 12A-C) represent planes above the mean wall height. At these locations, the  Figure 12A where only few roughness peaks are present. Moving to lower wall-normal locations, the number of roughness elements increases, as a larger number of elements is cut by a plane at constant wall-normal height ( Figure 12B,C). Correspondingly, the radial wall-jet is affected heavily by the influence of the multitude of roughness elements. This is especially visible in Figure 12C, where the presence of roughness is responsible for the appearance of radial streaks characterized by high radial velocity component. At wall-normal heights below the mean roughness height, a rather different situation is observed. Examples can be seen in Figure 12D-F. At these locations, the distribution of the time averaged radial velocity component is characterized by nonzero values only in the "valleys" of the surface roughness. Shifting the considered plane of observation toward smaller wall-normal heights, the likelihood of intersecting regions of fluid (i.e., the valleys of the surface roughness) becomes smaller and, therefore, the radial velocity distribution is zero almost everywhere ( Figure 12E,F). Then, it is not surprising to observe that the circumferential average of the radial velocity component is very close to zero at all wall-normal locations below the mean wall height, as shown in Figure 11. The same behavior cannot be well reproduced with the application of the present PFA, where the absence of wall-normal forcing results in the nonpenetration boundary condition being satisfied only at the bottom of the computational domain, that is, at x 3 ∕D = 0. Therefore, radially directed flow occurs beneath the mean wall height for all the PFA cases, determining the visible discrepancies observed in the comparison with the IBM-resolved cases.
Nonetheless, for x 3 ∕D > k 99 ∕D, a good qualitative agreement between the shapes of the radial velocity profiles is observed in Figure 11 at all radial locations for all the three cases. In particular, values and locations of mean radial velocity peaks are well estimated for a wide range of radial locations by the PFA-resolved solution. Inner and outer layer distributions along the radial direction are depicted in Figure 13. For the k 99 ∕D = 0.05 case, the agreement on the inner layer thickness between the PFA and IBM is good for r∕D < 3, while the discrepancy becomes more evident for larger radial locations. From the figure it is also noted that the thickness of the outer layer is also well estimated by the PFA for the case with k 99 ∕D = 0.05 for all radial locations r∕D < 2.5. For the other two cases having a larger k 99 , the estimation of the inner and outer layer thicknesses with the PFA is also reasonably good, but departures from the IBM-resolved solution start already at smaller radial distances (r∕D > 2 for both the inner and outer layer approximately). Inner and outer layer distributions suggest that the PFA-estimated radial wall-jets decay more quickly compared to the IBM-resolved cases. In fact, it is evident from the figure that, for all the cases, the thickness of the two layers is overpredicted by the PFA for increasing radial distances. Hence, comparing the PFA and IBM predicted wall-jets at a fixed radial location, the former appears thicker and slower (in the sense that the maximum radial velocity estimated by the PFA is lower than the IBM-estimated one). This reasoning suggests that, overall, the PFA introduces excessive total wall-stress, which is responsible for the faster decay of the wall-jet. Here, total wall-stress refers to the shear stress at the wall, which includes skin-friction and pressure drag contributions.
The analysis of the inner and outer layers gives information about how, in general, the presence of surface roughness affects the growth of the boundary layer along the impingement plate. In this respect, the behavior can be inferred by observing Figure 14, which depicts the inner and outer layer distributions along the plate for the examined rough wall cases and the smooth wall case. Only IBM-resolved cases are reported in the figure. In order to allow the comparison with the layers height of the smooth case, the height of the layers of the cases with roughness has been shifted downward by the respective mean wall height k m ∕D. It is stressed that this is a rather arbitrary choice, as a different zero-wall location could be chosen for the rough wall cases. For radial wall-jets over rough surfaces, Wu et al. [19] adopt the zero-plane displacement thickness proposed by Jackson [9]. However, in the current analysis, the choice of identifying the zero-wall height with the mean roughness height k m ∕D is supported by the fact that, for all the simulated rough cases, the mean radial velocity profile is zero almost everywhere for 0 ≤ x 3 ∕D ≤ k m ∕D, as already noticed in the analysis of Figure 11.
From Figure 14 it is evident that surface roughness has a significant effect on the development of both the inner and outer layers. In general, the two layers appear thicker for rough surfaces of larger size. This is a result of the greater skin-friction along the plate experienced on average by the radial flow over a rough surface. Nonetheless, the development of the inner layer for the rough cases presents the same characteristics observed for the inner layer development over the smooth plate ( Figure 14A). Very close to the jet axis, that is, for r∕D < 0.2, the outer layer of the rough wall cases tends to be much more thin with respect to the smooth case. In the range 0.2 < x 3 ∕D < 1, the outer layer thickness appears to all the three investigated rough surfaces. The visible scatter of the IBM radial forcing should be expected since the volume force distribution at different (x 1 , x 2 ) locations will greatly differ in both intensity and wall-normal height distribution. On the other hand, the PFA volume force is applied evenly at each wall-normal height x 3 and the averaging process leads to a smooth curve. The case with the smallest roughness size (i.e., k 99 ∕D = 0.05) displays the best agreement between the PFA and IBM radial forcing distributions. It is also emphasized that this is the case that displays the best agreement for the mean velocity statistics. In the other two cases, the overall shape of the radial forcing distributions remains similar for both the PFA-and IBM-resolved solutions, even though the one estimated from PFA decays faster with increasing radial distance. This symptom suggests that additional fine tuning is required for the two cases. In fact, it is recalled that all the three PFA-resolved cases adopt the same model constants and the only difference lies in the different roughness shape functions used to spread the volume forces in the wall-normal direction. Although the two methods introduce very similar mean radial forcing for all the cases, the two methods can still produce different total wall-stress, as suggested by the analysis of the development of the inner and outer layers. The reason is the different mean flow predicted by the two methods. A control volume analysis on the mean flow field helps in identifying this issue.
Consider the control volume depicted in Figure 16 which is centered at a distance r from the jet axis and has width Δr and height k max , where k max indicates the maximum wall-normal height at which the IBM and PFA introduce nonzero volume forcing). The total wall-stress w can be obtained by applying the integral momentum balance to the considered control volume. More precisely, In the equation, S ABCD represents the area of the control volume. Note that the integral over the side AB of the control volume is not included in Equation (8) because fluxes at the bottom of the computational domain are, in general, zero. For small Δr, the term on the left-hand side of the equation can be approximated by w Δr, and therefore the total stress at a certain radial location can be easily deduced from Equation (8). Nonetheless, the considered momentum balance shows that it is not possible to assess the total wall-stress distribution by merely taking into account the mean radial forcing distribution. Momentum, pressure, and viscous stresses fluxes through the sides of the control volume are determined F I G U R E 16 Control volume adopted for the estimation of the total wall-stress at a fixed radial distance from the jet axis. A black dashed line depicts the sides of the control volume. The base of the control volume lies on the bottom of the computational domain, at x 3 = 0, while the top side of the volume is placed at x 3 ∕D = k max , where k max represents the maximum height at which nonzero volume force is applied by either the immersed boundary method or the parametric forcing approach by the mean flow field. Therefore, an accurate estimation of the latter is required in order to properly predict the total wall-stress distribution. This qualitative line of reasoning suggests that the observed discrepancies in the PFA-resolved mean flow field have important consequences on the estimation of a key quantity such as the total wall-stress.

CONCLUSIONS
In this work, the PFA was used to introduce the effects of surface roughness on the mean velocity field statistics of a turbulent impinging jet. The examined flow configuration offers the chance to test the novel PFA (presented in [6]) in the context of spatially developing turbulent flows over rough surfaces. After the jet impingement, a wall-jet develops along the radial direction. To the authors knowledge, it is the first time that the PFA is applied to such type of nonequilibrium flows.
To test the efficacy of the method, reference DNS using an IBM strategy to resolve the surface roughness were carried out. The presented results show that the PFA is capable of providing a qualitative good description of the mean flow features. In this respect, it was shown that for x 3 ∕D > k m the PFA predicted mean flow resembles closely the one predicted using the IBM. In particular, mean radial velocity profiles and boundary layer development are very similar for the two methods. Such good correspondence is quite remarkable given that no particular tuning was applied beforehand to the PFA model constants but they were derived from a turbulent channel flow configuration. This proves the robustness of the approach and encourages future developments of the PFA. In this respect, different areas of improvement of the PFA can be easily identified. The most immediate improvement of the method would be the inclusion of a wall-normal volume force. In this study it was shown that, in complex flows where regions may be dominated by strong wall-normal velocity component, the absence of wall-normal forcing has severe consequences in the estimation of important integral quantities, such as the total wall-stress. While the need of an additional wall-normal force term is straightforward to understand, the proper modeling of such component is challenging and it is the focus of ongoing work. The main issue to address in this respect is the introduction of a physically meaningful force model that can be tuned a priori based on the geometrical features of the surface roughness. The possibility of fully characterizing the PFA model by the knowledge of the surface geometry alone is another important issue that needs to be addressed in the future. The method adopted in this study uses two model constants, K k and C d , that need to be tuned beforehand in order to assure the proper intensity of the introduced volume forces. Ideally, it should be possible to fully characterize the model constants starting from the roughness topology alone. This idea was already suggested in the original version of the PFA presented in Busse and Sandham [1]. The resulting method would be completely determined by the roughness geometry alone, thus avoiding the burden of tuning any model constant.
Performance Computing Center Stuttgart (HLRS) under the grant number zzz44198. The authors would like gratefully to acknowledge the financial support by the German Research Foundation (DFG) for measuring equipment under HBFG programme INST 121384/178-1 FUGG. Open access funding enabled and organized by Projekt DEAL.