Applicability of Frozen-Viscosity Models of Unsteady Wall Shear Stress

The simulation of unsteady fluid flows in extensive pipe or duct networks such as water supply, sewerage, oil and gas lines, and railway tunnels is nearly always undertaken using one-dimensional $(1\mathrm{D})\{x,t\}$ methods in which no explicit account is taken of lateral variations in a cross section. This is a practical necessity because the use of two-dimensional (2D) and three-dimensional (3D) methods would be prohibitively time-consuming. Fortunately, $1\mathrm{D}\{x,t\}$ approximations are sufficiently accurate for most practical purposes and, indeed, they can sometimes yield more accurate predictions than their 2D/3D counterparts. This may be so, for instance, when empirical resistance coefficients are known with good accuracy in pipes with irregular surface finishes. Nevertheless, there are also cases when the consequences of 2D/3D flow phenomena are of special importance. In such cases, it can be difficult to develop $1\mathrm{D}\{x,t\}$ approximations that provide adequate representations of the phenomena. One example is the detailed response of flows after a strong change of acceleration or deceleration. In particular, the response of wall shear stresses depends strongly on the relationship between turbulence timescales and bulk-flow timescales (e.g., Ghidaoui et al. 2002).

In this paper, attention focuses on a special case of unsteady flow in which an initially steady flow in a pipe is suddenly accelerated. The ultimate aim is to enhance the capabilities of $1\mathrm{D}\{x,t\}$ methods of analyzing such flows, but the simulations presented herein are focused on radial changes in time and are obtained using $2\mathrm{D}\{r,x,t\}$ or, when appropriate, $1\mathrm{D}\{r,t\}$ analyses. The $1\mathrm{D}\{x,t\}$ methods are needed by analysts wishing to study unsteady flows in extensive networks of pipes, and so on. Evidence is drawn from direct numerical simulation (DNS) studies that have been used elsewhere to study the flow in detail and to validate Reynolds-averaged Navier Stokes (RANS) models based on assumed models of turbulence. The DNS simulations are possible only at relatively low Reynolds numbers, but the RANS models can be used for much higher Reynolds numbers. They are used herein to investigate the overall flow behavior and to validate underlying assumptions in the methods used to enhance $1\mathrm{D}\{x,t\}$ approximations. Particular attention is paid to sustained periods of flow in which the turbulence structure responds much more slowly than the bulk velocity of flow.

Herein, attention is focused on accelerating flows of an incompressible fluid in a straight pipe of circular cross section. The presented results are applicable only in regions sufficiently far from boundaries for there to be negligible axial variation of the kinematic conditions. There are radial variations (in time and space), but no axial variations. Furthermore, for ease of interpretation, the accelerations are either (1) constant or (2) constant and then zero.

For completeness, it is useful to clarify the intended interpretation of the word acceleration. Except where the context requires otherwise, the word is used herein in a physical sense, implying an increasing velocity amplitude rather than in a strictly mathematical sense (in which, for example, a change from $-5$ to $-3\mathrm{m}/\mathrm{s}$ would be regarded as positive acceleration). The illustrative examples are chosen in such a manner as to avoid ambiguity in this respect and, in particular, to avoid flow reversals. Physically, the importance of the phenomena studied herein tends to be greatest when the Reynolds number of the initial condition is small so the focus of the paper is on flows accelerating from small to larger Reynolds numbers. Significant differences can exist between the behaviors of accelerating and decelerating flows (He et al. 2008; Ariyaratne et al. 2010), partly because the velocity profile in a decelerating flow is less stable than that in an accelerating flow (as is also the case for spatial accelerations and decelerations).

The response of turbulence and other flow properties to a suddenly imposed acceleration is described in detail in the following section using computational fluid dynamics (CFD-RANS) simulations based on the Launder-Sharma low–Reynolds number model of turbulence (Launder and Sharma 1974). Comparisons are then made with equivalent predictions based on other models of turbulence. These show strong differences in important aspects, but also show close agreement during early stages of the acceleration. The focus thereafter is on these early stages, and particular attention is given to the behavior of the wall shear stress, for which $1\mathrm{D}\{x,t\}$ models are ultimately required in engineering practice. It is shown that the use of prescribed radial distributions of the effective kinematic viscosity of the flow is especially well-suited to enabling certain consequences of inherently 2D behavior to be encapsulated in $1\mathrm{D}\{x,t\}$ outcomes. The true viscosity distributions cannot be known a priori, so the possibility of doing this with simplified distributions is explored. Simple applications of the $1\mathrm{D}\{x,t\}$ outcomes are presented, and implications for the extension of the analysis to more general flows are considered. The paper closes with concise conclusions.

Fig. 1 shows successive profiles of axial, local-mean components of velocity and shear stress in an axially uniform turbulent flow accelerating at a constant rate from an initially steady, smooth-walled pipe flow. The shear stresses are sums of the viscous and turbulent contributions, calculated from $\rho (\nu +{\nu }_t)\partial U/\partial y$, where $\nu$ and ${\nu }_t$ are the molecular and turbulent eddy viscosities, respectively. The figure shows predictions obtained using CFD based on the Launder-Sharma (LS) low–Reynolds number model of turbulence. The $1\mathrm{D}\{r,t\}$ numerical method is described in detail by He et al. (2008) and Ariyaratne et al. (2010) and needs not be restated here. It is shown in the next section that the detailed outcome is sensitive to the choice of turbulence model, but that is unimportant at this stage of the paper because attention is focused on qualitative behavior, not fine detail. More important, the sensitivity is small during the early time period when detailed use is to be made of the predictions. To facilitate comparisons with He et al. (2008), the fluid is water ($\rho =\mathrm{1,000}{\mathrm{kg}/\mathrm{m}}^3$, $\nu ={10}^{-6}{\mathrm{m}}^2/\mathrm{s}$), the pipe diameter is 50 mm, and the Reynolds number of the initial steady flow is $\mathsf{R}=\mathrm{5,000}$. At the instant $t=0$, the flow instantaneously begins to accelerate at a constant rate of ${0.9\mathrm{m}/\mathrm{s}}^2$, reaching $\mathsf{R}=\mathrm{150,000}$ at approximately $t=3.2\mathrm{s}$.

