© N. van den Bos et al., Published by EDP Sciences, 2023

1. Introduction

A particularly important topic in the field of nuclear safety is the behaviour of fuel rods. Fuel rods are submerged in a coolant, such as water in current-day reactors or liquid metal in advanced fast reactors, which flows axially along the fuel rods. While the axial flow leads to efficient cooling of the fuel rods, it can also cause Turbulence-Induced Vibrations (TIV) of the fuel rods, which originate from pressure fluctuations in the flow. This phenomenon plays a critical role in terms of nuclear safety, as it can cause structural effects such as fatigue problems, fretting wear, and stress corrosion cracking [1, 2]. The phenomenon has been studied since the start of nuclear reactor development in the 1950s, and it has been the root of several incidents [1, 3]. In the previous century, an emphasis was put on experiments to increase the understanding and knowledge of TIV, as well as using semi-empirical laws to establish a relation between the amplitude of vibration and relevant flow parameters [4, 5]. More recently, due to the larger availability of computational resources, the use of Fluid-Structure Interaction (FSI) simulation applied to fuel rods has received more interest. In one approach the turbulence is (partially) scale resolved (LES/DNS), but due to the excessive computational requirements, these simulations typically rely on a simplified approach for the fluid-structure interaction, at moderate Reynolds numbers [6–8]. Due to the large Reynolds numbers associated with the axial flow along fuel rods, LES/DNS is too expensive for industrial use of simulations of fuel rods or fuel assemblies. For reference, the LES calculations of De Ridder [9] took 2200 days in equivalent computational time to simulate 0.1 s of 1/10th of the length of a fuel rod. Extending this to a fuel assembly is not realistic. Furthermore, these simulations are typically one-way coupled. However, due to the large density ratios that are typically found in nuclear reactors, two-way FSI coupling is preferred, in order to capture the effect of the changing flow field on the fuel rods. A computationally much cheaper method is to use URANS for the fluid domain [10, 11], which would allow simulations of fuel assemblies. However, URANS causes severe underprediction in the vibration amplitude. This is because URANS only calculates forces resulting from the mean flow, as it does not resolve the instantaneous turbulent pressure fluctuations. To remedy this, the Pressure Fluctuation Model (PFM) was proposed by Kottapalli et al. [12], which assumes that the turbulence statistics from URANS are sufficiently accurate to construct synthetic turbulent pressure fluctuations to be superimposed on the Reynolds-averaged pressure. This would artificially reintroduce the turbulence-induced forcing to the structure, at a lower computational cost than scale-resolving approaches. From comparison to experimental data, it was found that the PFM shows vibration amplitudes in the same order of magnitude, however, it did not give the required accuracy yet [12]. In the current paper, an improved pressure fluctuation model, called AniPFM, is proposed. The new model tackles several limiting assumptions of the PFM, such as the assumption of isotropic turbulence, as well as the method of time correlation. With this, the AniPFM the accuracy in the prediction of vibration amplitudes through FSI-simulations of TIV. Figure 1 summarizes how the AniPFM interacts with the other components in FSI simulations.

Fig. 1. A flow-chart of the global interactions between the several numerical solvers.

2. Methodology

The AniPFM is based on the formulation given of the PFM of Kottapalli et al. [12] and Senthooran et al. [13], though major modifications have been made to incorporate anisotropy and time correlation, using approaches from Billson et al. [14], and Shur et al. [15]. Details of the new AniPFM are given in this section.

2.1. Pressure fluctuations

Similarly, to the method of Kottapalli et al. [12], the divergence operator is applied to the Reynolds averaged and instantaneous Navier–Stokes equations, to obtain the following Poisson equation:

$$\begin{array}{c}\hfill \frac{{\partial }^2{p}^{\prime }}{\partial x_i\partial x_i}=-\rho [2\frac{\partial \overline{u_i}}{\partial x_j}\frac{\partial u_j^{\prime }}{\partial x_i}+\frac{{\partial }^2}{\partial x_i\partial x_j}(u_i^{\prime }{u}_j^{\prime }-\overline{u_i^{\prime }{u}_j^{\prime }})].\end{array}$$(1)

Equation (1) shows that the pressure fluctuation p′ solely depends on the mean velocity u̅′_i, the Reynolds stresses u̅′_iu̅′_j, and the velocity fluctuations u̅′_i. Only the velocity fluctuations are unknown. The governing equation for computing the pressure fluctuations is given by equation (1) where the velocity fluctuations are modelled.

