Turbulence anisotropy effects on corner-flow separation: physics and turbulence modelling

Abstract The secondary motion caused by turbulence anisotropy is one of the crucial factors for determining the size of corner-flow separation in a side-wall interference flow field. Therefore, through a wall-resolved large-eddy simulation (LES) of a side-wall interference flow field, this study investigates the effects of the secondary motion on the corner-flow separation and explores the turbulence modelling that can reproduce the secondary flow motion. The momentum transport analysis using the LES results shows that the secondary vortex has twofold effects on delaying the corner-flow separation: the convective transport of the streamwise momentum towards the corner, and the enhanced production of turbulence by increasing the shear. Also, the vorticity transport analysis reconfirms that the secondary motion is caused primarily by turbulence anisotropy in the outer layer of the turbulent boundary layer. Furthermore, a quadratic constitutive relation (QCR) is proposed based on the analysis of the relationship between the Reynolds stress and velocity gradient. The proposed QCR consists of two quadratic terms and three constant parameters. The a priori analysis using the LES data shows that the proposed QCR represents the anisotropy of the Reynolds stress overall better than the existing QCR. Reynolds-averaged Navier–Stokes simulation using the proposed QCR with the Spalart–Allmaras turbulence model shows improvements in the prediction of the corner-flow separation compared to the results obtained by the existing QCR with the same turbulence model.