In the early stages of the acceleration (before $R\approx \mathrm{40,000}$ in the figure), the shear stress remains almost unchanged over most of the flow cross section except close to the wall (i.e., $y/R\ll 1$, where $y$ is the distance from the wall and $R$ is the pipe radius). During this period, the velocity increases almost uniformly over most of the cross section. After further acceleration ($R\approx \mathrm{80,000}$ in the figure), the shear stress has changed strongly over more than half of the radius, and the velocity profile differs markedly from its original shape. In the core of the pipe, however, the shear stress remains almost unchanged, and the velocity has continued to increase almost uniformly. Subsequently ($R\approx \mathrm{120,000}$), the conditions throughout the cross section are qualitatively similar to those in a steady flow. For instance, the shear stress distribution is approximately linear.

Fig. 2 shows corresponding histories of the velocity and shear stress at successively increased distances from the wall. In each case, there is clear evidence of a phenomenon propagating radially inward from the wall to the axis. For example, at small times, Fig. 2(a) shows the velocity increasing at almost the same rate at each depicted location. The histories then deviate from this trend, first at locations close to the wall and then at successively greater distances. At some locations, the velocity reduces when the disturbance is first encountered. At later times, however, it evolves smoothly and continuously at all locations. Also, there is a qualitative difference between the response at $y=R/16$ and the responses at greater distances.

The existence of a radially propagating phenomenon can also be seen in the shear stress histories [Fig. 2(b)], and this shows further evidence of the origin of the phenomenon. By inspection, a strong change is seen at $y=R/16$ (just after $t=1.0\mathrm{s}$) before it is seen at either $y=0$ or $y=R/8$ (at about $t=1.2\mathrm{s}$). That is, the disturbance originates at a finite distance from the wall, not at the wall itself. This is consistent with conclusions deduced by Laufer (1954) from hot wire measurements of air flow in a 250 mm diameter pipe. It is shown experimentally by He and Jackson (2000) and computationally by He et al. (2008) that this is a consequence of instability in the buffer layer beyond the dominantly viscous part of the boundary layer. More recently, He and Seddighi (2013) have associated such early responses to laminar-turbulent bypass transition. This effect is responsible for the different qualitative natures of the velocity histories shown in Fig. 2(a) at $y=R/16$ and $y=R/8$. In the core of the flow, the velocity increases linearly until the arrival of the radially inward propagating front caused by the instability initiated close to $y=R/16$. The inward propagation speed of the disturbance is shown by He and Jackson (2000) to be roughly proportional to the friction velocity of the initial flow and to be largely uninfluenced by the rate of acceleration. Close to the wall, the velocity behaves in a more laminar-like manner until the arrival of the outwardly propagating consequences of the instability.

The strong qualitative differences between the behaviors in successive stages of the overall flow excursion are especially pronounced in the acceleration history shown in Fig. 3(a). The general behavior is implicit in the velocity histories in Fig. 2(a), but the detailed behavior is clearer in Fig. 3. At very small times after the commencement of the overall acceleration, the local acceleration is almost identical at all $y>0$. At $y=R/16$, the rate begins to reduce significantly after about 0.1 s, but large changes are not seen in the graphs until much later. At about 1.2 s, the strong radially propagating front arrives at $y=R/8$ and causes a sudden, large reduction. It has similar, but less concentrated consequences at $y=R/4$ and $y=R/2$ at later times, and it reaches the pipe axis at about 2.0 s, where it “reflects” and begins to propagate radially outward. The amplitude of the reflection decays rapidly as it reaches successively greater radii.

Although the strong front has the greatest visual impact in Fig. 3, other features of the figure are also of interest. First, a small increase in acceleration is seen at $y=R/4$, $y=R/2$, and $y=R$ at approximately 1.3 s. Because the bulk acceleration of the flow is constant (prescribed condition), the increased acceleration in the core region will necessarily be accompanied by simultaneous decrease at other locations (as seen at $y=R/16$, for example). That is, on average, the inwardly progressing disturbance tends to increase the variation of acceleration over the cross section. Before the disturbance arrives at any particular radius, the flow at that radius and at all smaller radii (and hence greater $y/R$) accelerates at an almost common value. When the disturbance front arrives, the local acceleration initially reduces strongly and then evolves gradually toward a nearly steady value indicative of quasi-steady-velocity profiles. Eventually, the radial distribution of acceleration reflects the velocity profile shown at $R\approx \mathrm{120,000}$ in Fig. 1. This is much closer to a quasi-steady-flow profile than to the more uniform profile at $R\approx \mathrm{40,000}$.

The general behavior of the flow—in particular, the existence of the radially propagating (turbulent) viscosity front—is characteristic of suddenly accelerating flows. However, the fine details of the behavior will depend upon the exact nature of the acceleration. In the case considered above, the bulk acceleration is prescribed to be zero until $t=0$ and to be constant thereafter. As shown in Fig. 3(b), the amplitude of the required pressure gradient increases continuously. The rate of increase is initially small, but there is a rapid increase at about $t=1.3\mathrm{s}$ and an on-going increase thereafter. Because the prescribed bulk acceleration is constant, the whole of the variation in $\partial p/\partial x$ is attributable to changes in the wall shear stress. This is why the sudden rapid increase in the amplitude of $\partial p/\partial x$ at about $t=1.3\mathrm{s}$ coincides with the rapid increase in wall shear stress shown in Fig. 2(b).

