An extension of Thwaites’ method for turbulent boundary layers

Often in the early stages of the engineering design processes, viscous effects are estimated using boundary layer integral methods. Common elements in many such engineering flows are that boundary layers are turbulent and subjected to strong pressure gradients - both favorable and adverse. A pressure gradient effect can alter the transition location, or even cause the flow to separate, leading to significant changes in the flow patterns. Clauser (1954, 1956) studied turbulent boundary layers and defined a parameter, $\beta$ that quantifies the relative strength of the pressure gradient in relation to the skin friction across the boundary layer,

$$\beta=\frac{\delta^{*}}{\rho u_{\tau}^{2}}\frac{dP}{ds}$$

(1)

where $\delta^{*}$ is the displacement thickness, $dP/ds$ is the pressure gradient and $u_{\tau}$ is the friction velocity. This parameter has since been used in the analyses of equilibrium and non-equilibrium boundary layers. In fact, in the presence of strong, and/or prolonged favorable pressure gradients, the flow may relaminarize (Sreenivasan, 1982); thus regaining an ordered structure with a reduction in Reynolds stresses. On the other hand, in the presence of adverse pressure gradients, it is possible for the flow to eventually separate (Simpson, 1981, 1983), often associated with periods of rapid growth of the boundary layer. An adverse pressure gradient significantly energizes the outer flow (Knopp et al., 2021), and increases the momentum deficit compared to the freestream value. In fact, the effect of a pressure gradient is not only local but strongly dependent on its history (Bobke et al., 2017; Parthasarathy & Saxton-Fox, 2023; Prakash et al., 2023). Wei et al. (2017) derived the relationships between the turbulent shape factor and the quantities at the edge of the boundary layer of an adverse pressure gradient boundary layer and demonstrated that the boundary layer thicknesses, its growth rate, and the resulting skin friction are affected by upstream effects. Other experimental and computational investigations include the works of Harun et al. (2013); Skote et al. (1998); Knopp et al. (2015); Kitsios et al. (2017); Romero et al. (2022); these studies have indicated that the adverse pressure gradient energizes both the large-scale and the smaller-scale structures, while also modulating the impact of the large scales on the smaller scales.

In the limit of thin boundary layers, the growth rate of the momentum thickness is related to the pressure gradient as follows (via the Karman momentum integral, see White (2008) for more details),

$$\frac{d\theta}{ds}=\frac{C_{f}}{2}+(2+H)\frac{\theta}{U^{2}_{e}}\frac{dP_{e}}{ds}$$

(2)

where $\theta$ is the momentum thickness, and $C_{f}$ is the skin-friction along the streamwise coordinate. $H$ is the shape factor which is the ratio of the displacement and the momentum thickness lengths of the boundary layer. However, this equation, in itself, is not particularly useful as a predictive tool and requires the knowledge of the skin friction and shape function distributions, which may be unknown a priori.

For laminar boundary layers, Thwaites (1949), and later Curle & Skan (1957), and Dey & Narasimha (1990) fitted data for various flows, including those with pressure gradients and flows undergoing separation, to derive an approximate expression for the growth of the momentum thickness for a given inviscid (freestream) velocity distribution. This analysis coupled with the Falkner & Skan (1931) analysis in the viscous layer provides a complete description of the mean velocity profile of the corresponding laminar flow. However, a similar analysis for turbulent boundary layers is more complex and lacking. Weber (1978) developed a semi-empirical method for calculating the thickness of turbulent boundary layers with arbitrary pressure gradients for the AFOSR-IFP-Stanford conference (Coles & Hirst, 1968; Kline et al., 1969). However, his approach relied on implied correlations between the skin friction, shape factor, and an inviscid pressure gradient parameter, some of which may be unknown, especially in complex engineering flows. Das (1987), Das & White (1986), Kalkhoran & Wilson (1986) developed integral approaches for measuring skin friction in turbulent boundary layers. A drawback of Das (1987)’s and Das & White (1986)’s approaches included the use of multiple correlations while also not accounting for history effects (through appropriate streamwise derivatives). Kalkhoran & Wilson (1986) made equilibrium assumptions for the wake and the pressure gradient even for non-equilibrium boundary layers. A comprehensive review of these approaches can be found in Das et al. (2004).

From a practical standpoint, the effects associated with pressure gradients pose significant modeling challenges for near-wall models that are employed to predict wall-shear stress in scale-resolving methods. In the past, various modeling efforts have been pursued (Balaras et al., 1996; Wang & Moin, 2002; Duprat et al., 2011; Hickel et al., 2013; Park & Moin, 2014) to include the non-equilibrium effects due to pressure gradients. However, these models have not shown significant differences in their predictive performance in comparison to the equilibrium wall model (Cabot & Moin, 2000; Kawai & Larsson, 2012) across a wide range of flows. Park & Moin (2014) developed a non-equilibrium wall model while including the full two-dimensional RANS equations, which adds to the computational complexity. This model was shown to have limited improvements over the equilibrium wall model in the flow past a hump (Park, 2017). Duprat et al. (2011) did not account for convective terms in their non-equilibrium wall model while adjusting the RANS eddy viscosity to account for pressure gradients. Hickel et al. (2013) accounts for both convective, and pressure gradient terms, but inconsistently uses an equilibrium flow-based eddy-viscosity model.