Introduction
Accurate prediction of the separation pattern at a stall condition is crucial in ensuring the safety of aircraft because the separation causes loss of local lift, which potentially induces harmful pitch-up motion.One of the typical separation patterns on aircraft is corner-flow separation at the wing-body junction.Near the junction, the boundary layers over the wing and fuselage interact, causing a low-velocity region prone to separation.Therefore, substantial efforts have been made to understand and predict wing-body juncture flows, as summarized by Gand et al. (2010).More recently, the NASA Juncture Flow experiment (Rumsey, Neuhart & Kegerise 2016;Rumsey et al. 2022) was conducted to obtain validation data for computational simulations.They employed a simplified wing-body geometry, and measured velocity and Reynolds stress profiles in the corner region by laser Doppler velocimetry.Subsequent studies (Abe, Mizobuchi & Matsuo 2020;Abdol-Hamid et al. 2020;Eisfeld et al. 2022;Ghate et al. 2020;Iyer & Malik 2020;Lozano-Duran, Bose & Moin 2020;Rumsey et al. 2020) have utilized the data from this experiment to validate and improve simulation methodologies and turbulence models.
One of the crucial features of wing-body juncture flows is the interaction between the boundary layers over the wing and body.The side-wall interaction causes the second type of Prandtl's secondary flow (Bradshaw 1987), i.e. a flow motion caused by the anisotropy of the Reynolds normal stress.The secondary flow motion appears as a pair of longitudinal (i.e.streamwise) vortices at the corner.A typical example of this type of flow is the flow in a square duct, which has been investigated by many researchers (e.g.Pinelli et al. 2010;Vinuesa et al. 2014;Zhang et al. 2015;Pirozzoli et al. 2018).Most recently, Pirozzoli et al. (2018) conducted direct numerical simulations (DNS) of square duct flows at several Reynolds numbers up to the friction Reynolds number Re τ = 1.0 × 10 3 , and discussed the transport of velocity and vorticity.They showed that the secondary motion transports streamwise momentum towards the duct corner.Hence the secondary motion essentially delays the corner-flow separation.However, the effects of the secondary vortices on the corner-flow separation have rarely been discussed quantitatively in previous studies.The challenges to this kind of discussion are that square duct flows do not involve corner-flow separation, whereas realistic wing-body juncture flows are too complex for quantitative analysis of the momentum transport by the secondary motion.
Moreover, simulating corner-flow separation using turbulence models based on the Reynolds-averaged Navier-Stokes (RANS) equations also poses challenges.In particular, the turbulence model must reproduce accurately the anisotropy of Reynolds stress, which cannot be realized by classical turbulence models based on the Boussinesq approximation.Hence several types of turbulence models have been proposed to incorporate the effects of anisotropy.One approach is the Reynolds stress transport models (e.g.Launder, Reece & Rodi 1975;Speziale, Sarkar & Gatski 1991), which involve six partial differential equations to account for the six independent components of Reynolds stress.To reduce computational costs, several researchers (e.g.Gatski & Speziale 1993;Wallin & Johansson 2000) have also proposed algebraic stress models.In these models, the Reynolds stress equations are reduced to algebraic equations.Another approach, which is simpler than the models above, is the model based on the nonlinear expression of Reynolds stress (e.g.Speziale 1987;Spalart 2000).One of the most prevalent models in this approach is the quadratic constitutive relation (QCR) proposed by Spalart (2000) (often referred to as QCR2000).This model provides each Reynolds stress component by a simple algebraic equation consisting only of scalar eddy viscosity, strain, vorticity, and one constant parameter.More recently, modifications to QCR2000 have also been proposed (Mani et al. 2013;Rumsey et al. 2020).Mani et al. (2013) introduced a modification to QCR2000 that incorporates the effects of turbulence kinetic energy, which are not included in the original Spalart-Allmaras (SA) turbulence model.Additionally, Rumsey et al. (2020) investigated the validity of the QCR based on data from the NASA Juncture Flow experiment, and presented a modified version of the model where the model parameters vary in space.980 A21-2 Several existing studies have reported that the QCR improves the prediction of the corner-flow separation at the wing-body junction.Yamamoto, Tanaka & Murayama (2012) found that the SA turbulence model with QCR2000 shows the pitch-up motion of the DLR-F6 aircraft model at the high angle of attack.This trend is in good agreement with the past wind tunnel experiment (Rivers & Dittberner 2011), while the baseline SA model (without QCR2000) does not predict it.The remarkable difference between the two simulation results is the side-body separation (i.e. the corner-flow separation at the wing-body junction).Here, QCR2000 suppresses the side-body separation, resulting in the correct pitch motion at high angles of attack.Yamamoto et al. (2012) explained that the difference in the side-body separation occurs because the secondary motion convects momentum to the corner region.Dandois (2014) also reported a similar improvement in predicting the corner-flow separation of the DLR-F6 aircraft model.After these studies, QCR2000 has been used widely in RANS simulations of aircraft at high angles of attack (e.g.Rumsey 2018;Tamaki & Imamura 2018;Tinoco et al. 2018).
However, despite the prevalence of the model, the validity of the parameter of QCR2000 (referred to as C cr1 in the original paper) remains unclear.Spalart (2000) showed qualitative improvement in predicting the secondary motion of a square duct flow compared to the baseline turbulence model using the Boussinesq approximation.The square duct flow and the wing-body junction flow differ in many points, such as the existence of the corner-flow separation or the streamwise development of the boundary layer.Furthermore, although modifications to QCR2000 (Mani et al. 2013;Rumsey et al. 2020) have been proposed, their validation has been limited to comparisons to the measurements in wind tunnel experiments.Therefore, comprehensive validation using a database from high Reynolds number DNS or large-eddy simulation (LES) will be beneficial for a better understanding and improvement of the QCR.
Considering the situations stated in the preceding paragraphs, we conduct a LES of a simplified side-wall interference flow field imitating the wing-body junction of aircraft.The purpose of this simulation is twofold.The first purpose is to explain the effects of the secondary motion on the corner-flow separation based on the momentum budget analysis.Unlike the former studies on square duct or wing-body junction flows, the geometry simulated in this study involves corner-flow separation while retaining a simple geometry suitable for the budget analysis.The second purpose is to provide databases for turbulence modelling.(The obtained data are uploaded on the website https://www.klab.mech.tohoku.ac.jp/database.)Using the obtained database, we investigate the validity of the constitutive relations and propose a QCR that can reproduce the Reynolds stress distributions better.Also, to eliminate undesirable low Reynolds number effects, we conduct the simulation at a relatively high Reynolds number Re L ∼ 10 6 , where Re L is the Reynolds number based on the freestream velocity and the length from the leading edge.This Reynolds number is comparable to that of large-scale wind tunnel testing.Pirozzoli et al. (2018) have reported in their study on the square duct flow that the strength of the secondary motion saturates at high Reynolds numbers.Since a LES at a high Reynolds number requires massive computational resources, the database may be obtained first by the state-of-the-art supercomputer.
The structure of this paper is as follows.In § 2, the computational setting and methodologies of the LES are presented.Then § 3 describes the results from the LES and analyses based on the momentum and vorticity budgets, which will explain the effects and generation mechanism of the secondary motion.In § 4, we explore the modelling of the constitutive relation using the LES data, and propose a modified QCR.The developed model is validated in RANS simulations in § 5, and finally, § 6 provides the concluding remarks.

Geometry
We simulate the side-wall interference flow field shown in figure 1.The computational domain has a square cross-section whose size increases and decreases in the streamwise direction.The boundary layers on the bottom and side walls interfere, and the change in the cross-section induces a pressure gradient in the streamwise direction.The shape of the walls is designed so that a small separation bubble occurs at the corner of the diverging and converging sections.Here, the bottom wall coordinate y w is defined as where y 0 = −0.16L,H = 0.03L, x 0 = W = L (see figure 2a), and L is the length between the leading edge and the starting position of the change in the wall geometry.These parameters were determined through preliminary RANS analysis so that corner-flow separation occurs while the boundary layer away from the side wall does not separate.The wall starts from x = 0, and the laminar boundary layer is tripped at x/L = 0.1 by introducing the unsteady body force presented by Schlatter & Örlü (2012).The parameters for the tripping, such as the temporal and spatial cut-off scales, are the same as the baseline case in Schlatter & Örlü (2012), where we employ the laminar displacement thickness at the tripping location δ * 0 /L = 6.0 × 10 −4 .The side wall is set so that the geometry is symmetric about y = z.Also, a periodic boundary condition is implemented so that velocity vector (u, v, w) and pressure p at y/L = 0 match (u, −w, v) and p at z/L = 0, respectively.The Reynolds number based on the freestream velocity u ∞ and length L, i.e.Re L ≡ u ∞ L/ν ∞ , is set to 1.0 × 10 6 , where ν ∞ is the freestream kinematic viscosity.At x/L = 1.0, the Reynolds number Re L = 1.0 × 10 6 corresponds to Re θ ≡ u ∞ θ/ν ∞ ≈ 2.0 × 10 3 and Re τ ≡ u τ δ 99 /ν w ≈ 6.9 × 10 2 , where θ is the momentum thickness, δ 99 is the 99 % boundary layer thickness, u τ is the friction velocity, and ν w is the kinematic viscosity at the wall.

Computational grid
The computational grid is a structured curvilinear grid.As shown in figure 2, the grid within the y-z plane is orthogonal, while the streamwise grid lines are not always perpendicular to the y-z plane.In the turbulent flow region (x/L > 0.1), the grid spacing in the streamwise direction is set to 3.0 × 10 −4 L to 4.0 × 10 −4 L. In x/L < 1.0, the spanwise grid spacing, except for the near-wall region, is set to 2.0 × 10 −4 L. In the diverging and converging sections (1.0 < x/L < 2.0), the near-wall grid resolution is retained, while the outer grid in the y-z plane is stretched slightly as shown in figure 2(b).These grid spacings are determined based on the viscous wall unit at (x/L, y/L) = (1.0,0.0), i.e. a zero-pressure gradient location sufficiently away from the side wall.At this location, the streamwise and spanwise grid spacings in the wall viscous unit are approximately 16 and 8, respectively.Also, the wall-normal grid spacing at the wall is 1.5 × 10 −5 L, which is determined so that the spacing in the wall viscous unit does not exceed 0.9.In the diverging and converging sections, the grid spacing in the wall viscous unit becomes less than that described above because the skin friction decreases.

Numerical methods
The present simulation is based on the spatially filtered compressible Navier-Stokes equations.The fluid is an ideal gas with molecular viscosity following Sutherland's law.Here, we employ an implicit LES regime that has been well validated in our previous works (e.g.Kawai, Shankar & Lele 2010;Asada & Kawai 2018;Tamaki & Kawai 2023).The space is discretized based on the finite-difference method, where the spatial derivatives are evaluated by the sixth-order compact difference scheme (Lele 1992).Also, we introduce the tridiagonal compact filter (Lele 1992;Gaitonde & Visbal 2000) to remove high-frequency numerical oscillations.Since the filter plays a similar role to a subgrid-scale (SGS) turbulence model, the simulation does not employ any explicit SGS turbulence model.
For the time integration, we employ a second-order implicit scheme (Obayashi, Fujii & Gavali 1988;Izuka 2006) with five subiterations per time step.The magnitude of the residual reduces by more than two orders during the subiterations.Here, the time increment is set to t u ∞ /L = 2.4 × 10 −5 , corresponding to a viscous-scale time increment t + ∼ 0.1 and maximum Courant number approximately 8. According to the previous studies (Choi & Moin 1994;Kawai et al. 2010), this time increment ( t + ∼ 0.1) is sufficiently small for simulating wall turbulence.
The present LES is conducted on the Fugaku supercomputer provided by the RIKEN Center of Computational Science.The computational domain is divided into 60 × 8 × 8 subdomains, where the message-passing interface is used for inter-subdomain communications.Here, four subdomains are allocated per computational node (48 cores A64FX CPU).Furthermore, the computation for each subdomain is parallelized by OpenMP with 12 threads.

Statistical averaging
For statistical averaging, we conduct the simulation during a period T ave u ∞ /L = 144 (6 × 10 6 time steps) after an initial period Tu ∞ /L = 33.6,where T is the time.Previous studies on square duct flows (Pinelli et al. 2010;Vinuesa et al. 2014;Pirozzoli et al. 2018) reported that the statistical convergence of the secondary flow motion takes far longer than the convective time unit.In these studies, the statistical average is taken for T ave U b /h from 5.9 × 10 3 (Vinuesa et al. 2014) to 1.0 × 10 4 (Pirozzoli et al. 2018), where U b is the bulk velocity, and h is the half-height of the duct.Since the present simulation consists of developing boundary layers, we use the nominal freestream velocity u ∞ and the boundary layer thickness x/L = 1.0 (δ 99 /L = 0.0164; see § 3.1) for the normalization.Essentially, these values have the same meaning as U b and h in the duct-flow cases.The averaging period using these scales is calculated as T ave u ∞ /δ 99 = 8.8 × 10 3 , which has the same order of magnitude as values in the studies described above.
Also, when showing cross-sectional distributions of variables, we take averaging through the range ±0.01L in the x direction and by inverting around the y = z line.Appendix A summarizes the statistical convergence and the effects of the spatial averaging.

Overview
To validate the employed methodologies and computational grid, we first examine the velocity and Reynolds stress profiles of the fully developed turbulent boundary layer at x/L = 1.0,where the pressure gradient is almost negligible.Here, the obtained data are averaged over z/L ∈ [−0.01, 0.0], which is sufficiently away from the side wall.At this location, δ 99 /L = 0.0169 and θ/L = 0.00204 (i.e.Re θ = 2.04 × 10 3 ).Figure 3 compares the obtained mean streamwise velocity and Reynolds stress profiles in the wall viscous units to the DNS data at a similar Reynolds number (Re θ = 2.00 × 10 3 ) presented by Schlatter & Örlü (2010).Here, the overline and prime denote averaged and fluctuating   components, respectively.In the results, both streamwise velocity and Reynolds stress components are in excellent agreement with the DNS data.These results show that the spatial and temporal resolutions of the present LES are sufficiently high as wall-resolved LES and closer to DNS.Note that the density fluctuation is less than 1 % of the freestream value, suggesting that the flow is almost incompressible.
Figure 4 shows the distributions of ū, mean pressure coefficient C p , and mean skin friction coefficient C f over the periodic boundary plane (y/L = 0.0).As shown here, the flow decelerates by the expansion of the cross-section in 1.0 < x/L < 1.5, where an adverse pressure gradient occurs.Since the pressure gradient is mild, C f remains positive, i.e. the mean flow separation does not occur at the locations away from the side wall, though the flow involves the corner-flow separation as described in the following paragraph.Furthermore, C f in x/L < 1.0 shows reasonable agreement with the classical power-law flat-plate correlation C f = 0.027 Re −1/7 x (White 2006), where Re x ≡ u ∞ x/ν ∞ .Next, we examine the entire flow field, including the effects of the side wall.Figure 5 shows the distributions of ū, turbulence kinetic energy (TKE) K ≡ (u u + v v + w w )/2, and mean streamwise vorticity ω x ≡ ∂ w/∂y − ∂ v/∂z at several streamwise sections.x/L = 1.4 x/L = 1.5 x/L = 1.6 x/L = 1.4 x/L = 1.5 x/L = 1.6 x/L = 1.4 x/L = 1.5 The streamwise velocity distributions indicate that flow separation occurs at the corner near x/L = 1.4 and reattaches at x/L ≈ 1.6.In the diverging section (1.0 < x/L < 1.5), the TKE distributions (see figure 5b) show a significant increase.Also, in figure 5(c), a vortex pair is observed at the corner of each section.The size of the vortices increases as the flow field diverges and is distributed away from the corner.
Figure 6 shows the mean velocity distributions near the corner at streamwise cross-sections x/L = 1.0, 1.5 and 2.0.The contours of the streamwise velocity bulge into the corner, which is especially remarkable at x/L = 1.5.As will be indicated in § 3.2, this bulge occurs due to the momentum transport by the secondary motion and delays the corner-flow separation.Furthermore, the cross-sectional velocity contours show the secondary motion, i.e. a cross-plane flow occurs along the diagonal line toward the corner, and an outward flow occurs along the wall.At x/L = 1.0, the secondary motion is similar to the DNS square duct case Pirozzoli et al. (2018).The maximum cross-sectional velocity in the x/L = 1.0 plane is 0.017u ∞ , which agrees with Pirozzoli et al. (2018), who reported that the maximum cross-sectional velocity is approximately 2 % of the duct centreline velocity.Also, figure 7 shows the distribution of v at x/L = 1.0 with axes in the wall viscous units.Note that these wall viscous units are calculated using u τ at the centreline (y = 0).Pirozzoli et al. (2018) reported that the maximum cross-sectional velocity occurs at y + ≈ 10 and 50 z + 100, where they used the span-averaged u τ for calculating the wall viscous units.Although there are slight differences in the definition of u τ , the velocity distribution shown in figure 7 shows good agreement with their results.In the downstream (x/L = 1.5), the cross-sectional velocity is higher and more widely distributed than for x/L = 1.0.At x/L = 2.0, the secondary motion remains in a broader area compared to x/L = 1.0, which shows the history of the divergence of the flow field geometry.
Figure 8 shows the distributions of the Reynolds stress components at the streamwise cross-sections.At x/L = 1.0 and 2.0, only the near-wall peak is observed.At x/L = 1.5, the Reynolds normal stress components (u u and v v ) have a positive peak in the off-wall region.Note that w w , which is not included in figure 8, has a symmetric distribution with v v with respect to the y = z line.Also, the primary shear stress (u v ) has a negative peak at the same location.Compared to these components, the secondary shear stress v w has a relatively small magnitude overall.

Momentum transport
To clarify the influence of the secondary motion on the corner-flow separation, we investigate the transport of the streamwise momentum (more precisely, the streamwise velocity in a constant-density flow).The transport equation of the streamwise velocity is written as Here, C, P, R and V denote the convective, pressure gradient, Reynolds stress and viscous diffusion terms defined as respectively.Furthermore, to visualize the momentum transport within the plane, C, R and V within the y-z plane are rewritten as where the in-plane fluxes are defined as (3.4a-f ) Figure 9 shows each term of (3.1) and the vectors of the in-plane fluxes defined by (3.4a-f ) at x/L = 1.0, 1.5 and 2.0.In this analysis, we have confirmed that the residual (C + P + R + V) is essentially negligible in all the planes examined here; the magnitude of the residual is smaller than the minimum contour level (0.2) across the planes.Furthermore, the pressure gradient term (P) is smaller than the other components in all the planes.At x/L = 1.0 and 2.0, C and R balance away from the walls, while R and V balance in the near-wall region.The trends of each term are in qualitative agreement with the DNS square duct case by Pirozzoli et al. (2018).Moreover, the flux vectors (C y , C z ) show the convective transport of the momentum towards the corner, i.e. the effects of the secondary motion.Also, the flux vectors (R y , R z ) show that the Reynolds shear stress transports the momentum towards the wall, which is commonly known as the role of turbulence in the boundary layer.At x/L = 1.5, the magnitude of V in the near-wall region decreases compared to that at x/L = 1.0, and only C and R remain.Here, C and R along the diagonal line become more prominent compared to x/L = 1.0.Essentially, the enhanced convection induces the bulge of the streamwise velocity contours shown in figures 6(a-c).Also, as shown in the flux vectors in figures 9(g-i), the Reynolds shear stress transports the momentum from the diagonal line to the regions around (( y − y w )/L, (z − z w )/L) ≈ (0.01, 0.04) and (0.04, 0.01).This transport corresponds to the negative peak of the Reynolds shear stress shown   , (g-i) and ( j-l).For visibility, the vector length in (a-c) is halved (i.e. the magnitude corresponding to the unit vector length is twice as large as that in (g-i) and ( j-l)).Plots are for (a-c) C with vectors in figures 8(g-i).In these regions, the secondary motion lifts the flow from the wall, i.e. promotes flow separation.Conversely, the increased Reynolds shear stress suppresses flow separation in these regions.
To clarify the cause of the increase of the Reynolds shear stress at x/L = 1.5, we investigate the production of Reynolds shear stress.Figure 10 shows the production of u v near the corner at x/L = 1.3, 1.4 and 1.5.Note that only the component v v (∂ ū/∂y) of the production is shown here since the other components are almost zero.Since u v is overall negative, the production also has negative values.As shown in figures 10(a-c), the production peak position coincides approximately with the peak location of the Reynolds shear stress.Here, we further decompose the production term into v v and ∂ ū/∂y, as shown in figures 10(d-f ) and 10(g-i).These figures show that the cause of the large production is the enhanced velocity gradient at the corner.As shown in figures 10( j-l), the large velocity gradient occurs where the secondary motion distorts the streamwise velocity contours.Therefore, these results suggest that the secondary motion enhances the production of Reynolds shear stress by increasing the velocity gradient in the corner region.
In summary, the secondary motion has twofold effects of suppressing the corner-flow separation.First, the secondary motion convects momentum directly towards the corner, as shown in figures 9(a-c).Second, as an indirect effect, the secondary motion enhances turbulence production by increasing the shear.This enhanced Reynolds stress also transports momentum toward the wall, as shown in figures 9(g-i).

Vorticity transport
To investigate the generation mechanism of the secondary motion, we investigate the budget of the vorticity transport.The transport equation of the streamwise vorticity is written as where with ω y ≡ ∂u/∂z − ∂w/∂x and ω z ≡ ∂v/∂x − ∂u/∂y.Each term of (3.5) stands for the effects of convection, shear-driven vortex production, diffusion by the secondary Reynolds shear stress, transport by the turbulence anisotropy, and viscous diffusion, respectively.Furthermore, similar to the momentum budget, the terms other than S ω may be rewritten using in-plane fluxes as where the in-plane fluxes are defined as Figure 11 shows the terms on the right-hand side of (3.5).Here, S ω and the residual are essentially small at all sections and are therefore excluded from the figure.Note that figure 11 focuses on the corner region because the terms do not have noticeable distributions in the outer region.The term distributions at x/L = 1.0 are almost identical to the previous study on the square duct flow (Pirozzoli et al. 2018).Furthermore, the in-plane flux vectors, which have not been presented in the previous study, visibly show the transport by the turbulence anisotropy (A ω ) from the upper left to the lower right in the region slightly off the wall.This transport seems to be the dominant cause of the negative and positive vortex pair shown in figure 5(c).Even at x/L = 1.5 and 2.0, the anisotropy term is dominant compared to the other terms.These results indicate that the anisotropy of Reynolds normal stress (i.e. the difference between v v and w w in the current coordinate system) plays a crucial role in generating the secondary motion.

Constitutive relation for Reynolds stress
In this section, we investigate the constitutive relation between the velocity gradient and Reynolds stress tensors to develop a RANS-based turbulence model that reproduces the Reynolds normal stress accurately.For the modelling, the Reynolds stress tensor is divided into the deviatoric and TKE parts as where i and j (= 1, 2, 3) are dimensional indexes, and δ ij is Kronecker's delta.

Modelling the deviatoric part
Here, the velocity is assumed to be solenoidal since the flow considered in this study is at a low Mach number.By introducing a scalar kinematic eddy viscosity ν t , the deviatoric part of the Reynolds stress is written as where Ŝij is a second-rank tensor composed of the velocity gradient and other parameters.(d-f ) R ω with vectors (R ω,y , R ω,z ), (g-i) A ω with vectors (A ω,y , A ω,z ) and ( j-l) V ω with vectors (V ω,y , V ω,z ).
rewritten as where k, m and n are dimensional indexes, A (l) (l = 1, 2, 3, 4, 5) are parameters with the dimension of S −1 ij , and 3) has five parameters A (l) .Therefore, choosing proper A (l) and ν t reproduces the six independent components of the Reynolds stress tensor exactly.However, it is challenging to determine all six parameters in (4.3), which do not necessarily have similar magnitudes of sensitivity.Gatski & Speziale (1993) and Jongen & Gatski (1998) employ a simplified formulation, which may be rewritten in the form where α and β are also parameters with a dimension of S −1 ij .Jongen & Gatski (1998) showed that the expression defined by (4.4) is the best approximation when taking only two of the quadratic terms.We adopt this expression to parametrize the constitutive relation more easily than (4.3).Note that Sabnis et al. (2021) also adopted the same form of the quadratic terms.As Modesti (2020) showed, increasing the number of terms may improve the expression of the Reynolds stress components.However, increasing the number of terms also increases the complexity of the formulation and difficulty of the parameter study.Therefore, we include only two of the quadratic terms to retain the conciseness of the model.
For convenience, we redefine the above expression using the formulation like QCR2000 (Spalart 2000) as where u x is the magnitude of the velocity gradient, i.e. u x ≡ (∂ ūm /∂x n ) 2 .Also, C q1 and C q2 are the parameters to control the anisotropy, which are not necessarily constant in space.For example, the constant parameter pairs (C q1 , C q2 ) = (0.0, 0.0) and (C q1 , C q2 ) = (0.6, 0.0) give LCR and QCR2000, respectively.We seek optimal values of these parameters through a priori testing using the LES results.
To evaluate the validity of the constitutive relation, we employ the tensorial inner product introduced by Schmitt (2007).The tensorial inner product is defined as which represents the cosine of the angle between Ŝij and R ij .If σ R Ŝ = 1, then the product of Ŝij and ν t perfectly represents R ij by properly choosing a positive scalar for ν t .Note that we leave the model of ν t and focus on the validity of the constitutive relation.The evaluation of ν t depends closely on the employed baseline turbulence model (e.g.Spalart & Allmaras (1992) or k-ε models), which should be validated in a different context.First, we seek the optimum values of C q1 and C q2 using the LES result at x/L = 1.0 (i.e.upstream of the diverging section).Here, we pick three locations A, B and C, shown in figure 12(a).These locations represent the secondary vortex centre, off-wall location away from the side wall, and near-wall inner-layer location (y + ≈ 15), respectively.Figures 12(b-d) show σ R Ŝ for C q1 and C q2 varying in [0, 10], and figures 12(e-g) are close-up views of figures 12(b-d).At locations A and B, σ R Ŝ increases as C q1 increases from zero, and takes its maximum at C q1 ≈ 1.The dependency on C q2 is weak compared to that on C q1 .With C q1 = 1, 0.15 C q2 1.3 gives σ R Ŝ > 0.99.At location C, the optimum values for C q1 and C q2 are much higher than those at the other two locations, suggesting the strong turbulence anisotropy in the near-wall region.Based on the results in figure 12, we choose the parameters (C q1 , C q2 ) = (1.0,0.5).To investigate the generality of these values, we calculate σ R Ŝ over the planes x/L = 1.0, 1.5 and 2.0. Figure 13 shows the distributions of σ R Ŝ over these three planes.Here, we compare the candidate value pair (C q1 , C q2 ) = (1.0,0.5) to QCR2000 (C q1 , C q2 ) = (0.6, 0.0) and LCR (C q1 , C q2 ) = (0.0, 0.0).As shown in this figure, the parameter pair (C q1 , C q2 ) = (1.0,0.5) gives overall higher values of σ R Ŝ compared to QCR2000 and LCR.
Furthermore, we check the quantitative validity of the parameters.For this purpose, we need to determine ν t .ν t may be determined uniquely by assuming a two-dimensional (2-D) simple shear flow, i.e. a flow where the velocity gradient tensor consists of only the ∂ ū/∂y component.This assumption is almost valid in the region sufficiently away from the side wall.When the velocity gradient tensor consists of only the ∂ ū/∂y component, the components of the traceless stress tensor (4.2) become   parameters are constant in space, the results here suggest that the proposed parameter pair (C q1 , C q2 ) = (1.0,0.5) has generality to a certain extent in wall turbulence, except for the near-wall inner layer.Note that the formulation of the proposed QCR leaves room for setting (C q1 , C q2 ) as spatial variables, although only the constant parameter pair is employed in the following.For example, the inner-layer peaks of the Reynolds stress components can be reproduced by increasing (C q1 , C q2 ) locally.Appendix B describes the attempt to develop a correction to reproduce the inner-layer peaks.
The analysis above suggests that C q1 should be larger than the value adopted by QCR2000 (i.e.0.6).However, Spalart, Garbaruk & Strelets (2014) and Leger, Bisek & Poggie (2016) pointed out that increasing C cr1 of QCR2000, which is equivalent to C q1 /2 in this study, may cause non-physical vortices.The proposed model avoids this by introducing C q2 as explained in the following.Equation (3.5) shows that the generation of the streamwise vortex depends on the difference between the Reynolds normal stress components v v − w w ≡ −R yy + R zz .From (4.7), we obtain for which Sabnis et al. (2021) derived an equivalent equation.Equation (4.8) suggests that C q2 partially cancels C q1 .Therefore, introducing C q1 = 1.0 alone may lead to too strong secondary vortices, and the balance between C q1 and C q2 is important for predicting the secondary vortices accurately.