The general behavior illustrated in Figs. 1–3 is consistent with experimental evidence presented, for example, by Maruyama et al. (1976), He and Jackson (2000), and He et al. (2011). The flow is characterized in three phases. Phase 1 is a period of delayed turbulence response, Phase 2 is a period of strong turbulent response, and Phase 3 is a continuing period of gradual response broadly similar to a succession of quasi-steady states. Phase 2 is found to be complex. The sudden response is measured first in the axial components of turbulence, and a further delay occurs before corresponding changes occur in the radial and tangential components. These three phases can be seen clearly in CFD predictions such as those shown by He et al. (2008) using RANS modeling and by Jung and Chung (2012) using large eddy simulation (LES) modeling. They can also be seen in DNS simulations presented by Seddighi et al. (2011) for flows in a planar 2D channel.

Throughout this paper, attention is focused on the particular example presented above. The details of the overall flow behavior are specific to this particular case, but the resulting conclusions are not. He et al. (2008) studied the behavior of the corresponding flows for a range of values of the initial Reynolds number, the pipe diameter, and the rate of acceleration. It was found that the essential behavior during Phase 1 was similar in all cases even though significant differences existed in Phases 2 and 3. Subsequently, He and Ariyaratne (2011) took an important step forward, showing that the Phase 1 results collapse to a common curve when the data are presented in a suitable nondimensional form. As a consequence, the use of a single case herein is sufficient provided that attention is focused on the behavior during Phase 1.

Fig. 4(a) compares wall shear stress histories predicted using different turbulence models, namely, (1) LS (as in Fig. 2), (2) Chang, Hsieh, and Chen turbulence model (CHC; Chang et al. 1995), (3) $\gamma -R_{\theta }$ (Langtrry and Menter 2009), and (4) $v^2–f$ (Durbin 1991). All four of these are so–called low–Reynolds number models of turbulence and so may be expected to characterize the accelerating flow much better than conventional equilibrium models that are tailored to steady flows. The 2D $\{r,x,t\}$ predictions are based on the commercial CFD package Fluent 13.0 using an axisymmetric computational domain with a pipe radius of 25 mm and a length of 8 m. Standard velocity inlet and mass outflow boundary conditions are applied at the inlet and outlet, respectively. The flow is effectively axially independent about 5 m from the inlet, and the conditions presented herein are even further downstream—at $x=7.5\mathrm{m}$ ($L/D=160$).

The governing equations are discretized using a second-order upwind spatial scheme and a second-order implicit temporal scheme described in detail by Gorji et al. (2014). The presented results are based on a time step of 0.0002 s and a spatial mesh of $100\times 30$ ($\mathrm{radial}\times \mathrm{axial}$). The normalized wall distance ($y^{+}$) of the wall adjacent node is less than 0.6 throughout the early stages of flow, thereby complying with the low–Reynolds number turbulence models criterion and avoiding the need for empirical wall functions. Systematic sensitivity tests have been performed to ensure near independence on the time step and the spatial mesh, especially during the early stages of flow—in which attention is paid below to small differences between solution pairs, not only to absolute values in individual solutions. Further simulations using time steps of 0.002 and 0.00002 s show negligible difference, as do further simulations using 140 and 180 radial mesh points, although a mesh of 70 points shows small differences. The axial mesh size is of low importance because all presented results are in the region of axially uniform flow.

By inspection, each of the turbulence models predicts similar qualitative behavior, but there are strong differences in the timing and character of the radially moving front triggered by the delayed turbulence response to the acceleration. The existence of such strong differences is clearly undesirable, but it is an inevitable consequence of the state of the art in the development of models of turbulence suitable for use in CFD-RANS simulations. This deficiency has been discussed, for example, by Scotti and Piomelli (2002) and Cotton (2007). More recently, Gorji et al. (2014) have compared predictions obtained using various turbulence models with DNS predictions for suddenly accelerating flow in a planar 2D channel, showing that each model tends to have advantages during some phases of the flow, but disadvantages in others. For example, some predict the timing of key events quite well, but are less good at predicting rates of change. Likewise, others predict aspects of the unsteady behavior quite well, but are less good at representing the initial steady-flow and the subsequent quasi-steadily accelerating flow. Mathur and He (2013) demonstrated that the implementation of the Launder-Sharma model of turbulence in Fluent 13.0 (and earlier versions) is unsatisfactory. They demonstrated much more satisfactory agreement using an alternative implementation. The improved LS model is used herein.

It is not the purpose of the present paper to attempt to infer reasons for the differences between the various different models of turbulence. Instead, the purpose is to assess the extent to which it is possible to develop models that are much simpler than those used herein, but that are nevertheless suitable for underpinning $1\mathrm{D}\{x,t\}$ simulations of unsteady flows. At first sight, this objective appears to be implausible in the light of Fig. 4(a), but that is not actually the case. On the contrary, advantage can be taken of a highly positive feature of the figure, namely, strong similarities in the predicted behavior during the early part of the acceleration before the sudden turbulence response. This is highlighted in Fig. 4(b), which shows a subset of the same data at an expanded scale. By inspection, there is very little sensitivity to the assumed turbulence model in this early period of the acceleration.