Note that equation (1) assumes that the full spectrum of pressure fluctuations is computed, but in reality, only the fluctuations that can be resolved by the mesh are taken into account. Thus, the effect of unresolved turbulence on pressure fluctuations is not accounted for in this method.

2.2. Velocity fluctuations

Constructing the velocity fluctuations consists of three main parts, namely the construction of dimensionless velocity fluctuations, correlating the fluctuations in time, and finally scaling the fluctuations to match the input Reynolds stresses. The dimensionless fluctuations w_t (x) are created by a Fourier decomposition, as shown in equation (2). Here, q_n is the mode amplitude, σ_n is the direction vector, k_n is the wavenumber vector, and ϕ_n is a random phase shift with a uniform distribution, the subscript n denotes that it is for the nth mode.

$$\begin{array}{c}\hfill {\mathbf{w}}_t\left(\mathbf{x}\right)=\sqrt{6}\sum _n^N\sqrt{q_n}\left[{\sigma }_{n}cos\left({\kappa }_n·\mathbf{x}+{\varphi }_n\right)\right].\end{array}$$(2)

2.2.1. Turbulent kinetic energy spectrum

The amplitude for each mode is determined by the non-dimensionalized Von-Karman spectrum, shown in equation (3) [15]. Here, k_e is the eddy wavenumber, which is the wavenumber which contains the highest energy density. The variable k_η refers to the Kolmogorov wavenumber, and f_cut is a cut-off frequency based on the maximum resolvable wavenumber. As shown in equation (4), the amplitude is based on the relative weight of the energy of mode n. Note that even though the AniPFM can represent anisotropic turbulence, the energy spectrum is based on isotropic turbulence, thus it is still expected that the results near walls are not as accurate as in fully isotropic conditions.

$$\begin{array}{c}\hfill E\left(k\right)=\frac{{\left(k/k_e\right)}^4}{{\left[1+2.4{\left(k/k_e\right)}^2\right]}^{\left(17/6\right)}}exp\left(-{\left(12\frac{k}{k_{\eta }}\right)}^2\right)f_{}\text{cut}\end{array}$$(3)

$$\begin{array}{c}\hfill q_n=\frac{E_k\left(k_n\right)\mathrm{\Delta }{k}_n}{{\sum }_n^N{E}_k\left(k_n\right)\mathrm{\Delta }{k}_n}.\end{array}$$(4)

The eddy wavenumber is based on the RANS length scale as shown in equation (5), where C_l is a calibration coefficient set to 3.0, based on isotropic turbulence.

$$\begin{array}{c}\hfill l_e=C_l·l_t.\end{array}$$(5)

The cut-off function is shown in equation (6) [15]. For the cut-off length, several definitions were investigated. However, it was found that the cut-off length as defined by Shur et al. [15] (shown in Eq. (7)) gives the most accurate results.

$$\begin{array}{c}\hfill f_{}\text{cut}=exp\left(-{\left[\frac{4max\left(k-0.9k_{}\text{cut},0\right)}{k_{}\text{cut}}\right]}^3\right)\end{array}$$(6)

$$\begin{array}{c}\hfill l_{}\text{cut}=2min\left[max\left(h_y,h_z,0.3h_{}\text{max}\right)+0.1d_w,h_{}\text{max}\right].\end{array}$$(7)

2.2.2. Wavenumber and direction vector

As shown in equation (2), the velocity fluctuations are determined by a sum of wavenumbers. As the sum is finite, it requires the specification of a start and final wavenumber, as well as the distribution of discrete wavenumbers in the resolved wavenumber space. The starting wavenumber is based on a conservative dimensional analysis, shown in equation (8). In the case of a zero-mean velocity, the RANS length scale is used to determine the starting wavenumber. Additionally, the starting wavenumber must also be lower than the wavenumber with the maximum energy density. Finally, the geometrical considerations must also be taken into account. This is done by introducing a user length scale, which denotes the maximum length that can be captured given the geometry of the problem. These requirements are fulfilled by the expression for the starting wavenumber, as shown in equation (9).

$$\begin{array}{c}\hfill k_{}\text{start}=\frac{\epsilon }{max\left({\left|\left|u\right|\right|}^3\right)}\end{array}$$(8)