4.2.
Modelling the TKE part Some of the existing RANS-based turbulence models, such as the SA turbulence model (Spalart & Allmaras 1992), do not contain the TKE part.Previous studies (Mani et al. 2013;Rumsey et al. 2020) introduced an estimation of TKE using the structure parameter presented by Bradshaw (1967).The estimation used by Rumsey et al. (2020) is written as where C k = 1/(3a 1 ), and a 1 = 0.155 is Bradshaw's structure parameter (i.e.C k = 2.15).
Figure 15 shows the wall-normal profiles of the estimated TKE at the same locations as in figure 14.The estimation by (4.9) is in good agreement with the actual TKE in the LES at all the locations examined here, except for the region near the wall.Furthermore, the correction for the inner layer is presented in Appendix B, similar to the deviatoric part.

Implementation of RANS simulations
In this section, the constitutive relations presented in § 4 are validated through RANS simulations of side-wall interference flows.The proposed constitutive relation is summarized as which contains three parameter coefficients, C q1 = 1.0,C q2 = 0.5 and C k = 2.15.In a compressible flow solver, S ij is replaced by S ij − S kk δ ij /3, similar to the existing QCRs (e.g.Spalart 2000;Rumsey et al. 2020).Here, we calculate ν t in (5.1) using the SA turbulence model, although the proposed QCR may be combined with any turbulence model based on the Boussinesq approximation.For a comparison, simulations are conducted using QCR2000 (Spalart 2000), QCR2020 (Rumsey et al. 2020) and LCR.In the following, we denote the proposed QCR defined by (5.1) as QCR(r) and the QCRs combined with the SA turbulence model as SA-QCR(r), and so on.