It is apparent from the aforementioned efforts that further developments in the understanding of the growth of turbulent boundary layers under the action of a pressure gradient are needed for improving potential flow solvers and for near-wall models to account for pressure gradient effects. The objective of this paper is to investigate whether a simple extension to the original Thwaites method for turbulent flows can be derived based on the assimilation of a limited set of high-fidelity simulation data. The rest of this article is organized as follows, Section II describes the original Thwaites method for laminar flows, and Section III provides the proposed extension to turbulent flows. Section IV describes the simulation database used in this study. Section V describes the model coefficients and their validation in several different flows. Section VI discusses an application of the model in terms of the sensitivity of the separation parameter to flow perturbations. Concluding remarks are offered in Section VII.

In this section, a brief summary of the method of Thwaites (1949) for approximating the momentum thickness of laminar boundary layers is provided. Recall, that for thin boundary layers in an incompressible flow, the Von Karman integral equation ((White, 2008)) is written as,

$$\frac{d\theta}{ds}+(2+H)\frac{\theta}{U_{e}}\frac{dU_{e}}{ds}=\frac{C_{f}}{2}$$

(3)

where in the $n-s$ coordinate system, $s$ is the streamwise coordinate, and $n$ is the wall-normal coordinate. $\delta^{*}$ and $\theta$ are the displacement and momentum thickness lengths of the boundary layer, respectively, and $H=\delta^{*}/\theta$ is the shape factor. The following relations also hold under the thin boundary layer approximation,

$$\nu\frac{\partial^{2}U}{\partial n^{2}}\bigg{|}_{n=0}=\frac{\partial P_{e}}{\partial s}=-U_{e}\frac{\partial U_{e}}{\partial s}$$

(4)

where $P_{e}$ and $U_{e}$ are the values of the pressure, and the streamwise velocity at the edge of the boundary layer, and $n=0$ refers to the wall. Using these relations, Thwaites derived that

$$2Re_{\theta}\frac{d\theta}{ds}=2\left(\left(2+H\right)m+\chi\right)=L$$

(5)

where $m$ is the Holstein-Bohlen pressure gradient parameter, $m=\frac{\theta^{2}}{U_{e}}\frac{\partial^{2}U}{\partial n^{2}}\big{|}_{n=0}$, and $\chi=\frac{\theta}{U_{e}}\frac{\partial U}{\partial n}\big{|}_{n=0}$. It should be noted that in deriving Equation 5, the assumption of a laminar flow has not yet been invoked. For laminar boundary layers, however, Thwaites further postulated that $L=L(m)$, or that the effect of the pressure gradient primarily determines the growth rate of a laminar boundary layer. Specifically, the linear fit proposed in Thwaites’ original work is

$$L\approx L(m)^{Thwaites}\approx 0.45+6m$$

(6)

Rearranging the model fit, an approximate closed-form expression of the momentum thickness can be written as,

$$\theta^{2}(s)\approx 0.45\frac{\nu}{U^{6}_{e}}\int^{s}_{0}U^{5}_{e}(r)dr$$

(7)

Further, Thwaites fitted experimental data to suggest that a laminar flow may separate when $m\approx 0.09$. Other similar empirical thresholds have been proposed by Stratford (1959), and by Curle & Skan (1957). More details can be found in Horton (1968). Curle & Skan (1957), and later Dey & Narasimha (1990) suggested slightly different model coefficients that resulted in a better agreement with the Blasius and Falkner-Skan type solutions. These methods have been used in potential flow solvers to account for viscous effects in laminar flows.

Unlike in laminar flows, the inner and outer layers of a turbulent boundary layer respond differently to external pressure gradients. Further, a thin internal layer develops near the wall in turbulent flows due to the pressure gradient (Uzun & Malik, 2022; Parthasarathy & Saxton-Fox, 2023), and the amplitude of the small scales near the wall gets modulated due to the larger ones (Knopp et al., 2015; Romero et al., 2022). A general method for predicting the momentum thickness, $\theta$, similar to that described in the previous section is hence lacking for general, non-equilibrium turbulent flows.

Consider a zero pressure gradient turbulent boundary layer at a high Reynolds number. The growth rate of the momentum thickness of a turbulent boundary layer is higher than that of an equivalent Blasius laminar boundary layer ($\theta/s\sim Re_{s}^{-1/2}$). Reconciling this fact with the fact that Equation 5 still holds for turbulent flows, the approximation made by Thwaites in Equation 6, $L=L(m)$ must not hold for turbulent boundary layers. Since a boundary layer at a higher Reynolds number is less sensitive to the effects of pressure gradients (Vinuesa et al., 2018), the growth rate can not depend only on the pressure gradient (as it did in the laminar case). Therefore, it is proposed that the growth rate of the momentum thickness may also depend directly on the Reynolds number. Equation 5 can be re-expressed as,

$$2\frac{d\theta}{ds}=2\frac{(2+H)m}{Re_{\theta}}+\frac{2\chi}{Re_{\theta}}.$$