Strictly, the fact that all of the models yield similar predictions does not, in itself, guarantee that the predictions are meaningful. However, because a similar outcome was found by Gorji et al. (2014) for planar flows and included close agreement with DNS models, there is a high probability that the present predictions are reliable during the early period of the acceleration. It follows that use may reasonably be made of the outcomes in the development of approximate methodologies for use in $1\mathrm{D}\{x,t\}$ analyses. In particular, it is reasonable to use the CFD predictions as a basis for formulating methods of simulating evolving wall shear stresses in unsteady pipe flows—provided that the resulting models are intended for use only in periods of flow that are comparable with those considered in Fig. 4(b).

Although the ultimate aim of this paper is to assist analysts using $1\mathrm{D}\{x,t\}$ methods to deduce wall shear stresses such as those in Fig. 4(b), it is not expected that this will be done directly. Instead, it is acknowledged that methods of predicting steady-flow wall shear stresses are already widely available and are routinely built into commercial and other software packages. Instead of replacing these methods in their entirety, it is advantageous to retain them and to provide expressions that predict only additional contributions from the unsteadiness. That is, the target is to find adequate representations of an unsteady component of wall shear stress ${\tau }_{wu}$, defined as the difference between the actual stress ${\tau }_w$ and the corresponding stress in a steady flow with the same instantaneous flowrate ${\tau }_{ws}$, that is,

In a practical analysis, the quasi-steady and unsteady components ${\tau }_{ws}$ and ${\tau }_{wu}$ are evaluated independently and are then summed to give the required total ${\tau }_w$. Typically, the steady component is evaluated using expressions that have been developed and tested over many years. This two-step approach has the advantage of enabling 100% compatibility to be achieved when incorporating the analysis into software packages that already allow for quasi-steady shear stresses. Provided that the calculated unsteady component ${\tau }_{wu}$ tends to zero as the true influence of unsteadiness tends to zero, its use will not change predictions for steady flows. This is an important practical benefit because different steady flow approximations are used in various commercial packages—either explicitly or because of differences in the underlying numerical discretizations.

In principle, it is possible to develop expressions describing ${\tau }_w$ directly instead of describing only the unsteady component. In practice, however, this would downgrade accuracy because assumptions that are necessary to enable the development of analytical expressions that allow for unsteadiness would then also affect the accuracy of the quasi-steady component in end-user simulations. That would be a potential disadvantage in all flow regimes, and it would be especially unsatisfactory in regions of flow that are either steady or nearly so. Many such regions exist in typical analyses of unsteady flows in pipes. Fig. 5(a) shows variations of ${\tau }_{wu}$ predicted using each of the four turbulence models illustrated above. The values shown in the figure have been derived by undertaking independent simulations of steady flows with each model and subtracting the results from those shown in Fig. 4(a). However, instead of undertaking multiple independent simulations for a large number of truly steady flows, the quasi-steady component ${\tau }_{ws}$ has been obtained in the same manner as in Fig. 4 except for the use of an extremely small acceleration rate. To ensure convergence, the simulations were repeated until a further reduction by a factor of ten in the acceleration had negligible influence. This approach is highly convenient for the purpose of data manipulation, and it has the special advantage of ensuring 100% compatibility between the grid structures and the numerical algorithms used in determining ${\tau }_w$ and ${\tau }_{ws}$. This minimizes the magnification of rounding errors that are always possible when evaluating the difference between two values that are nearly equal—as is sometimes the case with ${\tau }_w$ and ${\tau }_{ws}$.

General unsteady flows commonly include periods of gently varying acceleration extending over large timescales corresponding to—or perhaps much greater than—those considered above. In such periods, the wall shear stresses tend not to vary greatly from values that apply in steady flows at the same instantaneous Reynolds numbers. That is, the unsteady component of shear stress ${\tau }_{wu}$ is only a small proportion of the overall stress ${\tau }_w$. In water-hammer-like flows, however, periods of gently varying flow such as this are punctuated by occasional changes that are both large and sudden. As a first approximation, each such change may be characterized as a rapid change in velocity that, if regarded as occurring linearly, is equivalent to an isolated, short-lived period of constant acceleration. Alternatively, the pulse may be regarded as a combination of an unending, constant, positive acceleration commencing at the beginning of the pulse and an unending, constant, negative acceleration (deceleration) of equal magnitude commencing at the end of the pulse. An example of this type of pulse is presented below. For present purposes, however, it is sufficient to note that the timescales associated with the passing of water-hammer wavefronts of this sort can be much smaller than those implied by Figs. 4(b) and 5(b). When this is so, the greatest influence of the unsteady wall shear stress can be modeled with reasonable accuracy without taking account of the delayed turbulence response in Phase 2 of the overall process.

Because the present analysis is developed for axially uniform conditions, it might appear illogical to apply its outcomes in studies of water-hammer-like flows. However, this is in fact reasonable. At any particular axial location, the change caused by any particular water-hammer wavefront approximates closely to a rapid change in axial velocity followed by a gradual adjustment of the radial velocity distribution. That is, the dynamics of the radial diffusion have little local influence on the wavefront itself (although the cumulative influence on a wavefront traveling large distances can be significant). Accordingly, for the purpose of modeling the radial changes, the influence of the waves at any location is closely equivalent to an applied bulk-acceleration history. For completeness, it is acknowledged that the existence of water-hammer wavefronts implies the existence of tiny density changes as well as velocity changes. However, these will have negligible influence in their own right on the processes that control the turbulent development. Likewise, no account is taken of other second-order phenomena such as radial acoustic modes in pressure wavefronts.