Diverging-converging side-wall interference flow
5.1.1.Computational settings First, the same flow field as in the LES is simulated.Figure 16 shows the overview of the computational grid for RANS simulations.Here, the grid spacing in the streamwise direction is 0.005L.The minimum grid spacing in the wall-normal direction is the same as in the LES, while the stretching ratio is approximately 8 %.The resulting grid dimensions are 461 × 91 × 91.We have confirmed the grid convergence of the result as summarized in Appendix C. The f t2 function of the SA turbulence model (Spalart & Allmaras 1992) is set to 1.0 in x/L < 0.1 and 0.0 in x/L 0.1 to reproduce the forced laminar-turbulent transition.
In the simulation, the convective term is evaluated by the SLAU scheme (Shima & Kitamura 2011) with the third-order monotone upwind-central scheme for conservation law (van Leer 1979).The van Albada limiter (van Albada, van Leer & Roberts 1982) is applied only to the SA model equation.In addition, the viscous term and the diffusion and source terms of the SA turbulence model are evaluated by the second-order central difference.The time is integrated using the implicit method (Obayashi et al. 1988;Izuka 2006) with a local time stepping, where the maximum Courant number is approximately 7. The spectral radius for the viscous term on the left-hand side is multiplied by (1 + C q1 + C q2 + C k ) to enhance numerical stability.

