Estimating turbulent kinetic energy and dissipation with internal ﬂow loss coeﬃcients

,


Introduction and background
Turbulent kinetic energy in internal flow systems can often only be found through computational fluid dynamics (CFD), because measuring turbulence quantities in internal flows is usually difficult.Reynolds averaged CFD of single components is presently tractable, but computationally expensive.Accurate CFD of complex piping systems in practice is difficult, if not prohibitively expensive.Further, one source of error in CFD is the specification of inlet turbulent kinetic energy and dissipation boundary conditions.Given these concerns, there is a need for a computationally tractable methodology to predict turbulent kinetic energy in internal flows.
There's also an interest in extending the Bernoulli equation to more general contexts.There are two approaches to generalization.The first approach is to find paths other than streamlines where the Bernoulli integral is constant under more general conditions, e.g., in viscous flows.Brutyan et al. [2] identify such paths for general steady viscous flows.The second approach is to determine how the Bernoulli integral varies along streamlines under more general conditions.Typically, this is only done for internal flows.Grose [3] and Synolakis and Badeer [4] develop viscous corrections to the Bernoulli equation.A more common approach adds a term for energy "lost" along the flow path in internal flows.The equation has been called the "mechanical energy balance", the "engineering Bernoulli equation", or the "extended Bernoulli equation".The models for the lost energy are ∆p = ζ • 1  2 U 2 for pipe fittings ("minor loss") and valves and ∆p = ( f L/d) • 1 2 U 2 for pipe segments ("major loss").The lost energy is attributed to viscous dissipation by Panton [5, p. 133] and Bird et al. [6, p. 204].However, the previous derivations of the loss term are only valid for steady flows, precluding the possibility of anything one would call turbulence.Pope [7, p. 125] notes that the (laminar) mean flow dissipation term the previously studies highlight is generally negligible.Consequently, how the lost energy is distributed in a turbulent flow is unclear.A relationship between the loss and turbulence quantities like the turbulent kinetic energy and dissipation exists, and it will be detailed in this paper.Past researchers have attempted to relate flow losses to turbulence quantities, without clear success.Most previous theories modeled turbulence reduction by screens as a function of the loss coefficient of the screen.Loehrke and Nagib [8, p. 5] suggested that the theories disagree with the data because the theories treat the screen as only a turbulence suppression device, and do not include any way for turbulence to be generated.Examining both screens and honeycombs, Scheiman and Brooks [9] agree with that conclusion.Groth and Johansson [10] made the same suggestion and showed experimentally that the turbulence intensities immediately downstream of the screen are higher than that upstream of the screen, but the turbulence decays further downstream of the screen.Groth and Johansson's measurements show that both turbulence generation and dissipation are factors.
For modeling screen turbulence, Baines and Petersen [11, p. 471R] proposed an equation similar the Bernoulli equation along a streamline with additional terms for turbulent kinetic energy and dissipation.No derivation was provided.The equation was justified by stating it represents energy conservation, however, this is implausible as 1 there are no terms for energy transfer between streamlines.No link between this equation and flow losses was made.Focusing on more general hydraulic resistances, Nikitin and Nikitina [12] developed models for the maximum RMS pressure p using the friction factor or loss coefficient in pipe and general hydraulic resistances.Nikitin and Nikitina used dimensional analysis, not energy conservation, to determine p , and the equations developed are valid only if the inlets have negligible turbulence.Independent from the attempts to model turbulence through screens or general hydraulic resistances, researchers interested in liquid jet breakup developed simplified models to estimate turbulent kinetic energy at the outlet of a nozzle.Early models were developed by Natanzon [13], Tsyapko [14], and Jackson [15, p. 111].More notable, however, is the model of Huh et al. [1,16] (abbreviated as "Huh's model" here).In Huh's model, the turbulent kinetic energy at the nozzle exit can be estimated with the equation where k 0 is the plane averaged turbulent kinetic energy at the nozzle exit, U 0 is the (plane averaged) velocity through the nozzle orifice, L 0 is the nozzle orifice length, d 0 is the nozzle orifice diameter, C d is the discharge coefficient of the nozzle, ζ c is the contraction loss coefficient, and c ≡ A in /A 0 = (d in /d 0 ) 2 is the area contraction ratio of the nozzle.See figure 1 for an illustration of a conical nozzle with a cylindrical orifice using the previous notation.This model has two major problems evident by inspection.First, in this model the nozzle exit turbulent kinetic energy is not a function of the nozzle inlet turbulent kinetic energy.This is inconsistent with the evidence that strong turbulent kinetic energy at the nozzle inlet can affect the stability of a liquid jet [17,18].Second, this equation implies that turbulent kinetic energy grows arbitrarily large as L 0 /d 0 decreases to zero, and that turbulent kinetic energy goes to zero as L 0 /d 0 grows arbitrarily large.Both are false.Klein [19,p. 246,fig. 3] reviews past measurements of turbulent kinetic energy as a function of development length.For small L 0 /d 0 , k approaches the inlet value, and as L 0 /d 0 increases to infinity (i.e., the flow becomes fully developed), k goes to a value determined by the friction factor for fully developed flow.The latter value is independent of the inlet turbulent kinetic energy.Further, if k grew arbitrarily large as L 0 /d 0 decreased, this would violate energy conservation.k increasing would also imply that the stability of a liquid jet would be worse for shorter nozzles, but typically the opposite is true [20].
These problems are caused by an error made in the derivation of equation 1. Huh et al. use a force balance in an attempt to relate the turbulent kinetic energy to the pressure drop across the nozzle orifice.Presumably the "turbulent stress" used in the model is the Reynolds shear stress ( uv if one divides by the liquid density), but this is not the same as the turbulent kinetic energy (k), what is desired in the model.The model implicitly assumes the two are the same or at least proportional.One can relate uv to the turbulent kinetic energy via the stress-intensity ratio, | uv /k |, found to be about 0.3 in many shear flows [21, pp. 116, 121, 126, 138].However, this option is not viable, as using wall friction requires taking the force balance at the orifice walls.Both uv and k are zero at smooth walls due to the no-slip and no-penetration conditions [7, p. 269].The stress at a smooth wall comes entirely from the viscous component * .Finally, the constant stress-intensity ratio approximation is not particularly accurate [22].
Further, k is determined solely by the nozzle orifice walls in Huh's model, explaining why inlet turbulent kinetic energy does not factor into the model.In addition to ignoring the effect of inlet k, this implies that the contraction does not change the turbulent kinetic energy.The effects of contractions on turbulent kinetic energy and anisotropy are well known [23][24][25].Rapid distortion theory (RDT) can model these effects.
The failures of previous models prompted the development of a more rigorous approach which also uses empirical loss coefficients.

