and the uniqueness of its solutions
Benoit Cushman-Roisin (benoit.r.roisin@dartmouth.edu)
Thayer School of Engineering, Dartmouth College, Hanover, NH 03755-8000, U.S.A. Tel.: +1-603-646 3248, Fax: +1-603-646 3856
Alastair D. Jenkins∗ (alastair.jenkins@bjerknes.uib.no)
Bjerknes Centre for Climate Research, Geophysical Institute, University of Bergen, All´egaten 70, 5007 Bergen, Norway. Tel.: +47-555 82632, Fax: +47-555 89883
Abstract. A non-local parameterisation of shear turbulence is proposed, which includes a dimensionless multiplicative constant as the sole tunable parameter. An- alytical and numerical solutions in the case of plane Couette flow exhibit sheared velocity profiles with logarithmic behaviour near the boundaries, and the classical logarithmic flow profile is reproduced for a semi-infinite domain. We also prove that the families of analytical solutions obtained are locally unique: if the velocity is a strictly-increasing function of the distance from the boundary, a small perturbation of the velocity profile must be of the same functional form as the basic flow.
Keywords: logarithmic boundary layer; non-local turbulence closure; turbulent Couette flow; transilient turbulence; uniqueness of solution.
∗ Corresponding Author
Couette flow; transilient turbulence; uniqueness of solution.
1. Introduction
In laminar flows, diffusion of momentum proceeds from molecular vis- cosity, which sets a connection between the velocity of a fluid particle and that of its immediate neighbours, and the the resulting term in the Navier–Stokes equations arises as the sum of spatial derivatives.
For planar unidirectional flow decelerating under the sole influence of viscosity ν, the velocityu(z, t) is governed by
∂u
∂t =ν∂2u
∂z2. (1)
In turbulent flow, however, fluctuations bring into contact fluid par- cels which would otherwise not be neighbours, and the interactions are no longer local. In such a non-local formulation, also termed anintegral closureortransilientmodel (Berkowicz and Prahm, 1980; Fiedler, 1984;
Stull, 1984, 1993; Nakayama et al., 1988; Nakayama and Bandou, 1995),
the spatial derivative of (1) is replaced by an integral over the entire domain:
∂u
∂t = Z h
0
a(ζ, z) [u(ζ, t)−u(z, t)]dζ, (2) witha(ζ, z) being a positive weighting factor measuring the importance of the momentum exchange between levels ζ and z.
We expect the weighting factor to be a function of the velocity difference|u(ζ, t)−u(z, t)| and distance |ζ−z| between the two levels.
Conservation of momentum is ensured if a(ζ, z) = a(z, ζ) for all ζ and z. Dimensional consistency requires a(ζ, z) to have dimensions L−1T−1, where Lis length andT is time. Ifa(ζ, z) depends only upon
|u(ζ, t)−u(z, t)| and |ζ−z|, dimensional analysis requires that it be equal to a constant multiplied by |u(ζ, t)−u(z, t)|/(ζ −z)2, and (2) becomes
∂u
∂t =A Z h
0
|u(ζ, t)−u(z, t)|
(ζ−z)2 [u(ζ, t)−u(z, t)]dζ, (3) where Ais a dimensionless constant.
PROPOSITION 1. Ifu(z, t) is differentiable with respect to bothzand t, the integral term on the right-hand-side of (3) respects conservation
of momentum, that is, (d/dt)R0hu(z, t)dz = 0.
Proof. Differentiability of u(z, t) implies the existence and bound- edness of ∂u/∂t and the integrand. Integrating (3) over the whole domain, and reversing the order of differentiation and integration on
the left-hand-side, we obtain d
dt Z h
0
u dz=A Z h
0
Z h 0
|u(ζ, t)−u(z, t)|
(ζ−z)2 [u(ζ, t)−u(z, t)]dζ dz. (4) The integrand on right-hand-side of (4) is antisymmetric with respect to the interchange ofζ andz. Interchanging the two integration variables is equivalent to reflecting the domain of integration [0, h]×[0, h] about the diagonal ζ =z. Hence the integrand is equal to minus itself, and the rate of change of total momentum (d/dt)R0hu(z, t)dz vanishes.
The integral term is also dissipative.
PROPOSITION 2. If (3) is satisfied, the total kinetic energy R0h12u2dz cannot increase with time.
Proof. Multiplying (3) byuand integration over the domain yields d
dt Z h
0
1
2u2dz = +A Z h
0
Z h
0
|u(ζ, t)−u(z, t)|
(ζ−z)2 [u(ζ, t)−u(z, t)]u(z, t)dζ dz
= −A Z h
0
Z h 0
|u(ζ, t)−u(z, t)|
(ζ−z)2 [u(ζ, t)−u(z, t)]u(ζ, t)dz dζ
= −A 2
Z h 0
Z h 0
|u(ζ, t)−u(z, t)|
(ζ−z)2 [u(ζ, t)−u(z, t)]2 dz dζ
≤ 0.
2. Stationary solutions
2.1. Semi-infinite domain
For h→ ∞we find that solutions of the form
u(z) =α+βlnz, (5)
whereαandβare arbitrary real constants, satisfy the time-independent form of (3):
Z ∞
0
|u(ζ, t)−u(z, t)|
(ζ−z)2 [u(ζ, t)−u(z, t)]dζ = 0. (6) Thus we may reproduce the mean flow in classical logarithmic turbulent boundary layer. Sinceα andβ are arbitrary, we may fit the usual form of the velocity profile
u(z) = (u∗/κ) ln(z/z0) (7) by setting β = u∗/κ and α = −(u∗/κ) lnz0, where u∗ = (τw/ρ)1/2 is the friction velocity, z0 is the aerodynamic roughness length, τw is the wall shear stress,ρis the fluid density, andκis von K´arm´an’s constant, whose experimentally-obtained value is approximately 0.4.
2.2. Turbulent plane Couette flow
For a finite domain, the corresponding family of solutions, to Z h
0
|u(ζ, t)−u(z, t)|
(ζ−z)2 [u(ζ, t)−u(z, t)] dζ= 0, (8) is
u(z) =α+βtanh−1 2z
h −1
=α+β 2 ln z
h−z (9)
[the arbitrary constants α and β being different from those in (5)], so that the velocity profile is logarithmic in the vicinity of each boundary.
2.3. Uniqueness of solutions
We have not so far ascertained whether there are solutions to (6) and (8) which differ from (5) and (9), respectively. It turns out that we can prove that if such other solutions exist, they cannot be obtained by continuous perturbations of (5) or (9).
Semi-infinite domain
PROPOSITION 3. A linear perturbation of the solution (5),
u(z) =α+βlnz+η(z), (10) with η(z) twice differentiable andu(z) strictly increasing or decreasing, satisfies Eq. (6) to O() if and only if η(z) is of the form γ1+γ2lnz.
Proof.Foru(z) strictly increasing or decreasing, Eq. (6) is equivalent to
Z z 0
u(ζ)−u(z) ζ−z
2
dζ = Z ∞
z
u(ζ)−u(z) ζ−z
2
dζ.
Substituting X = lnz, v(X) =u(z), andξ = ln|ζ/z|, we obtain Z ∞
0
[2v(X)−v(X −ξ)−v(X+ξ)] [v(X +ξ)−v(X −ξ)]
× e−ξ
(1−e−ξ)2 dξ= 0. (11) The formula (10) is equivalent to v(X) = α+βX+θ(X), where θ(X) =η(z) is a function to be determined. ToO(), we have
Z ∞
0
[2θ(X)−θ(X−ξ)−θ(X+ξ)] ξe−ξ
(1−e−ξ)2dξ= 0. (12) If we represent θ as a Fourier integral, θ(X) = R−∞∞ θ(k)eˆ ikXdk, we obtain
Z ∞
0
Z ∞
−∞
θ(k)eˆ ikX(1−coskξ)dk
ξe−ξ
(1−e−ξ)2 dξ= 0. (13) Reversing the order of integration, we then have
Z ∞
−∞
θ(k)ˆ
"Z ∞ 0
(1−coskξ) ξe−ξ (1−e−ξ)2 dξ
#
eikXdk= 0, (14) which we may write as
Z ∞
−∞
θc00(k)W(k)eikXdk= 0, (15) where cθ00(k) is the Fourier transform of the second derivative ofθ(X), andW(k) =k−2R0∞(1−coskξ)1−e−ξ
−2
ξe−ξdξis finite and greater than zero for all real k.
The Fourier transform of Eq. (15) iscθ00(k)W(k) = 0, so thatθc00(k)=0 for all real k. Taking the inverse Fourier transform, we deduce that θ00(X) = 0 for all real X, so that η(z) = γ1 +γ2lnz with γ1 and γ2
constant.
Turbulent Couette flow
For a finite domain, we may, without loss of generality, re-scalezso that the boundaries of the domain are z = −1 and z = +1. Equation (8) then becomes
Z 1
−1
|u(ζ, t)−u(z, t)|
(ζ−z)2 [u(ζ, t)−u(z, t)] dζ= 0, (16) which has solutions of the form
u(z) =α+βtanh−1z. (17)
PROPOSITION 4. A linear perturbation of the solution (17),
u(z) =α+βtanh−1z+η(z), (18)
with η(z) twice differentiable and u(z) strictly increasing or decreas- ing, satisfies Eq. (16) to O() if and only if η(z) is of the form γ1 + γ2tanh−1z.
Proof. Foru(z) strictly increasing or decreasing, Eq. (16) is equiva- lent to
Z z
−1
u(ζ)−u(z) ζ−z
2
dζ= Z 1
z
u(ζ)−u(z) ζ−z
2
dζ.
SubstitutingX = tanh−1z,v(X) =u(z), andξ =tanh−1ζ−tanh−1z, we obtain
Z ∞
0
[2v(X)−v(X −ξ)−v(X+ξ)] [v(X +ξ)−v(X −ξ)]
× 1
sinh2ξ dξ= 0. (19) Using the same reasoning as in the proof of Proposition 3, withv(X) = α+βX+θ(X), we find that the Fourier transformcθ00(k) of the second derivative of θ satisfies
Z ∞
−∞
θc00(k)W(k)eikXdk= 0, (20) with W(k) = k−2R0∞(1−coskξ) (sinhξ)−2ξ dξ also being finite and greater than zero for all real k. We deduce that θ00(X) = 0 for all real X, so that η(z) =γ1+γ2tanh−1zwith γ1 and γ2 constant.
Summary of uniqueness property
We have now proved that there are no strictly increasing or decreasing solutions of (6) and (8) which are close to the solutions given in (5) and (9), respectively, which are not also of the same form, so that
there is no continuous family of solutions containing the given solution and another type of solution. The question remains open whether other types of solution to the nonlinear equations (6) and (8) exist, although it appears to be physically unlikely.
The monotonicity condition appears to be a physically reasonable requirement for a stationary flow:
PROPOSITION 5. A non-constant solution to the either of the sta- tionary-flow equations (6) or (8) cannot have a finite absolute maxi- mum or minimum at any point in the domain or on the boundary.
Proof. In either case, the integrand at z = zm, where zm is the point at which u takes its absolute maximum or minimum value, will be bounded away from zero, so the equation cannot be satisfied.
3. Numerical discretisation
The boundary conditions of the plane Couette flow are u(z = 0) = 0 and u(z = h) = U. Non-dimensionalising z by h, u by U and t by h/AU, Equation (3) becomes
∂u
∂t = Z 1
0
|u(ζ, t)−u(z, t)|
(ζ−z)2 [u(ζ, t)−u(z, t)]dζ, (21)
which we solve numerically, starting with the non-dimensional laminar- flow profileu(z) =z, while imposing boundary conditionsu(z= 0) = 0 and u(z = 1) = 1. Time integration proceeds until a steady state is reached, i.e., until the integral term is made to vanish for every z.
Note that the model is totally deprived of every parameter, save those introduced by the numerical discretisation.
The solid line in Fig. 1a shows the numerical solution obtained with n= 100 intervals of discretisation. A profile of the form
u(z) = 1 2 −1
2
ln[(z+d)/(1 +d−z)]
ln[d/(1 +d)] , (22)
which corresponds to (9) with the addition of a roughness height d [d = 0.00442 in the present case] yields an almost perfect fit. The procedure was repeated many times with various levels of discretisation.
In all cases, almost perfect fits were obtained, by adjusting the value of d.
While Eq. (21) does not contain a single parameter, the numerical solution does by way of the grid size, ∆z= 1/n. A series of numerical experiments with various levels of numerical resolution (Fig. 2) reveals that the value of the roughness height d adopted by the solution is a direct function of the grid size ∆z. When molecular viscosity is added and the spatial discretisation is sufficient to resolve the viscous
Figure 1. (a) Velocity profile in a turbulent Couette flow withn= 100 (solid line) and fit of two logarithmic curves (dashed lines); (b) Same velocity profile re-plotted versus the logarithm of distance from the bottom and versus the logarithm of the distance from the top.
sub-layer, the solution becomes independent of the level of numerical discretisation (see Section 5 below).
Figure 2. Variation of the roughness height d with the numerical resolution n= 1/∆z.
4. Turbulent stress and connection with the von K´arm´an constant
With the integral term in Eq. 3, the momentum-exchange term does not appear as the spatial derivative of a stress (∂τ /∂z). However, a physical stress must exist, especially in the vicinity of the boundaries.
In the context of this formalism, the way to determine the stress τ(z, t) at levelzinside the domain is to consider the flux of momentum that crosses it, say from below to above (Stull, 1984, 1993). Consider a level z inside the domain (0< z < h) and two arbitrary levelsζ andζ0 one on each side of z (0≤ζ < z < ζ0 ≤h). The momentum exchange a(ζ0, ζ)[u(ζ0, t)−u(ζ, t)] betweenζandζ0contributes to the momentum flux through the intermediate level z. Summing all such exchanges, i.e., integrating over all possible values of ζ below z and of ζ0 above z, we obtain the overall momentum flux at level z:
τ(z, t)
ρ =A
Z z 0 dζ
Z h z
|u(ζ0, t)−u(ζ, t)|
(ζ0−ζ)2 [u(ζ0, t)−u(ζ, t)]dζ0, (23) where τ is the stress and ρ the fluid density. It can be shown that the z-derivative of this expression returns the integral of Eq. (3), and that
it yields a value independent of z in the case of the logarithmic profile (h→ ∞). We shall thus adopt this definition for the turbulent stress, which, we note, depends quadratically on the velocity.
For a semi-infinite domain, the steady-state solution (5) may be written asu(z) =Uln(z/d), and therefore the corresponding stress (23) is uniform in z. Indeed, withζ =az and ζ0=z/b,
τ(z)
ρ = A
Z z 0
dζ Z ∞
z
U2[ln(ζ0/ζ)]2 (ζ0−ζ)2 dζ0
= AU2 Z 1
0
Z 1
0
[ln(ab)]2 (1−ab)2da db
= 6AU2 X∞ n=1
1
n3 = 7.2123AU2 (24)
is found to be independent of z and thus equal to τw/ρ, the value at the wall.
In terms of the friction velocity defined as u∗ =pτw/ρ, laboratory experiments give U = u∗/κ, with κ = 0.4 being the von K´arm´an constant. Using (24), we then find the following relation between A and κ:
A= κ2
7.2123 = 0.02218 = 1
45.08. (25)
For κ= 0.41, an oft-quoted alternative value, A= 0.02331 = 1/42.90.
5. Inclusion of viscosity and comparison with laboratory data
We now combine the turbulent term of Eq. (2) with the original viscous term of (1) in order to overcome the dependency of the solution on the
level of numerical resolution, and our equation of motion becomes
∂u
∂t =A Z h
0
|u(ζ, t)−u(z, t)|
(ζ−z)2 [u(ζ, t)−u(z, t)]dζ+ν∂2u
∂z2, (26) which, after scalingu byU,z byh, andtbyh/U, can be re-written as
∂u
∂t =A Z 1
0
|u(ζ, t)−u(z, t)|
(ζ−z)2 [u(ζ, t)−u(z, t)]dζ+ 1 Re
∂2u
∂z2, (27) where Re = U h/ν is the Reynolds number. Boundary conditions are u(z= 0) = 0 andu(z= 1) = 1. For the value ofAgiven in (25) and for a series ofRevalues, numerical solutions were obtained with increasing levels of numerical resolution until convergence. Figure 3 displays the viscous, turbulent and total stress profiles in the case Re= 2500.
From every solution, two parameters were extracted:Sthe mid slope (S=∂u/∂z at z= 1/2) andu∗ the friction velocity (u2∗=Re−1∂u/∂z atz= 0). Then, for the purpose of comparison with existing laboratory data (Robertson and Johnson, 1970), a skin-friction coefficient cf was derived from the friction velocity according to
cf = 4τ(z= 0)
1
2ρU2 = 8u2∗. (28)
Figure 4 compares the present model results with the laboratory data reported by Robertson and Johnson (1970). We note that the model values predicted for the mid-slope and for the skin-friction coef- ficient (dashed lines) are significantly overestimated.
Figure 3. Viscous (dotted line), turbulent (dashed line), and total (solid line) shear stresses in a turbulent Couette flow when molecular viscosity is included (Re=hU/ν= 2500,n= 400).
Remembering that in the vicinity of each wall there exists a viscous sub-layer within which the flow is purely laminar, we then suppressed the non-local term involving points near each boundary by setting to zero the weighting factor a(ζ, z) wheneverζ orzfell within a distance of 5ν/u∗ [ = 5(Reu∗)−1 after scaling] from either boundary. The factor 5 is justified as the commonly-accepted value for the thickness of the
Figure 4. Variation of (a) mid-slope and (b) skin-friction coefficient with the Reynolds number. The data points are those reported by Robertson and Johnson (1970), while the dashed curves are the non-truncated model predictions and the solid curves are the truncated-model predictions.
viscous sub-layer (Hinze, 1975, p. 628), and the value of u∗ was sim- ply determined in every case from u2∗ =Re−1(∂u/∂z)z=0, which holds because the flow is laminar there.
The solid lines in Fig. 4 show the variation of the mid-slope and skin- friction coefficients with the Reynolds number. Obviously, the turbulent model does not include the transition from laminarity to turbulence.
Otherwise, the good agreement is very encouraging. The mid-slope values are still slightly overestimated, but the model prediction for skin-friction coefficient is excellent in the turbulent regime (Re>4000).
6. Conclusions
The parameterisation of shear turbulence as a weighted integral of velocity differences successfully reproduces the logarithmic profile of near-wall shear turbulence. To the authors’ knowledge, this is the first theory capable of doing so without the explicit introduction of a mixing length or other turbulent scale1. Questions arising as to whether solu-
1 The theory of A. E. Perry and collaborators [see, e.g., Perry and Chong (1982)]
reproduces the logarithmic velocity profile without resorting to a mixing length pro- portional to the distance from the wall, but requires prior assumptions on turbulent scales.
tions to the integral equations employed in this parameterisation are unique are partially answered by the uniqueness proofs in Section 2.3.
We show, under suitable conditions of differentiability and monotonic- ity, that any continuous perturbation of the logarithmic and inverse hyperbolic tangent solutions of (5) and (9), respectively, cannot change the solution to one of different functional form.
There is therefore a strong motivation to try to rationalise the in- tegral term from the Navier-Stokes equations according to a certain phenomenology, and to exploit the model further in order to evaluate its overall possibilities and limitations.
Acknowledgements
Gratitude goes to Professor Akihiko Nakayama of Kobe University (Japan) and to Dr. Vlado Malaˇciˇc of the Marine Biology Station in Piran (Slovenia) for numerous and enlightening discussions. We also wish to acknowledge the courtesy of Professor Roland Stull of the Uni- versity of British Columbia (Canada) for sending a set of his articles on transilient turbulence parameterisation. Support for this research was provided by the U.S. Office of Naval Research under grant No. N00014- 02-1-0065 to Dartmouth College, and by the Research Council of Nor-
way under projects 139815/720 and 155923/700. This is publication No. A0022 of the Bjerknes Centre for Climate Research.
References
Berkowicz, R. and L. P. Prahm: 1980, ‘On the Spectral Turbulent Diffusivity Theory for Homogeneous Turbulence’. J. Fluid Mech.100, 433–448.
Fiedler, B. H.: 1984, ‘An Integral Closure Model for the Vertical Turbulent Flux of a Scalar in a Mixed Layer’. J. Atmos. Sci.41(4), 674–680.
Hinze, J. O.: 1975, Turbulence. New York: McGraw-Hill.
Nakayama, A. and M. Bandou: 1995, ‘Incorporation of Nonlocal Effects in Two- Equation Models’. In: Proc. Int. Symp. Math. Modelling of Turbulent Flows, Tokyo. pp. 25–30.
Nakayama, A., H. D. Nguyen, and M. B. Daif: 1988, ‘Nonlocal Diffusion Model in Turbulent Boundary Layer’. In:Proc. Third. Int. Symp. Refined Flow Modelling and Turbulence Measurements. pp. 305–312, Int. Assoc. Hydraulic Res.
Perry, A. E. and M. S. Chong: 1982, ‘On the Mechanism of Wall Turbulence’. J.
Fluid Mech.119, 173–217.
Robertson, J. M. and H. F. Johnson: 1970, ‘Turbulence Structure in Plane Couette Flow’. Proc. ASCE, J. Engng Mech.96, 1171–1182.
Stull, R. B.: 1984, ‘Transilient Turbulence Theory. Part I: The Concept of Eddy- Mixing across Finite Distances’. J. Atmos. Sci.41, 3351–3367.
Stull, R. B.: 1993, ‘Review of Non-Local Mixing in Turbulent Atmospheres:
Transilient Turbulence Theory’. Bound.-Layer Meteorol.62, 21–96.