Results
Figure 17 compares the obtained streamwise velocity distributions over the streamwise cross-sections.The results using SA-QCR(r) predict the size of the separation bubble similar to the LES result (figure 5a).Here, SA-QCR2000 and SA-QCR2020 overpredict the corner-flow separation, although those results show slight improvement from the results using SA-LCR.Figures 18 and 19 compare the streamwise and cross-sectional velocity distributions to the LES data.The prominent bulge of the streamwise velocity contours at x/L = 1.5 is reproduced only by SA-QCR(r).Here, the distribution of the cross-sectional velocity (figures 19a-c) indicates two counter-rotating streamwise vortices.Although the magnitude of the cross-sectional velocity is slightly overestimated compared to the LES data, SA-QCR(r) well predicts the qualitative features of the in-plane flow motions.In contrast to SA-QCR(r), SA-QCR2000 and SA-QCR2020 show different trends of the velocity contours at x/L = 1.5 and 2.0.The in-plane velocity even has a different sign for these two results, suggesting that the overestimated flow separation causes a strong in-plane flow that conceals the secondary motion.Note that the results using SA-QCR2000 and SA-QCR2020 are non-symmetric about the corner.In these cases, the shape of the separation fluctuates during the iterations, and the flow field does not fully converge (see the residual history shown in figure 20).In contrast to these cases, the computation using SA-QCR(r) converges to a steady state.
The close-up views of the cross-sectional velocity contours at x/L = 1.0 are shown in figure 21 to validate the basic performance of the proposed QCRs for reproducing the secondary motion.The three RANS results consistently show that the secondary motion is distributed in a slightly wider area compared to the LES, which may be due to the difference in the boundary layer thickness at this location.When using SA-QCR(r), the negative peak value of the velocity is close to the LES.On the other hand, SA-QCR2000 and SA-QCR2020 slightly underpredict the magnitude of the in-plane velocity.
Figure 22 shows C p and C f distributions at the periodic boundary plane y/L = 0. Due to the differences in the bubble size, the C p and C f distributions in 1.0 < x/L < 2.0 differ between the cases.Here, the results obtained by SA-QCR(r) show good agreement with the LES data except for C f in the vicinity of the transition location (x/L = 0.1) and the recovery region (1.5 < x/L < 2.0).These slight discrepancies are assumed to be due to the baseline SA turbulence model because all the cases consistently show the same trend.The results using SA-QCR2000 and SA-LCR deviate from the LES data in 1.0 < x/L < 1.5, which suggests the influence of the different separation bubble size in these two cases.Furthermore, figure 23 compares the spanwise variation of C f at x/L = 1.0, 1.5 and 2.0.SA-QCR(r) shows overall reasonable agreement with the LES result at all the cross-sections in the figure.Also, although not shown here, the spanwise distributions of C p obtained by SA-QCR(r) show good agreement with the LES result.A minor difference between SA-QCR(r) and the LES remains in the dimple of C f slightly x/L = 1.4 x/L = 1.5 x/L = 1.6 x/L = 1.4 x/L = 1.5 x/L = 1.6 x/L = 1.4 x/L = 1.5 x/L = 1.6 x/L = 1.4 x/L = 1.5 x/L = 1.6 0.4 0.6 0.8 1.0 away from the side wall (( y − y w ) ≈ 0.04 in the x/L = 1.5 plane).This dimple is due to the flow departing from the wall shown in figures 19(a-c), which is slightly overestimated by SA-QCR(r).Moreover, we examine the reproducibility of the Reynolds stress components by the proposed QCRs. Figure 24 shows the distributions of the Reynolds stress components computed using SA-QCR(r).As shown here, the magnitude and qualitative features of the Reynolds stress components are in reasonable agreement with the LES data shown in figure 8.Note that v v at x/L = 1.5 (see figure 24e) appears to be different from the LES data.This discrepancy is noticeable because v v along the diagonal line y = z is underestimated.Except for this difference, the magnitude and the overall distributions of the Reynolds stress components are well predicted by SA-QCR(r).