To develop an analytical model of ${\tau }_{wu}$ that can be incorporated into a $1\mathrm{D}\{x,t\}$ analysis, it is necessary to make use of approximations that represent the consequences of turbulence as functions of bulk flow properties such as the mean velocity, the mean acceleration, and rates of change thereof. Furthermore, it is important to avoid excessive increases in computational requirements in comparison with those used in conventional $1\mathrm{D}\{x,t\}$ analyses. The first of these requirements prevents the adoption of methods such as those used in CFD turbulence models, and the second prevents the use of models that would require extensive calculations at every grid point in every time step of an analysis.

The best current methods of meeting these objectives stem from work on unsteady laminar flow, for which fully analytical solutions exist for some forms of accelerating pipe flows. Constant viscosity examples include the well-known solution presented by Szymanski (1932) for acceleration from rest under the action of a prescribed, constant pressure gradient. Corresponding solutions for nonconstant pressure gradients include Uchida (1956) and Das and Arakeri (2000). Vardy and Brown (2010, 2011) have shown that such flows can be analyzed relatively simply by the use of Finite Hankel Transforms, and they have also used this approach to study flows in which viscosity varies with time. Other authors, that is, Costa and Macedonio (2003) and Adegbie and Alao (2007), have considered viscosity-dependence on temperature, albeit in ways that effectively reduce to a prescribed time-dependence. Ng (2004), Vasudevaiah and Rajagopal (2005), and Massoudi and Phuoc (2006) have obtained solutions for flows in which the viscosity exhibits pressure-dependence. In these cases, however, account needs to be taken of the pressure itself, not only the pressure gradient. As a consequence, the viscosity varies along the pipe, thereby introducing an additional complication that is avoided herein.

These examples use analytical methods, thereby yielding solutions that of general applicability. However, the penalty for achieving this is a necessary restriction to simple flow conditions. Zielke (1968) recognized that the applicability of analytical methods could be greatly extended by using a convolution to represent the wall shear stress as a function of the acceleration history of the flow. He showed how such a relationship could be obtained for flows with constant viscosity and then used independently in numerical simulations of general, unsteady laminar flows. Vardy et al. (1993) showed how an extension of Zielke’s method could be used in a similar manner for turbulent flows, making use of a highly simplified representation of turbulent viscosity used by Wood and Funk (1970). Subsequently, Vardy and Brown (1995, 2003) and other authors (e.g., Zarzycki and Kudzima 2004) have published successive improvements to the model based on the increasingly detailed representations of turbulent viscosity distributions. Further papers have addressed efficient numerical implementations of methods used to approximate underlying weighting functions.

The prescribed turbulent viscosity distributions used in all of these cases have two features in common. First, they approximate the true radial distributions of the effective viscosity, and second, they treat the chosen distributions as constant (frozen). These simplifications are necessary to achieve convenient analytical relationships between the wall shear stress and the bulk flow history, but they inevitably limit the range of approximate validity of the resulting expressions. The implications of the simplifications are explored in the following section.

Fig. 6(a) shows the evolution of the effective kinematic viscosity ($\mathrm{laminar}+\mathrm{turbulent}$) for the case presented above based on the LS model. By inspection, the influence of the wave-like propagation of the delayed turbulence response to the initiation of acceleration is even stronger than that of the acceleration and shear stress shown in Fig. 4, and this has big implications for the development of methods that rely on predetermined effective viscosity distributions. First, the actual distribution shown in Fig. 6 is highly specific to the particular flow conditions that give rise to it. Second, even if the outcome were less sensitive to the particular flow, it would be very difficult to incorporate it into a general method of analysis. Furthermore, the true situation is even more complex than this description implies because the detailed behavior of the turbulent viscosity is highly specific to the particular turbulence model used in the simulations. Indeed, very different detailed responses are obtained with other turbulence models—as illustrated by the surface plots shown in Fig. 7. These differences are a clear indication that turbulence models for RANS computations have not yet advanced sufficiently to address this type of flow reliably, and as a consequence, it must be concluded that it would be pointless to base the proposed $1\mathrm{D}\{x,t\}$ methods on the detailed outcomes of any particular 3D model. However, a much more positive conclusion is possible if attention focuses on conditions preceding the strong response. In this period, as shown in Fig. 6(b), the turbulent viscosity varies very little—although there is evidence of a tendency for relaminarization close to the wall. Furthermore, the predicted values of the effective viscosity vary even more slowly than those of the shear stress [Fig. 2(b)]. This may be interpreted as strong evidence in support of the so-called frozen-viscosity models of unsteady skin friction developed, for example, by Vardy and Brown (2003).

Even during the period when frozen-viscosity assumptions are plausible, the proposed $1\mathrm{D}\{x,t\}$ models of unsteady wall shear stress cannot be based on detailed radial variations of turbulent viscosity because these are inherently unknown during $1\mathrm{D}\{x,t\}$ analyses. Instead, it is necessary to use approximate distributions of the turbulent viscosity, and so it is important to know which characteristics of the true distributions need to be reproduced faithfully and which may be modeled in a simplified manner without significant loss of accuracy. This matter is addressed in Fig. 8, which shows predicted variations of the unsteady component of wall shear stress based on various approximations to the true distribution of effective viscosity. Also shown, for comparison purposes, is the wall shear stress in the preexisting steady flow. The curves labeled ND Solve (numerical solution software; NDS) are $1\mathrm{D}\{r,t\}$ numerical solutions of the axisymmetric Navier-Stokes equations obtained using the software NDSolve in the commercial package Mathematica. In each case, the initial velocity distribution $u\{r,0\}$ is copied directly from the baseline CFD solution, as is the pressure-gradient history $dp/dx\{t\}$ throughout the simulation period. In the particular case of Fig. 8(a), the viscosity distribution throughout the simulation period is also copied, but this case is only for reference purposes. In the remaining cases, the viscosity distribution is chosen independently and is held constant for the duration of the simulation.