$$\begin{array}{c}\hfill k_{}\text{start}=max\left[min\left(k_{}\text{start},\frac{1}{2}{k}_e\right),\frac{2\pi }{l_{}\text{user}}\right].\end{array}$$(9)

The highest wavenumber is based on the cut-off wavenumber (see Eq. (10)), as the energy spectrum quickly goes to zero after this point. In between the lowest and highest wavenumber, the wavenumbers are logarithmically distributed.

$$\begin{array}{c}\hfill k_{}\text{end}=\frac{3}{2}{k}_{}\text{cut}.\end{array}$$(10)

The wavenumber vectors are defined similarly to [12, 14], as shown in Figure 2. The mathematical description is displayed in equation (11). Here, θ_n and ψ_n are uniformly distributed random variables, described by equation (12).

Fig. 2. Wave vector geometry of the th Fourier mode [14].

$$\begin{array}{c}\hfill {\kappa }_{\mathbf{n}}={\kappa }_{\mathbf{n}}{\left[sin{\theta }_{\mathbf{n}}cos{\psi }_{\mathbf{n}},sin{\theta }_{\mathbf{n}}sin{\psi }_{\mathbf{n}},cos{\theta }_{\mathbf{n}}\right]}^{\mathbf{T}}\end{array}$$(11)

$$\begin{array}{c}\hfill P\left({\psi }_n\right)=\frac{1}{\pi },P\left({\varphi }_n\right)=\frac{1}{2\pi },P\left({\theta }_n\right)=\frac{1}{2}sin\left({\theta }_n\right).\end{array}$$(12)

The wavenumber direction vector is determined from the wavenumber vector. From applying the divergence operator to equation (1), it can be found that k_n · σn = 0 must hold in order to adhere to the continuity equation. This is the only requirement. To achieve this, σn is defined as the normalised cross-product between a random vector ζ_n, and the wavenumber vector, as shown in equation (13). Since the dot-product of a vector with the cross product of the same vector is always equal to zero, the continuity condition is met.

$$\begin{array}{c}\hfill {\sigma }_n=\frac{{\zeta }_{\mathbf{n}}\times {\mathbf{k}}_{\mathbf{n}}}{\left|{\zeta }_{\mathbf{n}}\times {\mathbf{k}}_{\mathbf{n}}\right|}.\end{array}$$(13)

2.2.3. Time correlation and scaling

The procedure for calculating the non-dimensional fluctuations w_t (x) has been fully defined in the previous subsections. The next step is to create a space-time-dependent velocity signal v_t (x, t), from the purely space-dependent signal. Two main phenomena have to be taken into account when constructing a time-correlated velocity field. Namely, the convection of turbulent eddies and the decorrelation due to the production and dissipation terms.

A method similar to Billson et al. [14] is used to introduce this time dependency. First, the fluctuations are convected by solving for equation (14). Here, v_t^m−1 are the non-dimensional velocity fluctuations generated at timestep m − 1, and U_j is the Reynolds-averaged velocity as calculated by the accompanied URANS simulation. Then in the second part, a new solution v^m_t is calculated from a combination of the (convected) previous solution v_t^m–1, and a newly generated field w^m_t, as shown in equation (15).

$$\begin{array}{c}\hfill \frac{\partial {\mathbf{v}}_t^{m-1}}{\partial t}+U_j\frac{\partial {\mathbf{v}}_t^{m-1}}{\partial x_j}=0\end{array}$$(14)

$$\begin{array}{c}\hfill {\mathbf{v}}_{\mathbf{t}}^{\mathbf{m}}\left(\mathbf{x},\mathbf{t}\right)=\mathbf{a}{\mathbf{v}}_{\mathbf{t}}^{\mathbf{m}-\mathbf{1}}\left(\mathbf{x}\right)+\mathbf{b}{\mathbf{w}}_{\mathbf{t}}^{\mathbf{m}}\left(\mathbf{x}\right).\end{array}$$(15)

The coefficients a and b are defined in equation (16), where τ is the dissipation timescale determined from the URANS simulation, and f_τ is a modification factor for fine-tuning the correlation. For the latter, a value of f_τ = 17 is used, in line with what was used in [14] for the simulation of a 3D jet. Tests showed that this value also gives a satisfactory correlation in the simulations of interest. The coefficients a and b are defined such that the squared mean properties of v_t are still respected, i.e. v²_t = δ_ij.