Derivation of the turbulent Bernoulli equation and relationship between loss and turbulence quantities
We start with a standard derivation of the Bernoulli theorem, albeit from the Reynolds averaged equations rather than the instantaneous equations.We depart from the standard derivation when we take the new term for the work done by the Reynolds stress, and decompose that term into the production of turbulent kinetic energy and a turbulent flux.This returns an equation including terms for energy transfer between streamlines.Algebraic manipulation and integration over the volume of the internal flow path returns a relationship between flow losses, k, and ε.For brevity, the ensemble mean velocity in the Reynolds decomposition will be denoted with a capital U i rather than U i , and the fluctuating terms will be denoted with lowercase u i .* That is not true for the special case of rough walls, but in this case one would still need to know the viscous component of the wall stress, which is not generally known.The link between the plane average value of uv (to find k) and the value of uv at the wall is also unclear.
There is no simple link between τ w and k even for rough pipes.In the remainder of this paper, we assume that all fluctuations are zero at the wall.

Differential turbulent Bernoulli equation for a streamline
Start with the Reynolds averaged Navier-Stokes equation in Stokes form and consider the statistically stationary case (no change in the mean with time) [7, p. 23, p. 86]: Rewrite the previous equation so that only i j k U j ω k is on one side, and the remaining terms are on the other: Consider a vector s i which is in the direction of a streamline at a point in space.Multiply the previous equation by a differential element of this vector: The cross product term has no component in the direction of U i , which is the same as s i , so i j k U j ω k ds i is zero.Note that the local differential length element along a streamline can be related to the local time-averaged velocity vector by ds i = U i dt.This is done for the Reynolds stress term only: Now, note that the following decomposition is a consequence of the product rule: The second term is recognized as the production of turbulent kinetic energy.Substituting the expression for turbulent kinetic energy [21, p. 49] into the result above returns The previous equation can be used to see how turbulent kinetic energy changes along a streamline.Aside from the first and third terms, everything else is between streamlines.Later, when we average over the projected area of the pipe, the other terms go to zero at the walls as there is no transfer between the pipe and the walls due to the no-slip and no-penetration boundary conditions.These terms may need to be kept for inlets and outlets, however.
Substituting the result from equation 7 into equation 5 returns (after application of U i = ds i / dt to the turbulent kinetic energy derivative and some rearrangement)