Computational settings
To investigate the robustness of the proposed QCRs to different flow conditions, we simulate the supersonic square duct flow.The problem settings follow the definition presented in the NASA turbulence modelling resource (TMR; see https://turbmodels.larc.nasa.gov/).The Reynolds number based on the duct height D and the inlet velocity u ∞ is 5.08 × 10 5 , and the inlet Mach number is 3.9.In the calculation, the convection flux is computed by flux-difference splitting (Roe 1981), and the other computational methods are the same as in the problem in § 5.1.The computational grid is taken from the TMR and has dimensions 961 × 161 × 161 in each direction.As a verification, we have compared the computational result using SA-QCR2000 to that provided in the NASA TMR using FUN3D with the same turbulence model, and confirmed that the two results are almost identical.

Results
Figure 25 shows the obtained velocity and C f near the outlet (x/L = 50), where the flow is considered to be fully developed.Here, the velocity is normalized by the duct centreline  1.5 and (c, f,i,l) x/L = 2.0, with in-plane velocity vectors obtained by the RANS simulations: (a-c) SA-QCR(r) (proposed), (d-f ) SA-QCR2000, (g-i) SA-QCR2020, ( j-l) LES (reference).velocity u CL .The cross-sectional velocity v obtained by SA-QCR(r) has a slightly larger value than that from SA-QCR2000, and is close to SA-QCR2020. Figure 25     Furthermore, the streamwise velocity distribution along the diagonal (y = z) and centreline (2y/D = 1.0) are shown in figure 27.Along the diagonal line, the results obtained by SA-QCR(r) and SA-QCR2020 show good agreement with the experimental data, while that obtained by SA-QCR2000 deviates slightly from the other results.This corresponds to near the corner, suggesting that SA-QCR2000 underestimates the momentum transport towards the corner.Along the centreline, all the results show a slower streamwise velocity in the near-wall region, which may be due to the SA turbulence model.Despite the minor difference near the centreline, the obtained results confirm that SA-QCR(r) reproduces robustly the secondary motion at the different flow conditions.

Conclusions
We conducted a wall-resolved LES of a side-wall interference flow involving corner-flow separation.The objectives of this LES are first to clarify the influence of the secondary motion on the corner-flow separation, and second to explore turbulence modelling suitable for the side-wall interaction flow, such as the wing-body junction of aircraft.
For the first objective, the streamwise momentum budget analysis indicated that the secondary motion (i.e. the streamwise vortex) has direct and indirect effects on the corner-flow separation.The direct effect is the convection of the high streamwise momentum into the corner region.This convection causes the bulge of the streamwise velocity distribution and reduces the size of corner-flow separation.On the other hand, the indirect effect is enhanced turbulence production due to the increased shear.This enhanced turbulence transports the momentum towards the wall, which essentially reduces the corner-flow separation.
As indicated in the vorticity transport analysis, the secondary motion is generated primarily by the turbulence anisotropy.Therefore, regarding the second objective, we sought turbulence modelling based on the constitutive relation between the Reynolds stress and velocity gradient.We introduced a QCR (called QCR(r) in this study) with three parameters C q1 , C q2 and C k , which were adjusted by the LES data.The analysis revealed that even the constant-parameter QCR represents accurately the turbulence anisotropy over most of the computational domain.The notable feature of QCR(r) is that the formulation does not contain variables dependent on the baseline turbulence model.Therefore, the proposed QCR may be combined with any baseline turbulence model based on the Boussinesq approximation.
In the RANS testing, the proposed QCR combined with the SA turbulence model, SA-QCR(r), predicts the size of the corner-flow separation better than the existing models (i.e.SA-QCR2000 and SA-QCR2020).Also, the Reynolds stress distributions predicted by SA-QCR(r) show qualitatively good agreement with the LES data.The trend of the results is also consistent in the three-dimensional supersonic duct problem, which confirms the robustness of the proposed modelling for different Mach and Reynolds numbers.Although the universality of the model should be investigated further in the future by comparing it to other turbulence models (e.g. the Reynolds stress models), the proposed QCR may provide a simple way to reproduce the turbulence anisotropy and resulting secondary motions.difference between the averaged flow field for the period F + L and that for the period L. Note that the spatial averaging procedure investigated later is applied to both results.Figure 28 compares the Reynolds shear stress u v distributions at x/L = 1.5 obtained by averaging for these two different time periods.The results show that the difference between these two results is considerably smaller than the magnitude of the Reynolds shear stress.Also, we have confirmed that the other components of the Reynolds stress tensor or mean velocity essentially do not change by changing the averaging period.Therefore, we conclude that the averaging period F + L employed in this study is adequate, at least for the discussions in this study.
Moreover, we investigate the effects of spatial averaging in the streamwise direction described in § 2.4.Since high-order differential quantities tend to contain noise, we compare the vorticity budget described in § 3.3.Figure 29 compares the terms of the vorticity budget at x/L = 1.0 with and without the averaging.As shown here, the spatial averaging reduces the high-frequency noise but does not change the overall distributions of each term.

Appendix B. Inner-layer correction for the near-wall turbulence
Here, we demonstrate that the near-wall distribution of Reynolds normal stress may be reproduced more accurately by treating C q1 , C q2 and C k as spatial variables.Note that the following modelling is limited to attached fully developed boundary layers, and more intensive investigation will be needed for separated boundary layers.To model the inner layer, we assume a 2-D simple shear flow.From (4.7) and (4.9), the Reynolds stress components in a 2-D simple shear flow are written as It is known that u u , v v , w w and u v in the near-wall region are proportional to y 2 , y 4 , y 2 and y 3 , respectively (see Pope 2000, pp. 283-285).Therefore, C q1 + 1 6 C q2 + C k and − 1 3 C q2 + C k must be proportional to y −1 , while −C q1 + 1 6 C q2 + C k is proportional to y.To reproduce these wall-asymptotic behaviours, we assume the following function forms: The functions f b1 and f b2 are assumed to be functions of y + since the inner-layer behaviours of the Reynolds stress components are relatively robust for wide ranges of Mach and Reynolds numbers.For the inner part, we consider that f b1 is proportional to y −1 , such as where ε is a small constant to avoid division by zero (ε = 10 −12 in this study).Here, y + inner is a certain height for calibrating the near-wall behaviour, for which we choose y + inner = 10 as a representative location within the inner layer.Also, to ensure the asymptotic behaviour v v ∼ y −4 , −C q1,inner + 1 6 C q2,inner + C k,inner must be zero to eliminate the leading term with y −2 .By considering this condition and substituting the LES data at y + = 10 (at (x/L, z/L) = (1.0,−0.12)) into (B1), the inner layer parameters are calculated as (C q1,inner , C q2,inner , C k,inner ) = (8.5, 12.0, 6.5).Also, we found that the original blending function f b1 defined by (B4) diminishes too slowly in the outer layer.Therefore, we modify f b1 by introducing an additional damping function as where a + = 60 is an empirically chosen parameter.The additional damping function in (B5) works to decrease f b1 in y + a + , and does not change the near-wall behaviours.Furthermore, the parameters for the outer part, C q1,outer , C q2,outer , C k,outer , are equal to the values determined in § § 4.1 and 4.2, namely, 1.0, 0.5, 2.15.The function f b2 is determined considering v v ∝ y 2 , which requires f b2 ∝ y.Therefore, we introduce the function where b + = 25.
Figure 30 shows the Reynolds normal stress components calculated by the above method at (x/L, z/L) = (1.0,0.0).Here, the proposed method reproduces accurately the stress components, including their inner-layer peaks.The results also show that the stress components behave as intended, i.e. u u ∝ y 2 , v v ∝ y 4 and w w ∝ y 2 .
For implementation into RANS-based turbulence models, the use of y + ≡ u τ y/ν is inconvenient because it requires u τ at the nearest wall.It is difficult to search the nearest wall, especially in parallel simulation with domain splitting.To avoid the use of non-local quantities, we estimate u τ locally as where |S| ≡ 2S ij S ij .This estimation is based on the near-wall equilibrium assumption, i.e. the sum of the viscous and Reynolds shear stresses is equal to the wall shear stress τ w in the near-wall region.Figure 31 compares y + using (B7) to the actual y + .The estimation is very accurate in the inner layer (y + 30) when the pressure gradient is weak (i.e.x/L = 1.0 and x/L = 2.0), and still reasonable at x/L = 1.5.Note that |S| becomes zero in the freestream, resulting in y + → 0. Although this causes an unintended increase of the parameters C q1 and C q2 in the freestream with |S| = 0, the increase in the parameters essentially does not influence flows because |S| = 0 means that the viscous stress is zero.The estimation may also become inaccurate in a boundary layer with a strong pressure gradient, which is outside the scope of the inner-layer correction presented here.Moreover, to retain generality in three-dimensional fields, the distance from the nearest wall d is used instead of y.This d is readily available in most flow solvers because some turbulence models (e.g. the SA turbulence model) require it.It may be calculated by the level-set method (Sussman et al. 1999), which does not require direct searching of the nearest wall.We also note that the SA turbulence model does not satisfy the asymptotic behaviour of u v ∝ y 3 .Therefore, to reproduce the intended wall asymptotic behaviour, we have to replace the wall damping function f v1 of the SA model with a function satisfying f v1 ∝ y 2 , such as where A + = 17.We have implemented these inner-layer corrections to RANS simulations, and confirmed that their effects on the prediction of the corner-flow separation are minor.This result suggests that the inner-layer flow physics may not be crucial in forming the separation.Although not necessary for all types of flows, the inner-layer corrections shown here indicate the potential of the proposed QCR formulation for reproducing the strong turbulence anisotropy in the inner layer.

Figure 1 .Figure 2 .
Figure 1.Geometry of the side wall interference flow field.

Figure 4 .
Figure 4. Mean streamwise velocity, pressure coefficient and skin friction coefficient distributions over the periodic boundary plane (y/L = 0) obtained by the LES.

Figure 5 .
Figure 5. Overview of the statistically averaged flow field in the diverging and converging sections obtained by the LES: (a) mean streamwise velocity ū/u ∞ ; (b) TKE K/u 2 ∞ ; (c) mean streamwise vorticity ωx L/u ∞ .

Figure 6 .Figure 7 .
Figure 6.Mean velocity distributions over the streamwise cross-sections near the corner at (a,d) x/L = 1.0, (b,e) x/L = 1.5 and (c, f ) x/L = 2.0, where z w is the coordinate of the side wall.In (d-f ), negative contours are shown with black dashed lines (also applied to the following figures).Note that the scale of the in-plane velocity vectors varies with cross-section location.(a-c) Streamwise velocity ū/u ∞ .(d-f ) Cross-sectional velocity v/u ∞ with in-plane velocity vectors.

Figure 9 .
Figure 9. Streamwise momentum budget near the corner at (a,d,g, j) x/L = 1.0, (b,e,h,k) x/L = 1.5 and (c, f,i,l) x/L = 2.0.Each term is normalized by u 2 ∞ /L.The in-plane fluxes (3.4a-f ) are overlaid as vectors in (a-c), (g-i) and ( j-l).For visibility, the vector length in (a-c) is halved (i.e. the magnitude corresponding to the unit vector length is twice as large as that in (g-i) and ( j-l)).Plots are for (a-c) C with vectors (C y , C z ), (d-f ) P, (g-i) R with vectors (R y , R z ) and ( j-l) V with vectors (V y , V z ).
For example, Ŝij = S ij ≡ 1/2(∂ ūi /∂x j + ∂ ūj /∂x i ) for the standard eddy viscosity model based on the Boussinesq approximation (i.e. the linear constitutive relation, LCR).Also,Lumley (1970) introduced a general expression of the constitutive relation, which may be y w )/L ( yy w )/L ( yy w )/L (z -z w

Figure 12 .
Figure 12.Tensorial inner product σ R Ŝ (4.6) at probe locations in the x/L = 1.0 plane as a variable of parameters C q1 and C q2 .In (b-d) and (e-g), results at probe locations A, B and C are shown from left to right.(a) Probe locations (in-plane velocity vectors overlaid) with (b-d) σ R Ŝ, (e-g) σ R Ŝ (close-up).

Figure 15 .
Figure15.Wall-normal profiles of the estimated kinetic energy at (a) x/L = 1.0,(b) x/L = 1.5 and (c) x/L = 2.0, in the z/L = 0.0 plane.Symbols and lines denote the LES data estimation using (4.9), respectively.

Figure 16 .
Figure 16.The RANS computational grid at spanwise and streamwise cross-sections: (a) x-z plane at y/L = 0.0 (every 10 grid points are shown); (b) y-z planes (every 5 grid points are shown).
shows the C f distribution along the wall at x/D = 50.Near the corner (2z/D 0.2), the results obtained by SA-QCR(r) and SA-QCR2020 are almost identical.At the centreline, the result obtained by SA-QCR(r) shows a slight decrease of C f .Similar to the dimple of C f observed in figure 23, the slight decrease of C f may be due to the upward flow from the wall along 2z/D = 1.0 in figure 26.

Figure 25 .Figure 26 .
Figure 25.Plots of C f along the wall at x/D = 50.Symbols denote the reference experimental data (Davis & Gessner 1989).Lines are as in figure 22.

Figure 28 .
Figure 28.Reynolds shear stress u v /(u 2 ∞ ) distributions near the corner at x/L = 1.5 with different averaging periods: (a) averaged for period F + L; (b) averaged for period L; (c) difference between the two results.