$$\begin{array}{c}\hfill a=e^{-f_{\tau }\mathrm{\Delta }t/\tau }b=\sqrt{1-a^2.}\end{array}$$(16)

Finally, the last step is to scale the dimensionless fluctuations, such that they replicate the Reynolds stresses. Since 〈v_t〉² is equal to the Kronecker delta, the fluctuations need to be scaled with a tensor a_ij such that a_ija_ji = R_ij, as shown in equation (17). This can be realised with a Cholesky decomposition of the Reynolds stress tensor, which is shown in equation (18).

$$\begin{array}{c}\hfill {\mathbf{u}}_t(\mathbf{x},t)=a_{\text{ij}}{\mathbf{v}}_t(\mathbf{x},t)\end{array}$$(17)

$$\begin{array}{c}\hfill a_{\mathit{ij}}=\left[\begin{array}{ccc}\sqrt{R_{11}}& 0& 0\\ \frac{R_{21}}{a_{11}}& \sqrt{R_{22}-a_{21}^2}& 0\\ \frac{R_{31}}{a_{11}}& \frac{(R_{32}-a_{21}{a}_{31})}{a_{22}}& \sqrt{R_{33}-a_{31}^2-a_{32}^2}\\ \end{array}\right].\end{array}$$(18)

With this method, the AniPFM can reconstruct anisotropic Reynolds stresses. For flows with a constant pressure gradient, such as channel flows, linear Eddy viscosity models show isotropic Reynolds stresses. In order to improve the accuracy of these models, a correction is used to transform the isotropic tensor into an anisotropic tensor, based on the nonlinear Eddy viscosity model of Wilcox [16]. The Wilcox correction is shown in equation (19).

$$\begin{array}{c}\hfill \overline{u^{\prime }{u}^{\prime }}=\frac{8}{9}k,\overline{v^{\prime }{v}^{\prime }}=\frac{4}{9}k,\overline{w^{\prime }{w}^{\prime }}=\frac{6}{9}k.\end{array}$$(19)

2.3. Model overview

This concludes the full model of the AniPFM. In Figure 3, the full model is summarised. First, the non-dimensional velocity fluctuations are calculated, based on the energy spectrum. Then the time correlation is performed. After this, the velocity fluctuations are computed by scaling with the Cholesky tensor. Finally, p′ can be computed.

Fig. 3. Flow chart of the different computational steps of the proposed AniPFM.

3. Simulation set-up

Two cases were simulated, in order to validate the AniPFM. First, a Homogeneous Isotropic Turbulent (HIT) box was simulated, which was used to verify the AniPFM in isotropic conditions, since an isotropic energy spectrum was used as input. Second, a turbulent channel flow was simulated, with the purpose of testing the accuracy of the AniPFM versus the PFM in anisotropic turbulence. The setup of these two cases is discussed in this section. The AniPFM is implemented in OpenFOAM 8, and thus all subsequent simulations have been performed with the same software.

3.1. Homogeneous isotropic turbulent box

A box of L · L · L is created, which is discretized by a cartesian mesh of N · N · N cells. All boundaries are periodic, and the domain has a zero-mean velocity. The experiment of Comte-Bellot and Corrsin [17] has been replicated, as well as the DNS of Gotoh and Fukayama [18]. In Table 1, the details of the replicated set-up of Comte-Bellot and Corrsin are shown. The DNS of Gotoh and Fukayama was specified in dimensionless numbers. The simulation with Re_λ = 287 of Gotoh and Fukayama is replicated; the set-up is also shown in Table 1.

3.2. Turbulent channel flow

The results of the TCF are compared to the DNS results of Abe et al. [19]. The simulation with Re_τ = 640 is used for comparison. In this case, the flow between two infinitely long and wide stationary plates is simulated. The simulation domain used is equal to 6δ · 2δ · 3δ, where δ is the mid-channel height. The mesh has a size of N_x · N_y · N_z cells, which are kept as variables. The mesh distribution is uniform in the stream- and spanwise-direction, and it is geometrically expanding to the mid-channel plane in the wall-normal direction. The mesh in the wall-normal direction is configured such that y⁺ ≈ 1 for the first grid cell from the wall. The setup is summarised in Table 2, and an example mesh is shown in Figure 4. Simulations were performed with the k − w SST turbulence model, along with the Wilcox correction. However, for several simulations, no URANS calculation was done, but rather the mean flow properties of the DNS of Abe et al. were used as an input to the AniPFM. This was used to isolate any errors that originate due to the AniPFM, and not due to the input.