Extensive checks have been undertaken to ensure the validity of the resulting solutions. The most important of these are

In Fig. 8(a), the prescribed distribution of the effective viscosity is unmodified from the CFD output, and so the predicted shear stress should be identical to that obtained from the CFD analysis. In practice, the agreement is so close that it is difficult to distinguish between the two graphs in the figure. This provides validation that is not usually practicable for CFD analyses, and in particular, it confirms the suitability of the numerical grid structures used in both analyses. The original reasons for using different methodologies stemmed from different backgrounds of the various authors of the paper, but the cross-validation of their predictions considerably increases confidence in the robustness of the predictions.

Fig. 8(b) assesses the validity of using frozen-viscosity models of turbulence in this type of flow. That is, the prescribed viscosity distribution throughout the whole simulation is identical to that $t=0$ in Fig. 8(a). Once again, the agreement is very close, thus demonstrating that, during the period simulated, the evolving shear stress depends much more strongly on the initial turbulence state than on subsequent deviations from it. This validates the frozen-viscosity principle. That is, the primary response of the flow to the suddenly imposed change of pressure gradient occurs before turbulence diffusion can exert its influence. This is consistent with the evidence presented in Fig. 6.

Fig. 8(c) begins to assess the detail with which it is necessary to model the radial distribution of the effective viscosity. Instead of using the distribution obtained from the CFD analysis, the distribution input to NDSolve is a simple constant, namely, the molecular viscosity itself. That is, the evolving flow is analyzed as though it were laminar—although, of course, the prescribed initial velocity distribution would not be possible in a steady laminar flow. As a consequence, the flow state would change after $t=0$ even if there were no change in the imposed pressure gradient. Therefore, in the enforced accelerating flows studied herein, this is a second cause of change after $t=0$. However, except for very small imposed accelerations, the inertial cause is much stronger than the nonequilibrium state of the initial velocity profile, which can exert its influence only through diffusive means.

Once again, the predicted wall shear stress is very similar to that for the fully turbulent case. At first sight, this outcome might appear to merit surprise. Indeed, the authors did not initially expect such close agreement. However, with the benefit of hindsight after numerous simulations with alternative assumed viscosity distributions, the explanation has been found to be simple, namely, that the early-time response of the flow depends almost exclusively on the assumed distribution close to the wall. Remote from the wall, the whole flow accelerates at the same rate with negligible change of turbulence conditions or gradients of local mean-velocity, and so on. Very close to the wall, velocity gradients increase continually in time, but even here, the viscosity distribution remains almost frozen for a finite period. Thereafter, as predicted by He et al. (2008) using CFD-RANS and confirmed by Seddighi et al. (unpublished data, 2013) using DNS, local instability begins near the outer edge of the buffer layer and then diffuses both toward the wall and toward the axis. Until this diffusive process penetrates deep into the wall layer, diffusion therein is dependent on the preexisting viscosity distribution. For smooth-walled flows, this implies that it is governed by laminar viscosity alone.

In addition to showing why the early-time behavior near the wall in a strongly accelerating flow is influenced so strongly by the near-wall viscosity distribution and so weakly by the more remote viscosity distribution, this explanation shows why the assumed nonequilibrium velocity distribution has negligible small-time influence. It can exert influence only through the same diffusive process that is such a weak partner in the early response to the imposed overall acceleration.

This deduction has unfortunate implications for the first authors, who have developed convolution models of unsteady wall shear stress based on a bilinear approximation to the distribution of effective viscosity (see next section). As shown in Fig. 9, the bilinear shape represents the overall distribution much more realistically than a simple assumption of laminar flow, but it does not represent the wall region well. Such a distribution is used in the simulation presented in Fig. 8(d), and it strongly overestimates the unsteady component of wall shear. The agreement with the CFD predictions is substantially poorer than that based on other approximations, including even the wholly laminar distribution. This demonstrates an important limitation of the bilinear approach, but it also points the way forward to eliminating the deficiency, namely, by improving its representation of the true viscosity distribution close to the wall as is done in models developed by Zarzycki and Kudzima (2004), for instance.

It is useful to outline the process followed in the derivation of convolution models of the unsteady component of wall shear stress for use in $1\mathrm{D}\{x,t\}$ simulations. First, $2\mathrm{D}\{r,x,t\}$ analytical solutions are obtained in the Laplace domain based on the prescribed frozen distribution of viscosity, with the pressure gradient treated as a parameter. Independent solutions are obtained for unsteady flow and steady flow, and the difference is used to express the unsteady component (${\tau }_{wu}$) as a function of the pressure gradient. When desired, the Laplace inverse of this solution can be used to provide a relationship between ${\tau }_{wu}$ and the pressure gradient history. In practice, however, it is more useful to undertake additional development in the Laplace domain to enable the final result in the time-domain to be expressed as a relationship between ${\tau }_{wu}$ and the bulk-flow acceleration history. This formulation is beneficial in numerical simulations of unsteady flows because the numerical evaluation of spatial derivatives such as $\partial p/\partial x$ can be unreliable close to boundaries whereas temporal derivatives such as $\partial U/\partial t$ have no such limitations (provided that the boundaries are stationary).