Turbulent Bernoulli equation for a streamtube with no-slip and no-penetration boundary conditions
Equation 8 is valid along a particular streamline, but there are many terms which need to be modeled, limiting its utility.Consequently, we seek to develop a version for a streamtube applicable in internal flows where there are no-slip and no-penetration boundary conditions, like the "extended Bernoulli equation".Multiply all terms by ρ U l n l dA.n l is a unit vector which is normal to the integration area dA.The overall area integrated over is, for example, the cross sectional area of a pipe or pipe fitting.
For the first term, note that ∂ ∂ x i (• • • ) ds i is the change in the quantity inside of the parentheses along a streamline, so it can be written as d (• • • ).For the second term, note that ds i = t i ds, where t i is a unit vector tangent to the streamline.For the third term, note that dt = ds/ U m U m , where s is the arclength of the streamline.These changes result in The previous result is simplified by noting that U l / U m U m = t i : Now, we make the approximation t l n l = 1, which means that the streamlines are normal to the surface on each end of the streamtube.This is not satisfied in internal flows aside from long straight stretches of pipe, but it is approximately satisfied.A consequence of this approximation is that U i n i = U i U i ≡ U.After this, we have The viscous mean flow term can be decomposed as below, which follows from the product rule: The second term is the dissipation due to the mean flow, ε m , as defined by Pope [7, p. 124].Applying this decomposition to equation 12 and simplifying the result returns Now integrate equation 14 over the volume of the streamtube, note that ds dA ≡ d -V , and use the kinetic energy coefficient α [26, p. 218] (α ≥ 1, with α = 1 for only a uniform velocity profile if all velocities are forward): The averaging in the first term is averaging in the area normal to the streamlines, e.g., in a pipe flow under the approximation previously mentioned, this is the plane of the pipe.All of the terms in the second term aside from the dissipation terms go to zero at the walls after the divergence theorem, and applying the no-slip and no-penetration boundary conditions (U i and u i are zero at the boundaries).(If it is unclear, ∂ k ∂n = 0 at the wall due to k's quadratic dependence on the velocity, where n is the normal direction to the wall.)Unfortunately, the divergence theorem does not help at the inlet and outlet of the streamtube, as the terms of interest are not necessarily zero.However, as a hypothesis they may be zero or otherwise small compared against the others, or possibly well modeled similarly to the dissipation.The validity of this approximation may depend on where the control volume ends.Under this approximation the turbulent Bernoulli equation is The deviation from the typical Bernoulli equation is generally called "loss": The energy "loss" is decomposed into three components: 1. ∆k, which is mean flow energy converted into turbulent kinetic energy, 2. ( ρ/ ṁ) - V ε m d -V , which is energy dissipated by the mean flow, and 3. ( ρ/ ṁ) - V ε d -V , which is energy dissipated by turbulence.To reiterate a point made earlier, Panton [5, p. 133] and Bird et al. [6, p. 204] only include the term with ε m , which Pope [7, p. 125] notes is generally negligible in turbulent flows.

Differential loss decomposition for pipe flows
A differential version of the energy loss decomposition for pipe flows can be useful.Recognize that in a pipe flow ζ ≡ x 0 ( f (s)/d) ds.Consequently, the loss can be written as (integrating from point 0 to point x) Differentiating this expression with respect to x returns d dx Turbulence modeling Equation 17 is simple, but even with a known loss coefficient, k and ε can not be uniquely determined.A turbulence model is needed to estimate quantities of interest.Two simple models are detailed.To reiterate, we are are neglecting certain terms for simplicity.Future works should relax these approximations.

Dissipation fraction model
When modeling individual piping components, e.g., valves and fittings, the easiest model would simply assume a certain fraction of the loss energy is dissipated.Presumably typical values of this "dissipation fraction" could be tabulated alongside loss coefficients if it proves to be relatively universal for classes of valve and fitting geometries.
Following equation 17, we define the dissipation fraction as α ε can be assumed constant for a pipe system component.Note that α ε is only bounded below by zero.There is no obvious upper bound to the amount of dissipation.In some circumstances, e.g., where turbulence reduction is desired, the amount of dissipation will exceed ζ • 1 2 U2 , causing ∆k to be negative and α ε to be greater than one.

Dissipation scaling model
A common scaling used to estimate turbulence dissipation is ε = C ε k 3/2 /Λ, where Λ is the integral length scale [7, p. 244].Assuming that the integral length scale is proportional to a characteristic diameter, for example, the pipe diameter, the volume averaged dissipation can be modeled as This leads to the following equation for ∆k if k in the dissipation model is taken as the inlet value, k 1 , which simplifies the computation: for a single flow resistance between 1 and 2, where -V c is the characteristic volume of the component.Unlike the dissipation fraction approach, it is impossible to solve for the overall increase in the turbulent kinetic energy with a simple sum.The model requires solving for k at each node in a pipe system.

Models for specific internal flow situations Fully developed pipe flow
Equation 19 can be applied to the case of fully developed pipe flow.Here dk dx = 0 and f is also a constant, independent of the location.The dissipation scaling model leads to the equation † , which implies Tu FD ≡ 2k FD /3 To test this theory, an experimental database for both rough and smooth pipes was compiled, selecting only studies which measured all three velocity components [27][28][29][30][31][32][33] (17 points, 9 smooth, 8 rough).Fitting the theory with least squares returns Tu FD = 0.2458 f 1/3 (R 2 = 0.9089).Fitting a general power law returns Tu FD = 0.3655 f 0.4587 (R 2 = 0.9753).Neglecting measurement error, the exponent with 95% error is 0.4587 ± 0.0401.The scaling Tu ∝ f 1/3 is not consistent with this.A power closer to 1/2 seems justified.The scaling Tu ∝ f 1/2 would follow from the assumption that k 1/2 ∝ u τ .This discrepancy is likely due to the neglect of many terms when constructing the model.It is interesting to note that the power is between the two theories, suggesting reality blends both physics.
One further observation from the correlation is that the Blasius friction factor law for smooth pipes, f = 0.316 Re −1/4 , suggests that Tu FD decreases as Re ≡ U 0 d/ν increases, contrary to many expectations.This trend is consistent with experimental measurements at the centerline [28, p. 35, fig. 15].

Contracting nozzle
Contractions followed by cylindrical segments (e.g., figure 1) are the typical nozzle geometry in turbulent jet breakup.To offer an alternative to Huh's inaccurate model, we will apply RDT to the contraction and equation 19 to the cylindrical segment.We can then compute the overall turbulence intensity at the end of the contraction:

Estimation of inlet turbulent kinetic energy
Models of this variety typically require the inlet turbulent kinetic energy.There are a few possibilities to estimate this quantity: 1. Find the inlet turbulent kinetic energy from empirical measurements; 2. Select the inlet turbulent kinetic energy to match other measured quantities (less desirable); 3. Assume the flow is fully developed at the inlet; 4. Assume the inlet turbulent kinetic energy is zero, which may be acceptable for laminar or essentially quiescent entrances; or 5. Use standardized tabulated empirical values of turbulent kinetic energy.At present, these tables do not exist, but it seems plausible that one could measure and compile values of turbulent kinetic energy at the outlet for pumps and other common starting points.

Conclusions
While Huh's model is inaccurate, relating k to the pressure drop in the nozzle is valid and has a long history.Unfortunately, no simple trick like a force balance at the orifice walls can allow one to avoid turbulence modeling.Using standard turbulence modeling approaches we developed equation 24, which can be used in place of Huh's model.However, equation 24 has not been validated.We recommend instead using CFD to determine k for now.Earlier in this paper, the side boundaries of the streamtube were assumed to be no-slip.This assumption is false when one or more recirculation zones appear in the flow.Briefly, we'll show that recirculation zones bounded outside by no-slip boundaries have the same effect as no-slip boundaries.In figure 2, the solid line is a no-slip boundary, the dashed line is the boundary (dividing streamline) of the central streamtube (denoted with c), the dotted line from the left corner (point 1) to the large dot is the inlet to the top of the recirculation zone (denoted with rt; also the outlet of the bottom of the recirculation zone, denoted with rb), and the dotted line from the large dot to the reattachment point (point 2) is the rt outlet (rb inlet).Equation 15 can be decomposed into 0 = ∆ c + in/out,c + slip,c + no-slip,c for c, 0 = ∆ rt + in/out,rt + slip,rt + no-slip,rt for rt, and 0 = ∆ rb + in/out,rb + slip,rb + no-slip,rb for rb.∆ are the conservative terms.The "in/out" integrals are the inlet and outlet integrals.The slip and no-slip integrals are the integrals over surfaces with slip and no-slip boundaries, respectively.The terms no-slip,rt and slip,rb equal zero by construction.Because rt and rb form a loop, ∆ rt + ∆ rb = 0.And because the inlet of rt equals the outlet of rb, and vice versa, in/out,rt = − in/out,rb .The central streamtube can be connected to rt by noting that slip,c = − slip,rt .Combining these and rearranging, we find that 0 = ∆ c + in/out,c + no-slip,c + no-slip,rb for the central streamtube, indicating that the recirculation zone is equivalent to the no-slip boundary on its periphery.If multiple recirculation zones separate the central streamtube from the no-slip boundary, this procedure can be repeated multiple times with the same result.
19The plane averaged turbulence intensity defined as Tu ≡ 2k /3/U, not averaging over u /U as one might expect.Rather than plane averaging the turbulence intensity directly, a new turbulence intensity is formed from the plane averaged turbulent kinetic energy per the definition of the turbulence intensity.This is done so that Tu where Tu in is the turbulence intensity at the nozzle inlet and Tu c is the turbulence intensity at the end of the contraction.Using equation19for the nozzle orifice with the linearization k 3/2 ≈ k • k This model has not been experimentally validated.We are not aware of any nozzle turbulence experimental data where all of the model inputs (e.g., including inlet turbulence intensity, Tu in ) are available.However, in contrast to Huh's model, this model is at least qualitatively correct in the limits of L 0 /d 0 .Comparison with data compiled by Klein [19, p. 246, fig.3] suggests that the model likely has the wrong concavity in k for small to moderate L 0 /d 0 and also underestimates the development length.A better model is needed.More detailed nozzle turbulence experiments are also needed to validate both the contraction and cylindrical segment parts of any nozzle turbulence model.