Fig. 4. Example channel flow mesh of 40 · 64 · 30.

4. Results discussion

The results of both cases are elaborated upon in this section. In Table 3, the V&V done for all variables of the test cases is shown. In this section, the most important results from this verification are shown and discussed.

4.1. Homogeneous isotropic turbulent box

In Figure 5, the simulated velocity fluctuations of the HIT box of Comte-Bellot Corrsin are shown, with a mesh size of N = 128. The velocity fluctuations on the left are from the AniPFM, whereas the results on the right are simulated through means of a Large-Eddy Simulation [20]. From Figure 5, it was concluded that the AniPFM can qualitatively reconstruct similar flow structures as high-resolution methods.

Fig. 5. The instantaneous velocity fluctuations of the AniPFM, compared to LES results.

In Figure 6, the energy spectra generated by both the AniPFM and the PFM of Kottapalli et al. are compared to the experimental spectrum of Comte-Bellot & Corrsin. The AniPFM results show a very good resemblance with respect to the experimental results, right up to the cut-off wavenumber. Compared to the results of the PFM, the peak of the energy spectrum is better predicted, and there is no unphysical accumulation of energy near the cut-off wavenumber. The latter difference is due to the fact that a cut-off filter is used in the AniPFM, which makes sure that the unresolved energy is not redistributed over the resolved part of the spectrum.

Fig. 6. The energy spectrum of AniPFM (left) and the PFM (right), compared to the experimental values of Comte-Bellot and Corrsin [17].

The effect of the cut-off filter can also be found in the root-mean-squared pressure fluctuations. With a mesh of N = 64, the AniPFM predicted p′_rms within a 2.4% error with respect to the results of Goth et al., whereas the PFM predicted a 3.7% error. The AniPFM showed an error of 1.1% when the mesh was refined N = 128.

4.2. Turbulent channel flow

For the turbulent channel flow, there are several sources of errors that cause a discrepancy between the AniPFM results and experimental/DNS data. Therefore, it is important that each source is carefully evaluated.

There exists an uncertainty range in the results of the AniPFM, which is due to the random numbers that are used throughout the simulation. For the chosen temporal correlation scheme, the random numbers gave an uncertainty range of ±0.25% for p′_rms, for a confidence interval of 95.4%. The number of modes did not seem to have a large effect on p′_rms, however, with a lower number of modes, the uncertainty went up, as fewer random numbers were used per iteration. For the presented simulations, 1024 wave number modes were used.

The amount of energy that is resolved by the AniPFM depends on the fineness of the mesh. There is no sub-grid model that models the effect of unresolved velocity fluctuations on the pressure fluctuations. Thus, the statistics of the pressure fluctuations only converge if a large part of the velocity fluctuations is resolved. Since p′_rms is the quantity of interest, this quantity is evaluated for several meshes. This is shown in Figure 6, in which the results were obtained using the k − w SST turbulence model, together with the Wilcox correction. It can be seen that the meshes converge to a final solution when using finer meshes. At the second finest mesh, the results near the wall are within the given uncertainty range, thus deeming the solution converged.