Attention needs to be drawn to an important feature of this process, namely, that the same viscosity distribution is used in the independent developments for steady and unsteady flow. It is usually, but not necessarily, taken to be equal to that in the initial steady flow. This is an obvious choice for end-user applications in which the initial flow is steady. As shown above, the assumed frozen state is a good approximation for the fully unsteady flow, but it is a poor approximation for the associated steady flow analysis because the true viscosity distribution in a fully developed quasi-steady flow is Reynolds number–dependent, and therefore evolves as the acceleration proceeds. It is not possible to reflect the true behavior in the convolution models discussed above because the analytical methodology requires the use of identical viscosity distributions in the independent steady and unsteady flow analyses. Without this constraint, the method used to infer an unsteady component from differences between the steady and fully unsteady solutions would be invalid. Furthermore, it is not feasible to allow for simultaneous variations of viscosity in space and time in the Laplace domain analysis. Important consequences of this limitation of convolution models of unsteady friction are illustrated in the following section.

The need to use identical viscosity distributions in the steady and unsteady parts of the development has big implications for the interpretation of Fig. 8 above. In that figure, the CFD results can be assumed to characterize the true behavior. In particular

Likewise, the results presented for the prescribed-viscosity distribution simulations satisfy

As a consequence of these choices, differences between the CFD results and the corresponding results using prescribed-viscosity distributions may be used to infer the influence of the assumed effective viscosity distributions on the response of the wall shear stress to sudden accelerations. The figure should not, however, be regarded as an exact indicator of the performance of actual convolution models of unsteady friction based on prescribed viscosity distributions. This is because the analytical development of these is necessarily simplified in comparison with the numerical approaches used above. It may be described as

The most important consequence of using a unique, predetermined viscosity distribution in the steady flow simulations is that, in contrast with results in Fig. 8, the implied steady flow shear stress is necessarily incorrect and becomes increasing so as the acceleration proceeds. Furthermore, although the error may be only a small proportion of ${\tau }_{ws}$, it rapidly becomes a large proportion of ${\tau }_{wu}$, which is the sole output of the convolution models of unsteady friction. This greatly limits the range of applicability of the methods—although, fortunately, its importance is smallest at early times, and this is when ${\tau }_{wu}$ is greatest.

The importance of the method of choosing initial conditions is less easy to quantify. Clearly, initial conditions based on approximate viscosity distributions will be incorrect, and the errors will be greatest for the least satisfactory approximations. However, identical initial conditions are used for both analyses (steady and unsteady), so the effect is cancelled at $t=0$. The influence at greater times will depend upon the particulars of the flow history in the end-user application, but two general issues can be identified. To illustrate these, note that $\tau =\mu \partial u/\partial y$ and compare the consequences of this for (1) the bilinear approximation, and (2) the wholly laminar approximation. For the bilinear case, the dimensions of the distribution illustrated in Fig. 9 are (implicitly) chosen by the end-user to achieve the correct value of ${\tau }_{ws}$ for the chosen Reynolds number. However, as shown above, the assumed nonuniform viscosity distribution near the wall leads to an overestimation of the consequences of subsequent accelerations. In contrast, for a wholly laminar viscosity, it is not possible to achieve compatibility with the end-user’s initial conditions. If the wall shear stress (and hence the pressure gradient) is matched, the implied initial Reynolds number in the original analytical model will be far too great. Alternatively, if the Reynolds number is matched, the implied initial wall shear stress and pressure gradient will be far too small. Once again, however, these statements are true for both parts of the weighting function development, so the effect cancels at $t=0$. If it also cancels (or nearly so) as time increases, the consequences of the initial mismatch might be small. In that case, it would be reasonable to expect the performance of convolution models based on a wholly laminar approximation to outperform those based on a bilinear approximation, even though the individual components of the latter are a closer match to the true end-user condition.

Given this uncertainty, it must be concluded that there is no clear evidence to justify the use of convolution models based on the bilinear approximation in preference to that based on a wholly laminar viscosity distribution. Models such as those developed by Zarzycki and Kudzima (2004) allowing for the near-wall behavior shown in Fig. 9 should outperform the bilinear model, and they might also outperform the wholly laminar approach, especially at larger times.

In addition to assessing the validity of the underlying bases of convolution models, it is worth drawing attention briefly to inaccuracies that commonly arise when they are implemented in end-user software. Usually, to achieve a big reduction in the required CPU time, the weighting function is approximated by a series of exponentials, broadly in a manner originally proposed by Trikha (1975). However, it is important to ensure that the approximation is sufficiently accurate over the whole range of frequencies to be studied by the software. Furthermore, irrespective of whether exponential approximations are used, it is important, but not straightforward, to ensure that the influence of the most recent time step in the calculation sequence is evaluated accurately. This difficulty arises because the weighting function is infinite at its time origin. These issues have been addressed by, for example, Kagawa et al. (1983), Schohl (1993), and Vardy and Brown (2007).

Fig. 10 shows three numerical solutions obtained using a particular end-user software package for $1\mathrm{D}\{x,t\}$ unsteady flow simulations, namely, ThermoTun (www.ThermoTun.com) in which the bilinear model of wall shear stress has been implemented (see the following section). The package is usually used for simulating pressure waves in air, but it can also be used to simulate unsteady liquid flows. In these examples, the fluid and pipe properties are the same as those used above (i.e., $\rho =\mathrm{1,000}{\mathrm{kg}/\mathrm{m}}^3$, $\nu ={10}^{-6}{\mathrm{m}}^2/\mathrm{s}$, $D=50\mathrm{mm}$). In Fig. 10(a), a constant fluid acceleration is prescribed and is identical to that used above (i.e., $\partial U/\partial t=0.9{\mathrm{m}/\mathrm{s}}^2$). By inspection, the evolution of the predicted unsteady component of wall shear stress in the early stages of the acceleration approximates fairly closely to that in the simulations presented above. However, whereas the CFD values reach a maximum and then reduce and become negative, the $1\mathrm{D}\{x,t\}$ prediction increases monotonically and, at larger times than those shown in the figure, becomes asymptotic to a constant value. This qualitatively different behavior is an inevitable consequence of the limitation highlighted in the previous section, namely, the use of frozen-viscosity assumptions for the quasi-steady component of flow as well as for the overall flow. This limits the duration of the period of approximate validity of the solution.

