NOMENCLATURE
- Cd
-
drag coefficient
- Cf
-
skin friction coefficient
- Cp
-
pressure coefficient
- h
-
enthalpy
- H
-
flight height, km
- k
-
turbulent kinetic energy, m2/s2
- kr
-
equivalent sand grain height, mm
- k + r
-
dimensionless roughness height, k + r = kru τ/v
- L
-
flat-plate length, m
- M
-
mach number
- Re
-
reynolds number, 1/m
- S
-
strain rate magnitude, 1/s
- T
-
temperature, K
- uτ
-
friction velocity, uτ = (τw/ρ)0. Reference Busse and Sandham 5
- y
-
wall distance, mm
- ∆Cd %
-
relative rough-wall Cd derivation from the smooth-wall case
- ∆u
-
vertical velocity shift of the logarithmic profile
- γg
-
specific heat ratio, γg = 1.4
- μ
-
dynamic viscosity, μ = ρν
- ν
-
kinetic viscosity, m2/s
- νT
-
eddy viscosity, m2/s
- τw
-
wall shear stress, τw = ρν ∂u/∂y|w
- ω
-
specific turbulence dissipation rate
- o
-
reservoir conditions
- w
-
wall
- ∞
-
free-stream value
Subscripts
1.0 INTRODUCTION
In high-speed flight, the surfaces on rocket nozzles or re-entry nosecaps are always rough due to ablation. As the roughness height increases, the onset of laminar-turbulent flow transition gradually moves upstream(Reference Gianluca, Iaccarino and Shaqfeh1,Reference Fu and Wang2) . Once the boundary layer becomes turbulent, heat transfer can increase by a factor of 10. Therefore, the accurate and reliable prediction of the effects of surface roughness on fluid flow and heat transfer is critical for the design and optimisation of future space vehicles operating at sustained supersonic and hypersonic speeds. However, compared to the researches of the surface roughness effect on the aerofoils with icing or turbine blades, very few experimental investigations of the corresponding high-speed cases have been completed up to date(Reference Bowersox3).
There are various approaches to predict flow over rough surfaces. Hierarchically, these approaches involve, in order of descending requirement in computer resource, DNS (direct numerical simulation), LES (large eddy simulation), DEA (discrete element approach) and equivalent sand grain approach. Lu and Liou(Reference Lu and Liou4) and Busse and Sandham(Reference Busse and Sandham5) reviewed the progress and state of DNS and LES of flow around roughness elements, which is still limited to simple geometries at a research level. The idea of DEA is to modify the mean flow equations to account for blockage effects due to the presence of roughness elements as well as drag and heat flux on the roughness elements. This technique was introduced by Robertson Reference Robertson6 and further refined by Finson(Reference Finson7), Coleman et al(Reference Coleman, Hodge and Taylor8) and McClain et al(Reference McClain, Hodge and Bons9,Reference McClain, Collins, Hodge and Bons10) . However, it is difficult for DEA to be implemented into a general-purpose CFD (Computational Fluid Dynamics) code for they require altering the flow equations. There is no prospect of adequately handling non-uniform roughness effects in the Reynolds-Averaged Navier-Stokes (RANS) environment(Reference Aupoix11). The equivalent sand grain approach is thus adopted in the present study because it is the most popular and the only engineering approach.
The concept of equivalent sand grain roughness was introduced by Schlichting(Reference Schlichting12), to relate any rough surface to the sand grain size that would yield the same drag increase in the fully rough regime. The sand grain height, kr +, is introduced as a new parameter in the RANS turbulence model, either with wall functions(Reference Suga, Craft and Iacovides13,Reference Apsley14) or with near-wall resolution(Reference Knopp, Eisfeld and Calvo15-Reference Durbin, Medic, Seo, Eaton and Song22), to enhance turbulence in the wall region. However, the drag increase is essentially due to pressure forces on the roughness elements, while this approach only increases the frictional drag. Although it is unphysical and its performance depends on the turbulence model selected, the present approach has been widely accepted due to the simplicity to implement in RANS solvers as well as its comparably high efficiency in calculations(Reference Aupoix11).
In turbulence models with the equivalent sand grain approach, the rough surface is replaced by an effective, smooth surface, on which modified boundary conditions are imposed. For low-Re turbulence models, an effective displacement of the wall distance origin was introduced by Aupoix and Spalart for the SA (Spalart and Allmaras) one equation model(Reference Aupoix and Spalart18) and by Durbin et al for the two layer k–ε model(Reference Durbin, Medic, Seo, Eaton and Song22), respectively. The effective displacement is related to a hydrodynamic roughness length that is used to modify turbulence length scales as well as the boundary condition for turbulent quantities. The similar approach was also proposed to extend the k–ω type models to account for wall roughness. An early example is the Wilcox roughness modification(Reference Wilcox19). However, the Wilcox model requires a very fine mesh resolution and is not accurate for transitionally rough walls. The subsequent models by Seo(Reference Seo21) and Knopp et al(Reference Knopp, Eisfeld and Calvo15) give satisfactory results with near wall grid spacing similar to that for smooth walls. It is noted that in the aerodynamically smooth regime, the roughness elements are embedded in the viscous sublayer; whereas in the fully rough regime, the drag increase is only due to pressure forces on the roughness elements; the transitionally rough regime is between these two regimes, where both viscous and pressure forces on the roughness elements contribute to the drag. So far, the applications of these roughness modifications have been limited to low-speed flow regime(Reference Aupoix11).
Furthermore, within the RANS framework, a few recent studies propose roughness-induced transition models. They are based on an empirical correlation for the momentum thickness Reynolds number at which transition starts. The correlation in Stripf et al(Reference Stripf, Schulz, Bauer and Wittig23) depends on both kr + and density, while Boyle and Stripf(Reference Boyle and Stripf24) propose a simpler formula that only depends on k + r . Elsner and Warzecha(Reference Elsner and Warzecha25) introduced the roughness transition correlation by Stripf et al(Reference Stripf, Schulz, Bauer and Wittig23) into the γ-Reθt transition model proposed by Langtry and Menter(Reference Langtry and Menter26) for smooth walls, while Dassler et al(Reference Dassler, Kozŭlović and Fiala27) extended the γ-Reθt model to account for wall roughness in a different type. Both of them chose Wilcox's roughness boundary condition(Reference Wilcox19) for the fully turbulent boundary layer. Similarly, Ge and Durbin(Reference Ge and Durbin28) considered the wall roughness effect based on their own smooth-wall transition model(Reference Ge, Arolla and Durbin29).
Roy and Blottner reviewed and assessed turbulence models for hypersonic flows(Reference Roy and Blottner30). It is found that the Menter k–ω SST (Shear Stress Transport) model(Reference Menter31) with compressibility corrections(Reference Sarkar, Erlebacher and Hussaini32,Reference Sarkar33) preforms the best among a total of 18 one- and two-equation turbulence models. The compressibility corrections have been carefully examined for their effects on a hypersonic validation database. With the use of SST model with the compressibility corrections, the present authors have conducted high-speed simulations with a reasonably wide range of Mach number(Reference Fu and Wang2,Reference Wang and Fu34-Reference Wang, Xiao and Fu38) , such as supersonic and hypersonic flows past cones at small incidences and elliptic cones at zero incidences. The objective of this work is thus to simulate supersonic/hypersonic flows over a flat plate with transitionally and fully rough surface, by using the SST model with advanced compressibility as well as roughness corrections. Particularly, a modified roughness correction is proposed and validated by the present authors.
This article is organised as follows: Section 2 presents the mathematical formulation of the SST turbulence model including the modifications needed to represent the effects of sand-grain roughness as well as flow compressibility; Section 3 gives the numerical details, especially for boundary conditions; the comparison of the predictions with reference values from the literature for model validation is presented in Section 4; while Section 5 includes the numerical results on systematically refined grids to investigate the drag increase due to wall roughness for hypersonic flat-plate flows; the conclusions of this study are summarised in Section 6.
2.0 MATHEMATICAL FORMULATION
2.1 Menter's k–ω SST model
This model requires the solution of transport equations for k and ω:
$$\begin{equation}
\frac{{\partial \rho k}}{{\partial t}} + \frac{\partial }{{\partial {x_j}}}\left[ {\rho {u_j}k - \left( {\mu + {\sigma _k}{\mu _t}} \right)\frac{{\partial k}}{{\partial {x_j}}}} \right] = {\tau _{ij}}{S_{ij}} - {\beta ^*}\rho k\omega
\end{equation}$$
and
$$\begin{equation}
\frac{{\partial \rho \omega }}{{\partial t}} + \frac{\partial }{{\partial {x_j}}}\left[ {\rho {u_j}\omega - \left( {\mu + {\sigma _\omega }{\mu _t}} \right)\frac{{\partial \omega }}{{\partial {x_j}}}} \right] = {P_\omega } - \beta k{\omega ^2} + 2\left( {1 - {F_1}} \right)\frac{{\rho {\sigma _{\omega 2}}}}{\omega }\frac{{\partial k}}{{\partial {x_j}}}\frac{{\partial \omega }}{{\partial {x_j}}},
\end{equation}$$
where S is the mean strain rate, and F 1 is a blending function expressed as
$$\begin{equation}
{F_1} = {\rm{tanh}}\left\{ {{{\left\{ {{\rm{min}}\left[ {{\rm{max}}\left( {\frac{{\sqrt k }}{{0.09\omega d}},\frac{{500\mu }}{{\rho {d^2}\omega }}} \right),\frac{{4\rho {\sigma _{\omega 2}}k}}{{C{D_{k\omega }}{d^2}}}} \right]} \right\}}^4}} \right\},
\end{equation}$$
with
$$\begin{equation}
C{D_{k\omega }} = {\rm{max}}\left( {\frac{{2\rho {\sigma _{\omega 2}}}}{\omega }\frac{{\partial k}}{{\partial {x_j}}}\frac{{\partial \omega }}{{\partial {x_j}}},{{10}^{ - 20}}} \right),
\end{equation}$$
and the model constants(Reference Menter31) depend on F 1. The eddy-viscosity μT is calculated from
$$\begin{equation}
{\mu _t} = \min \left( {\frac{{\rho k}}{\omega },\frac{{\rho {a_1}k}}{{\Omega {F_2}}}} \right),
\end{equation}$$
where a 1=0.31, and F 2 is also a blending function:
$$\begin{equation}
{F_2} = \tanh \left\{ {{{\left[ {\max \left( {\frac{{2\sqrt k }}{{0.09\omega d}},\frac{{500\mu }}{{\rho {d^2}\omega }}} \right)} \right]}^2}} \right\}
\end{equation}$$
2.2 Velocity shift by wall roughness
A key parameter to characterise the roughness effects is the dimensionless roughness height
$$\begin{equation}
k_r^ + = {k_r}{u_\tau }/\nu ,
\end{equation}$$
where kr is the equivalent sand grain height, uτ = (τw/ρ)0.Reference Busse and Sandham5 the friction velocity based upon the wall shear stress, τw = ρν ∂u/∂y|w , the density ρ and the viscosity ν. Nikuradse proposed k + r = 3.5 and k + r = 68 as the limits of the transitional roughness regime(Reference Nikuradse39).
Near the wall, the flow is highly perturbed by the presence of the roughness elements. Nikuradse(Reference Nikuradse39) pointed out that, above the roughness sublayer, the logarithmic law is preserved but shifted. The velocity profile can be described by
$$\begin{equation}
\frac{u}{{{u_\tau }}} = \frac{1}{\kappa }{\rm{ln}}\left( {\frac{y}{{{k_r}}}} \right) + B,
\end{equation}$$
where κ of 0.41 is Karman constant, B is of 8 under fully rough conditions (k + r > 68), while for hydrodynamically smooth walls (k + r < 3.5), the classical log-law
$$\begin{equation}
\frac{u}{{{u_\tau }}} = \frac{1}{\kappa }\ln \left( {\frac{{y{u_\tau }}}{\nu }} \right) + 5.1
\end{equation}$$
is formally recovered by setting B = κ −1 ln (k + r ) + 5.1. Then, if ∆u, the vertical shift of the logarithmic profile, is defined in terms of u by
$$\begin{equation}
\frac{u}{{{u_\tau }}} = \frac{1}{\kappa }\ln \left( {\frac{{y + \Delta u}}{{\Delta u}}} \right)
\end{equation}$$
Equations (8) and (10) give
$$\begin{equation}
\frac{y}{{{k_r}}}{e^{B\kappa }} = \frac{{y + \Delta u}}{{\Delta u}}
\end{equation}$$
Using the approximation that ∆u/y is infinitesimal of higher order, the dimensionless velocity shift can thus be expressed as
$$\begin{equation}
\Delta {u^ + } = \frac{{\Delta u \cdot {u_\tau }}}{\nu } = {e^{ - B\kappa }}\frac{{{k_r} \cdot {u_\tau }}}{\nu } = {e^{ - B\kappa }}k_r^ +
\end{equation}$$
2.3 Roughness corrections
In the present approach, the rough surface is replaced by an effective, smooth surface, on which new boundary conditions are imposed. Under fully rough conditions, the log-layer solution k = uτ 2/Cμ 0.5, with Cμ = 0.09, extends to the effective wall origin, where the log-layer eddy viscosity νT = uτκ (y + ∆u) reduces to νT = uτ κ∆u. From the definition νT = k/ω, the boundary condition for ω should be
$$\begin{equation}
\omega = \frac{{{u_\tau }}}{{C_\mu ^{0.5}\kappa \Delta u}}
\end{equation}$$
Generally, the present work adopts the ω boundary condition proposed by Knopp et al(Reference Knopp, Eisfeld and Calvo15), as represented as
$$\begin{equation}
{\omega _w} = {\rm{min}}\left( {\begin{array}{*{20}{l}} {\frac{{{u_\tau }}}{{C_\mu ^{0.5}\kappa \Delta u{\phi _{r2}}}},}&{\frac{{800\nu }}{{{y_1}}}} \end{array}} \right),
\end{equation}$$
where ∆u equals to 0.03 kr according to Formula (12), y 1 is the grid point next to the wall and
$$\begin{equation}
{\phi _{r2}} = \min \left[ {\begin{array}{*{20}{c}} 1&,&{{{\left( {\frac{{k_r^ + }}{{30}}} \right)}^{{2
/ 3}}}} \end{array}} \right]\min \left[ {\begin{array}{*{20}{c}} 1&,&{{{\left( {\frac{{k_r^ + }}{{45}}} \right)}^{{1
/ 4}}}} \end{array}} \right]\min \left[ {\begin{array}{*{20}{c}} 1&,&{{{\left( {\frac{{k_r^ + }}{{60}}} \right)}^{{1
/ 4}}}} \end{array}} \right]
\end{equation}$$
Under transitionally rough conditions, Knopp et al(Reference Knopp, Eisfeld and Calvo15) use a linear blending function
$$\begin{equation}
{k_w} = \frac{{u_\tau ^2}}{{C_\mu ^{0.5}}}{\rm{min}}\left( {\begin{array}{*{20}{l}} 1&{\frac{{k_r^ + }}{{90}}} \end{array}} \right)
\end{equation}$$
for the k boundary condition. However, due to its first derivative discontinuity, we found that the use of Formula (16) might lead to sudden jumps of the skin friction distribution in streamwise direction. Therefore, the present authors propose a more continuous representation as expressed as
$$\begin{equation}
{k_w} = \frac{{u_\tau ^2}}{{C_\mu ^{0.5}}}{\left[ {tanh\left( {\frac{{k_r^ + }}{{60}}} \right)} \right]^2}
\end{equation}$$
2.4 Compressibility corrections
The above Menter's k–ω SST model was developed mostly for incompressible flows. To account for compressibility effects, it is necessary to re-examine the Favre-averaged equation for the turbulent kinetic energy k which can be written as
$$\begin{equation}
\bar{\rho }\frac{{Dk}}{{Dt}} = \bar{\rho }{P_k} - \bar{\rho }{\varepsilon _s} - \bar{\rho }{\varepsilon _d} + {T_k} + \overline {p^{\prime}d^{\prime\prime}} + M,
\end{equation}$$
where the terms on the right-hand side are the production term of turbulent kinetic energy
$\bar{\rho }{P_k}$
, the solenoidal (incompressible) dissipation
$\bar{\rho }{\varepsilon _s}$
, the dilatation (compressible) dissipation
$\bar{\rho }{\varepsilon _d}$
, the turbulent transport term Tk
, the pressure-dilatation term
$\overline {p^\prime d^{\prime\prime}} $
, and the mass flux variation M. The dilatation-dissipation and pressure-dilatation terms appear explicitly in k equation and directly affect the turbulence energetics. Here, we adopt the dilatation-dissipation model proposed by Sarkar et al(Reference Sarkar, Erlebacher and Hussaini32) as
$$\begin{equation}
{\varepsilon _d} = 0.6M_t^2{\varepsilon _s}
\end{equation}$$
The term Mt is the turbulent Mach number defined by Mt = (2 k)0.5/a, a being the speed of sound. The pressure-dilatation model proposed by Sarkar(Reference Sarkar33) reads
$$\begin{equation}
\overline {p^\prime d^{\prime\prime}} = 0.15\bar{\rho }M_t^2{\varepsilon _s} - 0.2\bar{\rho }{M_t}{P_k}
\end{equation}$$
In all, the SST model equations with compressibility corrections can be expressed as:
$$\begin{equation}
\frac{{\partial \left( {\bar{\rho }k} \right)}}{{\partial t}} + \frac{\partial }{{\partial {x_j}}}\left[ {\bar{\rho }{{\tilde{u}}_j}k - \left( {\mu + {\sigma _k}{\mu _t}} \right)\frac{{\partial k}}{{\partial {x_j}}}} \right] = \bar{\rho }{\tilde{\tau }_{ij}}{\tilde{S}_{ij}}\left( {1 - 0.2\bar{\rho }{M_t}} \right) - {\beta ^*}\bar{\rho }k\omega \left( {1 + 0.75M_t^2} \right)
\end{equation}$$
and
$$\begin{eqnarray}
&&\frac{{\partial \bar{\rho }\omega }}{{\partial t}} + \frac{\partial }{{\partial {x_j}}}\left[ {\bar{\rho }{{\tilde{u}}_j}\omega - \left( {\mu + {\sigma _\omega }{\mu _t}} \right)\frac{{\partial \omega }}{{\partial {x_j}}}} \right]\nonumber\\
&&\quad = {P_\omega } - \beta \bar{\rho }k{\omega ^2}\left( {1 - 0.15M_t^2} \right) + 2\left( {1 - {F_1}} \right)\frac{{\bar{\rho }{\sigma _{\omega 2}}}}{\omega }\frac{{\partial k}}{{\partial {x_j}}}\frac{{\partial \omega }}{{\partial {x_j}}}
\end{eqnarray}$$
3.0 NUMERICS
All of the simulations presented here are performed using an in-house code where the three-dimensional compressible Navier-Stokes equations are solved by using Roe's implicit, finite volume, upwind algorithm. By means of the monotone upstream-cantered schemes for conservation laws interpolation of the primitive variables, the quantity in the inviscid fluxes is obtained. The viscous flux terms are calculated by a second-order central difference.
Although the proposal described in Section 2 adopts simple Dirichlet boundary conditions, it depends on the friction velocity uτ . Its determination requires the numerical evaluation of the normal derivative ∂/∂y of the tangential velocity component, ut , at the wall. uτ is calculated with second-order approximation:
$$\begin{equation}
{u_\tau } \simeq \sqrt {\nu \frac{{2{u_{t,2}} - 0.5{u_{t,3}}}}{{2{y_2} - 0.5{y_3}}}}
\end{equation}$$
The subscripts 2 and 3 stand for the first two grid nodes away from the wall and ut ,1 = y 1 = 0.
For all test case geometries, gird independence has been achieved. The primarily created “baseline grid” is coarsened and refined in each direction by increasing and decreasing the number of grid points by a factor of 2, respectively. Figure 1 shows the computational domain and mesh for the Ligrani and Moffat flat-plate case(Reference Ligrani and Moffat40). The grid node distribution in the longitudinal (x) direction was non-uniform with clustering of nodes near the leading and trailing edge of the plate using one or two-sided stretching functions by Vinokur(Reference Vinokur41). In all the grids, 3/4 of the nodes in the longitudinal direction are located on the plate and 1/4 are located upstream of the plate. The number of grid nodes in the longitudinal direction is fixed for the different Reynolds numbers.
Figure 1. The computational domain and baseline mesh for the Ligrani & Moffat flat-plate case(Reference Ligrani and Moffat40).
In the wall-normal (y) direction, the number of grid nodes depends on the Reynolds number and the near-wall grid line distances. Especially, the grids for the direct application of the no-slip condition are obtained from a one-sided stretching function using different values of the stretching parameter per Reynolds number. For all cases, the first cell y + value is of 0.3 with 100 nodes inside the boundary layer. The investigation by Knopp et al(Reference Knopp, Eisfeld and Calvo15) led to the choice of these values.
In all cases, steady-state solutions have been obtained on both baseline and refined meshes, which show negligible difference and are therefore judged to be mesh independent. All cases are run to full convergence, determined based on a drop in residuals of typically five orders of magnitude, as well as a flattening of all residuals indicating that machine accuracy has been reached.
4.0 MODEL VALIDATION
The new roughness proposal has been validated for a reasonably wide range of rough-wall cases involving incompressible turbulent flows past flat plates and pipes under zero pressure gradients (PG). The test cases used by Knopp et al are adopted here to ensure that the present modification either improves or doesn't contaminate their original model. Results labelled as “DLR” correspond to those predictions using the roughness correction, Formula (16), proposed by Knopp et al(Reference Knopp, Eisfeld and Calvo15). Those of the present model, Formula (17), are labelled “New” or not labelled. It is noted that the curvature effect is neglected here and will be considered in further works.
4.1 Experiments by Ligrani and Moffat(Reference Ligrani and Moffat40)
The test case by Ligrani and Moffat(Reference Ligrani and Moffat40) is a flat-plate turbulent boundary layer flow over spherical roughness elements. The flat plate length L = 5 m, the kinetic viscosity ν = 1.5E-5 m2/s. The roughness height is held constant, with the corresponding equivalent sand grain roughness size kr of 0.79 mm, as confirmed in an earlier study using fully rough velocity-profiles information(Reference Ligrani and Moffat40). By altering the free-stream velocity, transitionally rough conditions are obtained, as listed in Table 1.
Table 1 Flow conditions in the experiment by Ligrani and Moffat(Reference Ligrani and Moffat40). The values for k + r are measured at x/L = 0.356
The predicted velocity profiles in the logarithmic region, as shown in Fig. 2, give good agreement with the theoretical relation, Equation (8), in transitionally rough conditions. It is also seen that with the increase of dimensionless roughness heights, the velocity shift of the logarithmic profile increases. Figure 3 shows the computed skin friction coefficient, Cf , with comparison to experimental data. It indicates that compared to DLR results, ours are closer to the measurements. Note that the uncertainty regarding the experimental data for Cf is about ±10%(Reference Ligrani and Moffat40). For this case, the new proposal can achieve reasonable results for both Cf and ∆u, which was thought to be difficult for the k–ω type models with roughness corrections(Reference Boyle and Stripf24).
Figure 2. Calculated velocity profiles at x/L = 0.356, for different dimensionless roughness heights measured at x/L = 1.
Figure 3. Skin friction (Cf) distributions predicted by the present and the DLR15 models, with comparison to experimental data(Reference Ligrani and Moffat40).
4.2 Experiments by Nikuradse(Reference Nikuradse39)
Here we consider the experiments by Nikuradse(Reference Nikuradse39) for fully developed flows in pipes of various roughness heights. Tables 2 and 3 provide a summary of the parameter variations for high and medium Re (Reynolds number) regimes, respectively. For R/kr larger than 252 in Table 2 as well as R/kr larger than 30.6 in Table 3, transitional roughness values can be observed.
Table 2 Overview of flow parameters at high Reynolds numbers in the Nikuradse experiment(Reference Nikuradse39). R denotes the pipe radius
Table 3 Overview of flow parameters at medium Reynolds numbers in the Nikuradse experiment( Reference Nikuradse 39 ) . R denotes the pipe radius
The tables also compare the present and DLR predictions for the roughness Reynolds number k + r , at x/L = 1, with the experimental data. It is noted that the prediction accuracy of Cf directly depends on the agreement in k + r . Good agreement is achieved between calculations and measurements for all the cases. Figure 4 depicts the predicted velocity profiles in the logarithmic region. It is seen that the predicted shift of the log-law matches the measurements fairly well, except for R/kr = 507 at medium Re and R/kr = 60 at low Re. Nevertheless, the new model gives almost identical predictions with the DLR results.
Figure 4. Velocity profiles at medium (left) and low (right) Re, predicted by the present and the DLR15 models, with comparison to experimental data(Reference Nikuradse39).
4.3 Experiment by Blanchard(Reference Blanchard42)
Blanchard(Reference Blanchard42) measured zero pressure gradient flow over a rough surface with the equivalent sand grain roughness of 0.85 mm. With reference to the work by Knopp et al(Reference Knopp, Eisfeld and Calvo15), the inflow velocity is set to 45 m/s and the reference length L = 0.8 m. Figure 5 shows the computed skin friction coefficient, Cf , with comparison to experimental data. It is seen that although there are some gaps between calculations and measurements, the consistent distribution trends.
Figure 5. Skin friction (Cf ) distributions predicted by the present and the DLR(Reference Knopp, Eisfeld and Calvo15) models, with comparison to experimental data(Reference Blanchard42).
4.4 Experiments by Hosni et al(Reference Hosni, Coleman and Taylor43, Reference Hosni, Coleman, Garner and Taylor44)
Hosni et al(Reference Hosni, Coleman and Taylor43,Reference Hosni, Coleman, Garner and Taylor44) experimentally studied the turbulent boundary layer flow over a rough surface of length L = 2.4 m composed of hemispheres of diameter l 0 = 1.27 mm. The equivalent sand grain roughness is of 1.09 mm. The investigation by Knopp et al(Reference Knopp, Eisfeld and Calvo15) led to the choice of this value. The test case MSU1 with inflow velocity of 12 m/s, gives the transitional roughness value, k + r = 40, at x = 2 m, while the case MSU2 with inflow velocity of 58 m/s, shows the fully rough condition, k + r = 200, at x = 2 m.
Figure 6 depicts the predicted skin friction coefficient, Cf , with comparison to measurements. For both MSU1 and MSU2 cases, the present and the DLR models give similar Cf distributions. It is noted that for both cases the uncertainty in Cf is estimated to be 10%-12%(Reference Hosni, Coleman and Taylor43,Reference Hosni, Coleman, Garner and Taylor44) .
Figure 6. Skin friction (Cf ) distributions For both MSU1 and MSU2 cases predicted by the present and the DLR(Reference Knopp, Eisfeld and Calvo15) models, with comparison to experimental data(Reference Hosni, Coleman and Taylor43,Reference Hosni, Coleman, Garner and Taylor44) .
5.0 RESULTS OF SUPERSONIC/HYPERSONIC FLAT-PLATE FLOW CASES
In this section, our validated method is applied to the supersonic/hypersonic flow past a flat plate with transitionally and fully rough surfaces. The flow conditions considered here refer to the typical atmospheric flight conditions with several flight height levels as well as Mach numbers. At first, we investigate the roughness effect on drag. Then, the effects of flight heights and Mach numbers on drag are studied within the roughness heights measured in real flight. The adiabatic wall boundary conditions are imposed here.
5.1 Roughness effect
The numerical experiments are set up with the inflow conditions as M ∞ = 6, H = 26 km (Re ∞ = 4.3 E6/m, T ∞ = 222.5 K), 0° angle-of-attack. The roughness heights simulated are 20 μm, 200 μm, 700 μm and 1,000 μm. Within the baseline grid, the streamwise and wall-normal directions are resolved by 241 and 161 points, respectively.
Figure 7 shows the roughness height effect on the computed velocity shift of the logarithmic profile, which is pretty weak compared with the low-speed case. Especially, the velocity profile at kr = 20 μm almost superpose with the one on the smooth wall. Therefore, dimensionless roughness heights need to be calculated, as given in Fig. 8. It is seen that the kr values of 20 μm and 200 μm correspond to the hydrodynamically smooth walls (k + r < 3.5), while the others to the transitional roughness regime (3.5 < k + r < 68)as defined by Nikuradse(Reference Nikuradse39) for incompressible channel flows. This can mostly be attributed to the dramatic increase of boundary layer thickness, δ, in hypersonic flows, as estimated by(Reference Anderson45)
$$\begin{equation}
\frac{\delta }{x} \propto \frac{{M_{}^2}}{{{\mathop{\rm Re}\nolimits} _x^2}}
\end{equation}$$
Figure 7. Calculated velocity profiles at x/L = 0.356, for different dimensionless roughness heights, the right zooms in the left.
Figure 8. Dimensionless roughness height distributions.
The δ variations with different roughness heights (kr ) are given in Table 4. The boundary layer edge is defined as the location where the local total enthalpy achieves 1.005 times the stagnation enthalpy at a certain profile. It verifies that the roughness height is pretty small compared with boundary layer thickness, though its absolute value is already beyond the one for real aircraft surface that is usually smaller than 100 μm. Moreover, Table 4 indicates that roughness causes additional increase of δ.
Table 4 Variations of boundary layer thickness (δ) with different roughness heights (kr ) at x/L = 0.356
Figure 9 depicts the skin friction distribution along the flat plate. It rises with the increase of the roughness height. Since there is no pressure contributed in it, the drag can be integrated from Fig. 9, as given in Table 5. For comparison, it is more reasonable to consider the relative rough-wall drag derivation from the smooth-wall case, as defined as
$$\begin{equation}
\begin{array}{*{20}{l}} {\Delta C_d^{^\% }\frac{{{C_{d,{\rm{rough}}}} - {C_{d,{\rm{smooth}}}}}}{{{C_{d,{\rm{smooth}}}}}},}&{{C_d} = \frac{D}{{{\rho _\infty }U_\infty ^2A}}} \end{array}
\end{equation}$$
Figure 9. Predicted skin friction (Cf ) distributions for different roughness heights.
Table 5 Variations of the drag coefficient (Cd ) as well as the relative rough-wall Cd derivation from the smooth-wall case (∆Cd %), with the equivalent sand grain height kr
Here D is the drag force and A the reference area of 1 mReference Fu and Wang2. It is seen from Table 5 that for the selected inflow conditions, the maximum drag increase of 4.33% is achieved at kr = 1,000 μm. In general, the roughness effect on drag is pretty weak in such hypersonic flat-plate flows.
5.2 Effects of flight heights and Mach numbers
Actually, the roughness height value for real aircraft surface (the nose region excluded) are usually smaller than 100 μm. Therefore, the roughness heights of 20 μm and 200 μm are selected to study the effects of flight heights and Mach numbers on drag with surface roughness. The Mach numbers simulated are 4, 5, 6 and 7; the flight heights considered are 20, 24, 28, 32 and 36 km. Table 6 shows the unit Reynolds numbers for different flight heights at M ∞ = 4. It is seen that with the flight height increase, the air density decrease is more intensive than the dynamic viscosity decrease, resulting in the reduction of unit Reynolds number. The computational meshes are generated corresponding to different Reynolds numbers, as discussed in Section 3.
Table 6 Variations of flow parameters with flight heights at x = 1 m, M ∞ = 4 and kr = 100 μm
Figure 10(a) depicts the flight height effect on the computed velocity shift of the logarithmic profile at x = 1 m, M ∞ = 4 and kr = 100 μm. It is seen that all curves superpose in the near-wall region. If the characteristic length and velocity are chosen as L and U∞ , respectively, the dimensionless velocity profiles (see Fig. 10(b)) shows that with the flight height decrease (Re increase), the velocity increases at the same wall distance.
Figure 10. Calculated velocity profiles, at x = 1 m, for different flight heights. (a) u += u/uτ , y += yuτ/ν; (b) U = u/U∞ , Y = y/L. M = 4, kr =100 μm.
Figure 11(a) plots the near-wall behaviour for turbulent kinetic energy, k, at x = 1 m and M ∞ = 4. It demonstrates that the new roughness model remedies the near-wall gradient of k appearing for the Wilcox roughness modification(Reference Wilcox19). It is seen that the flight height change has little effect on k in both viscous sublayer and log-layer regions, which can be attributed to the corresponding small k + r given in Table 6. In contrast, the flight height effect on skin friction distribution is visible, as shown in Fig. 11(b).
Figure 11. Predicted profiles of dimensionless turbulent kinetic energy (k +) at x = 1 m (a) and skin friction (Cf ) distributions (b) for different k + r that correspond to different flight heights at M ∞ = 4. k + = 2k/uτ 2.
Figure 12 compares flow parameters concerned, at M ∞ = 7, for different flight heights and for different roughness heights. Duo to k + r inside the hydrodynamically smooth-wall regime for all cases (see Fig. 12(d)), slight difference can be observed between the cases at the same flight height amplitude. With the flight height decrease, significant changes in both velocity and turbulent kinetic energy profiles occur, while the skin friction distributions shift downwards. It can be concluded that the flight height effect is much stronger than the effect of roughness heights considered in the present study.
Figure 12. Flow parameter predictions for different flight heights and for different dimensionless roughness heights at M ∞ = 7. Profiles of dimensionless velocity (a) and turbulent kinetic energy (b) at x = 1 m; skin friction (c) and dimensionless roughness height (d) distributions. k + = 2k/uτ 2.
We then investigate the Mach number effect by maintaining the same flight height. Table 7 shows the flow parameters at H = 36 km and kr = 100 μm. The variation of Mach number is purely due to the change of inflow velocity. It is seen from Figs 13 and 14(a) that with the Mach number increase, the logarithmic profile shifts upwards and the skin friction decreases. Figure 14(b) illustrates that the turbulent kinetic energy trends to zero in the near-wall region, because of the very small k + r for these cases (see Table 7).
Table 7 Variations of flow parameters with Mach numbers at x = 1 m, H = 36 km and kr = 100 μm
Figure 13. Calculated velocity profiles at x = 1 m, for different Mach numbers, the right zooms in the left. H = 36 km and kr = 100 μm.
Figure 14. Predicted skin friction (Cf ) distributions (a) and profiles of dimensionless turbulent kinetic energy (k +) at x = 1 m (b) for different Mach numbers. k + = 2k/uτ 2, H = 36 km and kr = 100 μm.
By integrating the skin friction distributions along the flat plate, we obtain the drag coefficients Cd for the 60 cases in total, as given in Table 8. It can be concluded that Cd rises with the increase of flight height H (at constant M ∞) while with the increase of Mach number M ∞ (at constant H); the derivation between rough-wall and smooth-wall conditions, ∆Cd %, achieves the maximum value of 0.79% for the case with M ∞ = 7, H = 36 km and kr = 100 μm. It is noted that greater resolution at Mach numbers and flight heights values can be interpolated to those already listed in Table 8.
Table 8 Variations of drag coefficient (Cd ) as well as its relative derivation from the smooth-wall case (∆Cd %), with the dimensionless roughness height kr and the flight height H, at M ∞ = 4 (a), 5 (b), 6 (c) and 7 (d)
6.0 CONCLUSION
In this study, a new extension for the SST k–ω turbulence model to account for surface roughness as well as flow compressibility effects has been presented which allows for the simulation of supersonic/hypersonic flows over rough surfaces at the same grid resolution requirements as for smooth walls. The new roughness modification gives slightly improved predictions in skin friction for low-speed cases compared to the roughness extension by Knopp et al(Reference Knopp, Eisfeld and Calvo15).
Then, using the current method, aerodynamics of supersonic/hypersonic flows along a flat plate of finite length is investigated numerically with Mach number up to 7 and flight height up to 36 km. It is found that the drag coefficient derivation between rough-wall and smooth-wall conditions, achieves the maximum value of at kr = 1,000 μm, M ∞ = 6, H = 26 km; within the roughness measured in real flight (kr < 100 μm), the roughness height change has little effect on drag compared to the variations of either flight heights or Mach numbers. This can mostly be attributed to the dramatic increase of boundary layer thickness in hypersonic flows. Hence, the dimensionless roughness height becomes pretty small, approaching the hydrodynamically smooth wall condition.
ACKNOWLEDGEMENTS
This work has been funded by the National Key Basic Research Program of China (2014CB744801), the National Natural Science Foundation of China for the Grants 11572177, 11572176, 51376106 & 11272183, and the Tsinghua University initiative Scientific Research Program (2014z21020).