The modelling error was found by performing a simulation of the turbulent channel flow on the finest mesh from Figure 7, but with the RANS input variables (k, ε, u̅_i and u̅′_iu̅′_j taken from DNS data. The replicated Reynolds stresses and p′_rms are shown in Figure 8. The Reynolds stresses are very closely approximated. The small underestimation is because the velocity fluctuations are not fully resolved. Nevertheless, it was found that these unresolved fluctuations had no effect on p′_rms. From Figure 8 (right), it can be seen that near the midchannel plane, the AniPFM very closely approximates p′_rms. This is because here the isotropic energy spectrum very closely approximates the actual energy spectrum. Near the wall, the energy spectrum is not approximated as accurately, due to the larger anisotropy in the flow. Thus, a larger discrepancy in p′_rms was found near the wall, with a maximum error of 4.4%.

Fig. 7. The root-mean-squared pressure fluctuations along the wall-normal coordinate for various meshes, versus the DNS results of Abe et al. [19].

Fig. 8. The Reynolds stress profiles (left) and the RMS pressure fluctuations (right) along the wall-normal coordinate, versus the DNS results of Abe et al. [19]. DNS data is used as input.

When the k − w SST turbulence model, together with the Wilcox correction was used as URANS input for the turbulent channel flow, an error of roughly 10% was found near the wall. This increase in error is due to the large difference in the Reynolds stress tensor between the selected RANS model and the DNS data. In comparison, the PFM showed an underestimation of roughly 47% at the wall. Thus, the AniPFM shows a sharp improvement in the prediction of p′_rms for anisotropic flows, with respect to the PFM.

5. Conclusion

In this paper, a new pressure fluctuation model was proposed, called AniPFM. This model improves the prediction of pressure fluctuations when using a URANS approach, which can be useful in particular for turbulence-induced vibration prediction. Several aspects of the AniPFM were adjusted with respect to the previous PFM of Kottapalli et al., namely the energy spectrum cut-off filter, the replication of anisotropic Reynolds stresses, and the method for time correlation.

Two test cases were performed to assess the performance of the AniPFM: namely, a homogeneous isotropic turbulent box, and a turbulent channel flow. From the HIT case, it was found that the prediction for p′_rms was slightly improved due to the addition of the cut-off filter. From the turbulent channel flow, it was found that the introduction of anisotropy gave a slightly higher model error. This is because the energy spectrum is still based on isotropic turbulence. Nevertheless, the AniPFM showed an improvement in the prediction of p′_rms at the wall, compared to the PFM of Kottapalli et al.

While the initial results give an optimistic view, there are still several challenges that need to be tackled. Further validation must be done, both for fluid-only cases, as well as for relevant use cases. While there is no consensus if the AniPFM two-way coupled approach would be an improvement with respect to the one-way coupled LES approach for determining fuel pin vibration, the research does show to be of interest for strongly coupled cases, such as full fuel assembly FSI simulations. For this, both the computational costs and accuracy of the AniPFM must be compared to LES or Hybrid URANS/LES FSI simulations, to evaluate if the AniPFM is beneficial to use. Finally, there are still assumptions in the AniPFM which have to be evaluated, such as the use of an isotropic energy spectrum scaled to meet the Reynolds stress requirements. The construction of this spectrum is based on integral properties of the flow, with assumed isotropy, meaning that the energy distribution might not be well represented in areas of high anisotropy.

The next step is to apply the model to fluid flow cases that resemble nuclear fuel rod operating conditions more closely. While validation is still ongoing, the AniPFM shows potential for accurate, computationally cheap simulations of turbulence-induced vibrations in industrial nuclear applications.

Conflict of interests

The authors declare that they have no competing interests to report.

Acknowledgments

The authors would like to acknowledge the contribution of the PIONEER research program funded by the ministry of economic affairs at the Nuclear Research & consultancy Group and the Technological University of Delft for supporting the research project and providing the necessary resources. The authors would also like to thank Dr. Tajfirooz and Dr. Dovizio for their feedback and helpful discussions on the AniPFM implementation.

Funding

The work described in this paper is funded by the Dutch national research program PIONEER sponsored by the ministry of economic affairs.

Data availability statement

The generated and analysed data associated with this paper is available, and it can be requested from the author K. Zwijsen.

Author contribution statement

All the authors were involved in the preparation of the article. All the authors have read and approved the final article.

References

All Tables

All Figures

	Fig. 1. A flow-chart of the global interactions between the several numerical solvers.
In the text

	Fig. 2. Wave vector geometry of the th Fourier mode [14].
In the text

	Fig. 3. Flow chart of the different computational steps of the proposed AniPFM.
In the text

	Fig. 4. Example channel flow mesh of 40 · 64 · 30.
In the text

	Fig. 5. The instantaneous velocity fluctuations of the AniPFM, compared to LES results.
In the text

	Fig. 6. The energy spectrum of AniPFM (left) and the PFM (right), compared to the experimental values of Comte-Bellot and Corrsin [17].
In the text

	Fig. 7. The root-mean-squared pressure fluctuations along the wall-normal coordinate for various meshes, versus the DNS results of Abe et al. [19].
In the text

	Fig. 8. The Reynolds stress profiles (left) and the RMS pressure fluctuations (right) along the wall-normal coordinate, versus the DNS results of Abe et al. [19]. DNS data is used as input.
In the text