In the curve labeled Ramp Accel in Fig. 10(b), the initial stages of the flow are identical to those in Fig. 10(a), but the acceleration ceases abruptly, thereby creating finite acceleration pulse and hence a finite ramp in the velocity history. The duration of the ramp $\mathrm{\Delta }{t}_{\text{pulse}}$ has been chosen as 0.1 s, so that it is within the period for which reasonable agreement is shown in Fig. 10(a) for the bilinear behavior under consideration. Immediately after the acceleration ceases, the predicted value of ${\tau }_{wu}$ begins to reduce rapidly in a mirror-image manner to the original increase. At any instant during this period, the imposed flow may be regarded as the sum of (1) the original uniform acceleration, and (2) a second uniform acceleration of equal amplitude but opposite sign, commencing 0.1 s after the first. As a consequence, the net value of ${\tau }_{wu}$ can be inferred from the uniform acceleration case. For example, the value shown for the RampAccel curve at the instant $t=0.4\mathrm{s}$ is the difference between the values shown for the UniformAccel curve at 0.4 s and at 0.3 s. This is the basis of the method used in end-user software to implement the method of predicting the unsteady component of the wall shear stress ${\tau }_{wu}$. The instantaneous velocity is regarded as a sum of successive finite acceleration pulses, each of which can be treated as in the preceding description, and the instantaneous value of ${\tau }_{wu}$ is the sum of the individual contributions. In general, each acceleration pulse will be of a different amplitude from the others, but that is not a major complication because the data shown in Fig. 10 can be expressed in a universal form (for any particular pipe) by scaling the vertical axis by the acceleration. In fact, the generalization can be extended to allow for multiple pipes, but that is of little importance for the purposes of this paper.

This process involves important approximations that significantly limit its range of validity. As indicated above, the accuracy of the predicted contribution of each pulse depends upon the suitability of the assumed frozen viscosity distribution and upon the consequences of using this distribution for the evaluation of the steady-flow component ${\tau }_{ws}$. In end-user software, these constraints are even more significant than is obvious from the above illustration because, in principle, different frozen-viscosity assumptions should be made for each acceleration pulse in the summation of pulses representing the historical changes in velocity. It is therefore fortuitous that, for smooth-walled flows that are the main focus herein, the frozen-viscosity can be represented with sufficient accuracy by uniform laminar flow.

Fig. 10(b) includes a curve labeled StepChange. For this case, the assumed acceleration is ten times greater than that used for the RampAccel curve, but the duration of the pulse is reduced by a factor of ten so that the overall velocity step is the same in both cases. By inspection, this causes a much larger value of ${\tau }_{wu}$ at small times, but the values at larger times are comparable with those for the slower acceleration. The duration of the pulse is so short that it cannot be seen easily in the figure, but it nevertheless has a major influence on the amplitude of ${\tau }_{wu}$ at very small times. Indeed, if the step were instantaneous, the predicted value of ${\tau }_{wu}$ at $t=0$ would be infinite. This is a further source of complications in end-user software—although, with suitable care, it can be overcome in a relatively straightforward manner because the numerical integrations are undertaken over finite time intervals, and the predicted ${\tau }_{wu}$ is finite at all instants $t>0$.

All of the examples presented above are based on smooth-wall flows. This is an important class of flows, and it has the big advantage of enabling accurate representations of conditions very close to wall surfaces. However, a large number of practical applications of unsteady flow involve rough-wall flows—either fully or partially rough. In these flows, the influence of turbulence extends as far as the wall, and so the use of an effective viscosity distribution with a laminar wall region is not appropriate. Instead, the effective viscosity at the notional wall surface exceeds the laminar viscosity, and its value increases with increasing distance from the wall. This condition is more similar to the bilinear approximation shown in Fig. 8 than to the smooth-wall approximations based on a viscous layer adjacent to the wall. Thus, although it is suggested above that Zielke’s laminar flow expressions should be used in preference to the Vardy-Brown smooth-wall relationships, this warning does not apply to the corresponding rough-wall relationships (Vardy and Brown 2004, 2007).

Irrespective of whether flows approximate to smooth-wall or rough-wall conditions, the importance of unsteady components of wall shear stresses tends to reduce with increasing Reynolds number. This is because the quasi-steady component tends to increase approximately with the square of the velocity whereas the unsteady component increases approximately linearly with the acceleration. Duan et al. (2012) give further details, with special reference to the implications for the overall damping of water-hammer flows. A related issue is the practicality of using small-scale models (physical or numerical) to infer behavior at large scale (or vice versa). When rapid accelerations exist, as in pressure waves such as those in water-hammer applications, damping caused by unsteady friction is unlikely to scale in a Reynolds-number independent manner.

Three characteristically different phases of turbulence response to a suddenly imposed constant acceleration of an initially steady, smooth-wall pipe flow have been described with particular reference to implications for existing methods of modeling unsteady wall shear stress in $1\mathrm{D}\{x,t\}$ software using frozen-viscosity hypotheses. The main conclusions from the work may be summarized as follows:

The authors gratefully acknowledge funding from the UK Engineering and Physical Sciences Research Council (EPSRC) through grants EP/G068925/1 and EP/G069441/1 for some of the work reported in this paper.