(8)

A Taylor’s series expansion of $2{d\theta}/{ds}$ in terms of $m/Re_{\theta}$ and $1/Re_{\theta}$ gives,

$$2\frac{d\theta}{ds}=C_{Re,\infty}+C_{m}\frac{m}{Re_{\theta}}+C_{c}\frac{1}{Re_{\theta}}+C_{m,2}\left(\frac{m}{Re_{\theta}}\right)^{2}+C_{c,2}\left(\frac{1}{Re_{\theta}}\right)^{2}+C_{m,c}\left(\frac{1}{Re_{\theta}}\right)\left(\frac{m}{Re_{\theta}}\right)+\cdots,$$

(9)

It has been well established that for turbulent boundary layers, the presence of strong pressure gradients leads to “history effects”, or that two boundary layers experiencing different upstream pressure gradients may respond differently to the action of the same pressure gradient even when the Reynolds number ($Re_{\theta}$ or $Re_{\tau}$) is held constant (Bobke et al., 2017; Devenport & Lowe, 2022; Vinuesa et al., 2017; Nickels, 2004). This is because under the influence of the pressure gradient, the near-wall, and outer layer scales of motion respond at different spatial, and temporal scales (Knopp et al., 2021). To account for many of these effects, a comprehensive database of turbulent boundary layers is used in this work. All of the data used is from numerical simulations for ease of accessing continuous data, for computing streamwise derivatives and integrals. Relevant details of the datasets are given below. Table 1 summarizes the variation in the Reynolds number and the Clauser parameter for the boundary layers considered. It is noteworthy that Goc et al. (2023) reported similar skin-friction Reynolds numbers across the wing of a transonic aircraft model as those observed in this database.

• •

  ZPG Eitel: Eitel-Amor et al. (2014) performed wall-resolved large eddy simulations (LES) of a spatially developing zero pressure-gradient boundary layer up to $Re_{\theta}=8300$.

• •

  APG Bobke: Bobke et al. (2017) performed wall-resolved LES of five different boundary layers at varying Reynolds numbers and strength of pressure gradients. Two of these five boundary layers were maintained at nearly constant Clauser parameters, $\beta=1,\;\mathrm{and}\;2$ while the other three boundary layers followed a power law for the streamwise velocity at the edge of the boundary layer, $U_{e}\sim(x-x_{0})^{m}$ with $m=-0.13,-0.16,\;\mathrm{and}\;-0.18$.

• •

  Airfoils: Vinuesa et al. (2018) and Tanarro et al. (2020) simulated the flow over two different NACA airfoils, namely NACA 4412 and 0012 respectively at varying chord-length-based Reynolds number and angles of attack. In this work, we consider three particular cases, NACA 0012 airfoil at $Re_{c}=0.4\mathrm{Mil}$, $\alpha=0.0^{o}$ deg., and NACA 4412 airfoil at $Re_{c}=0.1\mathrm{Mil},\;1\mathrm{Mil}$, $\alpha=5.0^{o}$ deg.

• •

  FPG/APG SBSE: The problem of predicting smooth body separation using both Reynolds averaged methods (RANS) and scale-resolving methods in CFD has remained a challenge. Uzun & Malik (2022) performed a quasi-DNS of the flow over a smooth Gaussian bump as it experienced both a favorable pressure gradient, followed by an adverse pressure gradient, and thereafter a turbulent separation. These results were in agreement with the experiments of Williams et al. (2020). Later, Agrawal et al. (2022); Whitmore et al. (2022) performed wall-modeled large-eddy simulations of this flow and accurately predicted the flow separation.

• •

  FPG/APG BJ: In aerospace applications, such as the flow over the transonic Common Research Model (CRM) aircraft (Rivers & Dittberner, 2010), weak compressibility effects affect the flow pattern due to the shock-formation, and the consequent flow separation at high Mach numbers. Uzun & Malik (2019) performed high-fidelity simulations of the experiments of Bachalo & Johnson (1986) that studied the flow over a hump on a cylinder experiencing shock and adverse pressure-gradient induced flow separation.

This work presents an extension of the method of Thwaites for determining the momentum thickness for a turbulent boundary layer under the action of pressure gradients. In the limit of large Reynolds numbers, a linear model is hypothesized, such that the Reynolds number and pressure gradient dependence on the growth rate of the momentum thickness can be superposed. A fit for the model coefficients is found using recent high-fidelity simulation data for various boundary layers ranging from Reynolds number, $Re_{\tau}\approx[100,4000]$ and Clauser parameter, $\beta,\approx[-6,10]$. The model fit predicts the growth of the momentum thickness, even up to the point of separation. The fits utilized in the model can then be used to derive a condition for imminent separation at high Reynolds numbers, with similar threshold values as compared to previously proposed values based on empirical data. An application of this extension of Thwaites method is demonstrated by estimating the importance of history effects on separation. Specifically, the importance of upstream history effects for the Boeing speed-bump flow is quantified. For a transonic bump flow, it is shown that only local perturbations just upstream of the shock location cause significant changes to the separation location.