HYDRODYNAMIC STABILITY MATM028 Matt Turner October 2011 Books: 1. Drazin: “Introduction to hydrodynamic stability” 2. Drazin & Reid: “Hydrodynamic stability” 3. Schmid & Henningson: Stability and transition in shear flows Aims of the Course 1. To explain phenomena such as why rising smoke changes from a smooth flow into an erratic flow or why a smooth river flow emerging from a river mouth can ultimately produce vortices in the sea/ocean. 2. For boundary layer flows close to the surface of a rigid body (e.g. flow over a wing), we explain how the flow can be laminar (i.e. parallel to the surface) near to the front edge, but further downstream the flow can be turbulent. 3. Understand fluid flow stability, i.e. whether small disturbances in the flow grow (UNSTABLE) or decay (STABLE) in either space or time.
i
Contents 1 Important concepts from Hydrodynamics
1
2 Kelvin–Helmholtz Instability (An inviscid instability) 2.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Normal modes . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 1 dimensional disturbances . . . . . . . . . . . . . 2.2.2 2 dimensional disturbances . . . . . . . . . . . . . 2.3 Kelvin–Helmholtz stability using normal modes . . . . . 2.3.1 Special cases . . . . . . . . . . . . . . . . . . . . . 2.4 Effect of surface tension on K–H instability . . . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
4 . 4 . 7 . 7 . 8 . 9 . 11 . 12
3 Inviscid Stability 3.1 General theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Piecewise Linear Base Flows . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Continuous velocity profiles . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Criteria for inviscid instability . . . . . . . . . . . . . . . . . . 3.3.2 Numerical calculation of eigenvalues c(K) for continuous U (y) 3.3.3 Critical layers in inviscid flows . . . . . . . . . . . . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
13 13 15 15 17 21 21 23 24
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
4 Viscous stability 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.1 Boundary conditions for viscous flow . . . . . . . . . . . . . . 4.1.2 Strictly parallel base flows . . . . . . . . . . . . . . . . . . . . 4.2 Orr–Sommerfeld Equation . . . . . . . . . . . . . . . . . . . . . . . . 4.2.1 Squire’s Transformation . . . . . . . . . . . . . . . . . . . . . 4.3 Squire’s Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Formal theory of eigenvalues of OSE . . . . . . . . . . . . . . . . . . 4.4.1 Finding eigenvalues numerically using Chebyshev polynomials 4.4.2 Summary of numerical results for OSE . . . . . . . . . . . . . 4.5 Viscous stability analysis for Poiseuille flow . . . . . . . . . . . . . . . 4.5.1 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . 4.5.2 Asymptotic theory for ODEs . . . . . . . . . . . . . . . . . . . 4.5.3 Asymptotic theory of Orr–Sommerfeld equation for large R . . 4.5.4 Summary of solutions of the Orr–Sommerfeld equation . . . . 4.5.5 Calculating neutral stability curve . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
29 29 29 29 30 31 33 34 35 36 36 37 38 40 44 47
5 Boundary layers 5.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . 5.2 General theory for R 1 . . . . . . . . . . . . . . . 5.3 Uniform flow past a flat plate . . . . . . . . . . . . . 5.3.1 Notes on Blasius solution . . . . . . . . . . . . 5.4 Flow past flat plate with arbitrary slip velocity U (x) 5.4.1 Falkner–Skan solutions . . . . . . . . . . . . . 5.4.2 Numerical results for Falkner–Skan . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
49 49 49 51 52 55 57 58
ii
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
. . . . . . .
5.5 5.6
5.4.3 Summary of Falkner–Skan profiles . . . . . . . . . . . . . . . . . . . . . . . . 59 Stability of boundary layers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 Comparison of experiment and theory for neutral curve for flat plate . . . . . . . . . 62
6 Laminar/Turbulent transition on a flat 6.1 Introduction . . . . . . . . . . . . . . . 6.2 Stability theory . . . . . . . . . . . . . 6.3 Receptivity theory . . . . . . . . . . . 6.4 Leading edge receptivity . . . . . . . . 6.5 Receptivity due to roughness element . 6.6 Importance of stability properties . . .
plate . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
67 67 67 69 70 70 72
A MATLAB Worksheets A.1 Kelvin–Helmholtz instability . . . . . . . . . . . . . . . . . A.1.1 Kelvin–Helmholtz general stability . . . . . . . . . A.2 Inviscid Stability . . . . . . . . . . . . . . . . . . . . . . . A.2.1 Inviscid stability of mixing layer . . . . . . . . . . . A.2.2 Unbounded piecewise linear mixing layer . . . . . . A.2.3 Two part piecewise linear boundary layer flow . . . A.2.4 Blasius boundary layer solution . . . . . . . . . . . A.2.5 Stability analysis for U = tanh y - numerical results A.3 Viscous Stability . . . . . . . . . . . . . . . . . . . . . . . A.3.1 Eigenvalues for Poiseuille flow. . . . . . . . . . . . . A.3.2 Asymptotic solution to ODE . . . . . . . . . . . . . A.4 Boundary Layers . . . . . . . . . . . . . . . . . . . . . . . A.4.1 Falkner–Skan velocity profiles . . . . . . . . . . . . A.5 Transition . . . . . . . . . . . . . . . . . . . . . . . . . . . A.5.1 Summary of Goldstein (1983) JFM . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
74 74 74 76 76 77 79 83 85 90 90 92 94 94 96 96
iii
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
1
Important concepts from Hydrodynamics
Here we review a few concepts from hydrodynamics that you should know!!! or at least have seen. Definitions We define the fluid velocity as u(x,t)=u(x,t)i+v(x,t)j+w(x,t)k, where x=(x, y, z) in Cartesian coordinates. A flow is said to be incompressible if ∇ · u = 0, where ∇ ≡ ∂/∂xi + ∂/∂yj + ∂/∂zk. The vorticity of a flow field is denoted by ω = ∇ ∧ u. A flow is said to be irrotational if ω = 0. If ω = 0, then we can find a scalar potential φ(x, t) such that u = ∇φ, where φ is known as the velocity potential. If the flow is incompressible and irrotational then ∇2 φ = 0. For 2–D flows which are incompressible, we can define a streamfunction ψ(x, y) such that u = (ψy , −ψx ). Navier–Stokes equations For incompressible flows ρ Du =− Dt
F, ∇p + µ∇2 u + Effect of Effect of External forces pressure viscosity e.g. Gravity
where
Du ∂u = + (u · ∇)u, Dt ∂t is the material or convective derivative. Ignore external forces and non–dimensionalise the equation. Assume the problem has a characteristic scale for fluid velocity, U , and a length scale L. e.g. a ball moving through the air x = Lˆ x, dimensional non − dimensional n pos vector posn vector ∇ = ∂/∂xi + ∂/∂yj + ∂/∂zk = 1
ˆ, u = Uu
1ˆ ∇, L
U L Figure 1: Characteristic length and velocity scales
ˆ U2 ∂u 1ˆ µU ˆ ˆ, ρ U + (ˆ u · ∇)ˆ u = − ∇p + 2 ∇ˆ2 u ∂t L L L t and p are still dimensional quantities. Divide by ρU 2 /L ˆ L ∂u ˆ u = − 1 ∇p ˆ + µ ∇ˆ2 u ˆ. + (ˆ u · ∇)ˆ 2 U ∂t ρU ρU L Non–dimensionalise p and t by t=
Lˆ t U
p = ρU 2 pˆ,
so
ˆ ∂u ˆ u = −∇ˆ ˆ p + 1 ∇ˆ2 u ˆ, + (ˆ u · ∇)ˆ ˆ Re ∂t where Re = ρU L/µ is the Reynolds number. If Re 1 =⇒ viscous term is small. The limit Re → ∞ known as inviscid limit. Streamlines Lines instantaneously parallel to velocity of particles at all points. For steady flows – streamline is path taken by fluid particle. For 2D flows – streamlines are lines ψ = constant. Inviscid flows For steady flows with gravitational field Bernoulli’s law relates pressure and velocity 1 p + ρu2 + ρgz = constant along streamline. 2 For unsteady flows which are irrotational 1 ∂φ p + ρ(∇φ)2 + ρ + ρgz = constant throughout fluid. 2 ∂t 2
u
Figure 2: Streamlines in the flow.
3
2 2.1
Kelvin–Helmholtz Instability (An inviscid instability) Formulation
Consider 2 fluids ρ = ρ2 ρ = ρ1
and U = U2 i z > 0 and U = U1 i z < 0 z
ρ
y
U2
2
x
z=0 U1
ρ1
Figure 3: Kelvin–Helmholtz instability setup U1 and U2 are constants. If we ignore viscosity, then this flow is incompressible and satisfies the Navier–Stokes equations if the pressure is p0 − ρ2 gz p0 − ρ1 gz
z > 0, z < 0.
But we don’t see this flow in experiments. e.g. Vortices develop at the interface between the two U2
U1
Figure 4: Kelvin–Helmholtz instability experiment fluids. Now we perturb the flow slightly, e.g. the velocity potential becomes φ1 = U1 x + φ01 (x, t) φ2 = U2 x + φ02 (x, t) Position of interface between the 2 fluids z = Z(x, y, t). 4
z < Z, z > Z.
φ=φ 2 z=0 z= ε Z
φ=φ 1
Figure 5: Kelvin–Helmholtz instability with perturbed conditions.
Assume the disturbance is restricted to the neighborhood of the interface. i.e. u → U2 i as z → ∞, so
∇φ02 → 0 as z → ∞ ∇φ01 → 0 as z → −∞
(2.1.1)
Irrotational and incompressible flow =⇒ ∇2 φ1 = ∇2 φ2 = 0, so
∇2 φ01 = 0 ∇2 φ02 = 0
(2.1.2)
A fluid particle on the interface must stay on the interface u2 z= ε Z
Figure 6: Fluid particle remains on the interface D(Z) , Dt ∂ = (Z) + u · ∇(Z), ∂t ∂Z = + (U2 i + ∇φ02 ) · ∇Z, ∂t ∂Z ∂Z = + U2 + ∇φ02 · ∇Z, at z = Z. ∂t ∂x
u2 = ∂φ2 ∂z ∂φ02 ∂z ∂φ02 ∂z
Similarly on the other side of the interface ∂φ01 ∂Z ∂Z = + U1 + ∇φ01 · ∇Z, ∂z ∂t ∂x Known as the kinematic condition on the interface. 5
at z = Z.
(2.1.3)
Final condition is that pressure is continuous across the interface. ∂φ1 1 2 (∇φ1 ) + + gz z < Z, p1 = C 1 − ρ 1 2 ∂t 1 ∂φ2 2 p2 = C 2 − ρ 2 (∇φ2 ) + + gz z > Z, 2 ∂t where C1 and C2 are constants in space. Pressure is continuous, 1 ∂φ1 1 ∂φ2 2 2 C1 − ρ 1 (∇φ1 ) + + gz = C2 − ρ2 (∇φ2 ) + + gz . 2 ∂t 2 ∂t
(2.1.4)
Now consider small perturbations (i.e. 1) and linearise the equations. NOTE: ∂φ1 + O(2 ), φ1 (z = Z) = φ1 (z = 0) + Z ∂z z=0 = φ1 (0) + O(), and
∂φ1 ∂φ1 = + O(). ∂z z=Z ∂z z=0
Linearising equations (2.1.3) ∂φ02 ∂z ∂φ01 ∂z
= =
∂Z ∂t ∂Z ∂t
+ U2 ∂Z + O() ∂x ∂Z + U1 ∂x + O()
) .
(2.1.5)
Next linearise (2.1.4) (∇φ1 )2 = (U1 i + ∇φ01 )2 , = U12 + 2U1 i · ∇φ01 + O(2 ), ∂φ0 = U12 + 2U1 1 + O(2 ), ∂x 0 0 1 1 ∂φ1 ∂φ1 ∂φ02 ∂φ02 2 2 + + gZ = C2 − ρ2 U2 − ρ2 U2 + + gZ . C1 − ρ1 U1 − ρ1 U1 2 ∂x ∂t 2 ∂x ∂t Hence equating the coeffs of O(0 ) : O() :
1 1 C1 − ρ1 U12 = C2 − ρ2 U22 2 2 0 ∂φ1 ∂φ01 ∂φ02 ∂φ02 ρ1 U1 + + gZ = ρ2 U2 + + gZ ∂x ∂t ∂x ∂t
at z = 0. (2.1.6)
Can now solve set of eqns (2.1.5) and (2.1.6) together with (2.1.1) and (2.1.2) to solve for φ01 (x, y, z, t), φ02 (x, y, z, t), and Z(x, y, t).
NB φ01 = φ02 = Z = 0 is a solution. To look for a non–zero solution =⇒ Eigenvalue problem. 6
2.2 2.2.1
Normal modes 1 dimensional disturbances
Consider f (x, t) = cos(kx − ωt), = Re ei(kx−ωt) , k and ω are real constants. t=0:
f(x,0)
1
x −1 x=2 π/k Figure 7: Wave disturbance, Wavelength= 2π/k, Wavenumber= k. t = t0 > 0 : Hence f represents a wave travelling in x−direction with speed ω/k. 1 x −1 x= ω t 0/k
Figure 8: Wave disturbance at t = t0 . Temporal behaviour x = 0 Temporal frequency ω/2π, Period 2π/ω. If we now consider ω complex, ω = e = f (x, t) = = i(kx−ωt)
ωr + iωi ei(kx−ωr t )eωi t , Re ei(kx−ωt) , eωi t cos(kx − ωr t).
If ωi > 0, f (x, t) is a wave growing in amplitude with t. If ωi < 0, f (x, t) decays with t. If ωi = 0, f (x, t) stays the same amplitude. 7
1 t −1 t=2 π/ ω Figure 9: Wave disturbance at x = 0. 2.2.2
2 dimensional disturbances
Consider f (x, y, t) = ei(κ·x−ωt) , where x = xi + yj, κ = ki + lj, f = ei(kx+ly−ωt) . At t = 0, plot lines with f = 1 So, kx + ly = 2nπ. Interpret f (x, y, t) as wave propagating in
k
Figure 10: Wave disturbance at t = 0 in 2D. direction κ. As before, if ω is complex then wave grows in amplitude with t if ωi > 0. Alternative notation ω , k f = ei(kx+ly−ckt) . c =
8
U2
ρ2 z=0
U1
ρ1
Figure 11: Kelvin–Helmholtz instability
2.3
Kelvin–Helmholtz stability using normal modes
Consider a perturbation of the interface z = Z(x, y, t) composed of normal modes either as a sum X Z(x, y, t) = ck,l ei(κ·x−ωt) , k,l
or as an integral Z Z Z(x, y, t) =
c(k, l)ei(κ·x−ωt) dk dl.
We have linearised the governing equations, so can consider each normal mode separately. Then find how ω depends on κ (i.e. on k and l). Let Z = ei(κ·x−ωt) , κ = ki + lj. φ01 (x, y, z, t) = f1 (z)ei(κ·x−ωt) φ02 (x, y, z, t) = f2 (z)ei(κ·x−ωt)
z < 0, z > 0.
Substitute φ01 into linearised (2.1.2), ∇2 φ01 = 0, 2 ∂ ∂2 ∂2 + 2 + 2 f1 (z)ei(κ·x−ωt) = 0, 2 ∂x ∂y ∂z d2 f1 i(κ·x−ωt) 2 2 (ik) f1 + (il) f1 + e = 0, dz 2 df1 − K 2 f1 = 0, 2 dz where K = (k 2 + l2 )1/2 > 0. Using the linearised boundary condition (2.1.1), ∇φ01 → 0 as z → −∞, ∴ f1 (z) = A1 e−Kz + B1 eKz , = B1 eKz . Similarly φ02 = f2 (z)ei(κ·x−ωt) , f2 (z) = A2 e−Kz , 9
A2 and B1 will be determined from eqns (2.1.6) and (2.1.5). Subst into (2.1.5) ∂φ01 ∂Z ∂Z + U1 , = ∂z z=0 ∂t ∂x f10 (0)ei(κ·x−ωt) = −iωei(κ·x−ωt) + U1 ikei(κ·x−ωt) B1 K = i(kU1 − ω). From eqn for φ02 , −A2 K = i(kU2 − ω). Equation (2.1.6) gives ∂φ02 ∂φ02 ∂φ01 ∂φ01 + + gZ = ρ2 U2 + + gZ ρ1 U1 ∂x ∂t ∂x ∂t ρ1 (U1 ikB1 − iωB1 + g) = ρ2 (U2 ikA2 − iωA2 + g) , i(kU1 − ω)ρ1 B1 − i(kU2 − ω)ρ2 A2 = (ρ2 − ρ1 )g.
at z = 0,
Substitute for A2 and B1 to give a quadratic for ω (ρ1 + ρ2 )ω 2 − 2k(U1 ρ1 + U2 ρ2 )ω + k 2 (ρ2 U22 + ρ1 U12 ) + (ρ2 − ρ1 )gK = 0, (check the arithmatic). Now solve for ω(k, l). ω=
S 1/2 k(U2 ρ2 + U1 ρ1 ) ± , ρ1 + ρ2 2(ρ2 + ρ1 )
where S = B 2 − 4AC
B = −2k(U1 ρ1 + U2 ρ2 ) etc.
Rearranging ω = αk ± β
p γK − k 2 ,
where U2 ρ2 + U1 ρ1 , ρ1 + ρ2 √ ρ1 ρ2 β = |U1 − U2 | ρ1 + ρ2 (ρ21 − ρ22 )g . γ = ρ1 ρ2 (U1 − U2 )2
α =
If γK − k 2 ≥ 0, roots ω are both real. Hence, disturbance neither grows or decays (known as NEUTRALLY STABLE). If γK − k 2 < 0 ω = αk ± β(k 2 − γK)1/2 i. Hence two possible behavious, one mode grows with time and one decays. If any mode grows with time then disturbance will eventually become large and system is said to be UNSTABLE to the particular k and l of normal mode. 10
2.3.1
Special cases
1: ρ1 < ρ2 Recall that K = (k 2 + l2 )1/2 > 0, so if ρ1 < ρ2 , then γ < 0. γK < 0 < k 2 , ∴ Unstable to all disturbances.
Heavy
Light Figure 12: ρ1 < ρ2 : Expect instability. 2: ρ1 = ρ2
ω = αk ± βki Unstable to all disturbances (i.e. to all k). Growth rate = βk = 2πβ/λ where λ is wavelength (λ = 2π/k). Shortest waves are most unstable. 3: If U1 = −U and U2 = +U
ω = ±U ki Zero phase speed, growth rate = U k = 2πU/λ. (see this flow again in §3). 4: ρ1 > ρ2 Unstable if k 2 > γK k 4 > γ 2 K 2 = γ 2 (k 2 + l2 ) i.e. unstable if 1 k2 > γ 2 1 + 2
s
4l2 1+ 2 γ
!
In (k, l)−plane Stable for small k, i.e. for long wavelengths. See appendix A.1.1 for MATLAB worksheet.
11
l STABLE UNSTABLE k
k= γ
Figure 13: ρ1 > ρ2 n p
2
p
1
surface
Figure 14: Surface tension at a fluid interface.
2.4
Effect of surface tension on K–H instability
The surface tension at a fluid interface satisfies the equation p1 − p2 = σ∇ · n, where n is normal to the surface and σ is coefficient of surface tension, where p1 and p2 are pressures on either side of the interface. For surface z = Z(x, y), the surface is ψ = z − Z = constant. ∇ψ = (−Zx , −Zy , 1), ∇ψ (−Zx , −Zy , 1) . n = = |∇ψ| (1 + 2 Zx2 + 2 Zy2 )1/2 Linearising for small n = (−Zx , −Zy , 1) + O(2 ), ∇ · n = −Zxx − Zyy + 0 + O(2 ), p1 − p2 = −σ(Zxx + Zyy ). Hence this modifies pressure equation (2.1.6).
12
y
U(y)
x
Figure 15: Channel Flow
3 3.1
Inviscid Stability General theory
Consider the stability of basic flow U = U (y)i. Stability to disturbances of the form ei(kx+lz−kct)
c=
ω , k
in the absence of gravity. NB Change of notation: Here base velocity is a function of y and disturbance is a wave in the (x, z)– plane cf §2, where U was dependent on z and disturbance in (x, y)–plane. Inviscid Navier–Stokes Eqn: Let total velocity be uT , ∂uT + uT · ∇uT = −∇pT , ρ ∂t also ∇ · uT = 0
(Incompressible flow).
Consider small perturbation to the base flow uT = U + u,
pT = P + p,
∂u 2 ρ + U · ∇U + U · ∇(u) + u · ∇U + O( ) = −∇P − ∇p. ∂t
O() :
ρ
O(0 ) : ρU · ∇U = −∇P, U = U (y)i, ∂ ρ U (U (y)i) = −∇P, ∂x ∇P =⇒ P = constant. = 0
∂u + U · ∇u + u · ∇U ∂t
13
= −∇p.
Take ρ constant ∂u ∂u + U (y) + u · ∇(U (y)i) = −∇ ∂t ∂x
p . q
Let u = ui + vj + wk and q = p/ρ. ∂u ∂u ∂U +U +v ∂t ∂x ∂y ∂v ∂v +U ∂t ∂x ∂w ∂w +U ∂t ∂x
∂q , ∂x ∂q = − , ∂y ∂q = − . ∂z = −
(3.1.7) (3.1.8) (3.1.9)
Also, incompressible, so continuity equation gives ∇ · u = 0, ∂u ∂v ∂w + + = 0. ∂x ∂y ∂z Take
(3.1.10)
∂ ∂ ∂ (3.1.7) + (3.1.8) + (3.1.9), ∂x ∂y ∂z ∂ ∂ ∂U ∂v ∂v ∂U (∇ · u) + U (∇ · u) + + = −∇2 q, ∂t ∂x ∂y ∂x ∂x ∂y
Using ∇ · u = 0 2U 0 (y)
∂v = −∇2 q. ∂x
Eliminate q (i.e. pressure) from (3.1.8) ∂ ∂ ∂ 2v ∂v ∂ (∇2 v) + U (∇2 v) + 2U 0 + U 00 = − (∇2 q), ∂t ∂x ∂x∂y ∂x ∂y where ∂ (U vx ) = U 0 vx + U vxy , ∂y ∂2 (U vx ) = U 00 vx + 2U 0 vxy + U vxyy , 2 ∂y Hence
∂ ∂ ∂v (∇2 v) + U (∇2 v) = U 00 . ∂t ∂x ∂x
(CHECK)
(3.1.11)
NB This is an equation which gives v, to find u and w we need extra equations. Consider vorticity in y−direction ∂u ∂w η = (∇ ∧ u)2 = − , ∂z ∂x can show (exercise) ∂ ∂ ∂v +U η = −U 0 . (3.1.12) ∂t ∂x ∂z 14
Hence, solve (3.1.11) to find v, v is forcing to (3.1.12) – hence η, u and w can then be obtained from continuity (3.1.10). Exercise: ∂ 2v ∂η ∂ 2u ∂ 2u − , + = ∂x2 ∂z 2 ∂z ∂x∂y and similar equation for w. NB– This is a spatial equation – no effect on stability. Return to (3.1.11) and use normal mode analysis. Write v(x, y, z, t) = v˜(y)ei(kx+lz−ckt) (i.e. ω = ck), d2 v˜ 2 2 2 ∇v = −k v˜ + 2 − l v˜ ei(kx+lz−ckt) , dy 2 d v˜ 2 = − K v˜ ei(kx+lz−ckt) , K 2 = k 2 + l2 , 2 dy in (3.1.11) −ick
2 d v˜ d2 v˜ 2 2 − K v˜ + ikU − K v˜ = U 00 ik˜ v, dy 2 dy 2 2 d v˜ 2 (U − c) − K v˜ − U 00 v˜ = 0. 2 dy
Given U (y) we have a second order ODE for v˜. Need 2 boundary conditions on v˜(y). Bounded flow – e.g. in channel −1 < y < 1. Normal velocity on wall is zero v˜(−1) = v˜(1) = 0. Unbounded flow – e.g. Kelvin–Helmholtz. Condition is that the disturbance is localised v˜ → 0 as |y| → ∞. NB If U = c when y = yc , we have a singular point (coeff. of highest derivative in ODE is zero). Return to this later in §3.3.3. ∂ For remainder of this section concentrate on solution of ODE for v˜. Writing D = ∂y (U − c)(D2 − K 2 )˜ v − U 00 v˜ = 0, known as Rayleigh’s equation.
3.2 3.2.1
Piecewise Linear Base Flows Theory
Consider flows such as
U (y) =
U0 y b
U0
Or more complicated 15
0
(3.1.13)
y
b
0
U0
U
Figure 16: Simple piecewise linear flow.
y
III y
2
II y
1
I 0
U
Figure 17: More complicated piecewise linear flow. In each region, U depends linearly on y (i.e. U = αy + β for some α, β) =⇒ U 00 (y) = 0. Subst. into Rayleigh Eqn (D2 − K 2 )˜ v = 0, in each region =⇒ v˜(y) = AeKy + Be−Ky , or = E cosh Ky + F sinh Ky, with constants A, B or E, F different in each region. Need matching condition between each region. Consider point y = ys where U (y) has different linear form, for y < ys and y > ys . U (y) and/or U 0 (y) is discontinuous at y = ys . 16
Rayleigh Eqn
d dy
(U − c)˜ v 00 − U 00 v˜ = (U − c)K 2 v˜, d˜ v 0 (U − c) − U v˜ = (U − c)K 2 v˜. dy
Now
Z
ys +
dy, ys −
and let → 0, Z ys + ys+ d˜ v 2 0 = K lim (U − c)˜ v dy , (U − c) − U v˜ →0 dy ys − ys− = 0. +
Letting [·] denote [·]yys− , s
[(U − c)˜ v 0 ] = [U 0 v˜].
Need 2nd matching condition. Let f = (U − c)˜ v 0 − U 0 v˜, f v˜0 U 0 v˜ d v˜ = − = . (U − c)2 (U − c) (U − c)2 dy U − c Again
R ys + ys −
dy and let → 0.
Z ys + f v˜ = dy → 0, 2 U −c ys − (U − c)
as → 0 (assuming U (ys ) 6= c). Hence 2 matching conditions
v˜ = 0 and U −c
[(U − c)˜ v 0 ] = [U 0 v˜].
In most examples considered U (y) is continuous, i.e. [U ] = 0, hence
3.2.2
[˜ v ] = 0, d˜ v [U 0 ]˜ v (ys ) = . dy U (ys ) − c
Examples
1–One part boundary layer U (y) =
U0 y b
U0
v˜(0) = 0 and v˜ → 0 as y → ∞. 17
0
y II b I 0
U0
U
Figure 18: One part boundary layer. In I v˜I = A eKy − e−Ky . In II v˜II = Be−Ky . Matching conditions at y = b, [˜ v] = 0 :
A eKb − e−Kb
= Be−Kb , A e2Kb − 1 = B.
v˜(b)[U 0 ] , [˜ v] = U (b) − c [˜ v 0 ] = −KBe−Kb − KA eKb + e−Kb , = −Ke−Kb B + A e2Kb + 1 , Be−Kb − Ub0 v˜(b)[U 0 ] = . U (b) − c U0 − c 0
So K B + A e2Kb + 1
=
U0 B . b(U0 − c)
Eliminating B K
U0 (e2Kb − 1) e2Kb − 1 + e2Kb + 1 = , b(U0 − c) Kb(U0 − c)2e2Kb = U0 e2Kb − 1 , U0 e2Kb − 1 U0 − c = , 2Kbe2Kb 1 − e−2Kb c = 1− . U0 2Kb 18
c is real for all values of K (i.e. for all k and l), i.e. Neutrally stable. c = f (X) U0
X = 2Kb, 1 − e−X f (X) = 1 − . X
f(X) 1
X Figure 19: f(X). So 0 < c < U0 and c is larger for larger K i.e. short waves travel fastest. 2–Mixing layer U (y) =
U0
U0 y b
−U0
y>b |y| < b , y < −b
I b −U0
U0 II −b III
Figure 20: Piecewise mixing layer. B. conditions on v˜(y) v˜ → 0 as |y| → ∞. 19
vI = P e−Ky , vII = Qe−Ky + ReKy , vIII = SeKy .
I: II : III : [˜ v ] = 0 at y = b
P = Q + Re2Kb .
(3.2.14)
S = R + Qe2Kb .
(3.2.15)
[˜ v ] = 0 at y = −b At y = b [˜ v 0 ] = −KP e−Kb − −KQe−Kb + KReKb , = −Ke−Kb P − Q + Re2Kb , P e−Kb − Ub0 v˜(b)[U 0 ] = . U (b) − c U0 − c Using [˜ v0] =
v˜[U 0 ] , U0 − c
and eliminating P using (3.2.14) Q = Re
2Kb
c −1 . 2Kb 1 − U0
Similarly at y = −b 2Kb
R = Qe
c −1 . 2Kb 1 + U0
(Check this) Eliminate Q and R 1 = e−4Kb = c2 = U02 =
c c e 2Kb 1 − −1 2Kb 1 + −1 , U0 U0 c2 2 2 4K b 1 − 2 − 4Kb + 1, U0 2 1 e−4Kb 1− − , 2Kb 4K 2 b2 2 1 e−4X 1− − = f (X), 2X 4X 2 4Kb
where X = Kb. If f (X) < 0 then c is imaginary =⇒ Unstable. If f (X) > 0 then c is real =⇒ Neutrally stable. As X → ∞, f → 1. As X → 0 1 1 1 (4X)2 f (X) = 1− + − 1 − 4X + , X 4X 2 4X 2 2 = 1 − 2, = −1. 20
f(X) 1
0.639
X
−1 Figure 21: f(X). For small X f < 0 i.e. Unstable for small K. For large X f > 0 i.e. Neutrally stable for large K. A MATLAB handout for the unbounded mixing layer can be found in appendix A.2.2. 3–Two part boundary layer 0 < y > 21 2by 2(1 − b)y + (2b − 1) 12 < y < 1 U (y) = 1 y>1 See the MATLAB worksheet in appendix A.2.3.
3.3 3.3.1
Continuous velocity profiles Criteria for inviscid instability
Consider U (y) with −1 ≤ y ≤ 1. Rayleigh’s equation (U − c)(D2 − K 2 )˜ v − U 00 v˜ = 0,
v˜(±1) = 0,
ט v ∗ (complex conjugate) v˜∗
U 00 d2 v˜ 2 2 − K |˜ v | − |˜ v |2 = 0, dy 2 U −c 1 Z 1 Z 1 2 ˜ v d˜ v d˜ v∗ ∗d v ∗ d˜ v˜ dy = v ˜ − dy, dy 2 dy −1 −1 −1 dy dy Z 1 2 d˜ v dy by b.c.. = 0− −1 dy
Hence Z
1
−1
2 ! Z 1 d˜ v U 00 2 2 K |˜ v | + dy = − |˜ v |2 dy, dy −1 U − c Z 1 (U − c∗ )U 00 2 = − |˜ v | dy. |U − c|2 −1 21
(3.3.16)
Taking imaginary part Z
1
0=− −1
ci U 00 |˜ v |2 dy, |U − c|2
so for unstable wave, ci > 0, =⇒ Z
1
−1
U 00 |˜ v |2 dy = 0, |U − c|2
(3.3.17)
Therefore U 00 (y) must change sign in interval −1 ≤ y ≤ 1. This is known as RAYLEIGH’S CONDITION, a necessary condition for inviscid instability is that the basic flow has an inflection point. i.e. ys such that U 00 (ys ) = 0. Taking the real part of (3.3.16) Z 1 (U − cr )U 00 2 |˜ v | dy < 0. (3.3.18) |U − c|2 −1 If U 00 (ys ) = 0 and U (ys ) = Us 1
Z (Us − cr )
−1
U 00 |˜ v |2 dy = 0 by (3.3.17). |U − c|2
(3.3.19)
Hence (3.3.18)-(3.3.19) =⇒ Z
1
−1
(U − Us )U 00 2 |˜ v | dy < 0. |U − c|2
Hence U 00 (y) (U (y) − U (ys )) , must be negative somewhere in the flow. This is known as FJORTOFT’S CRITERION- another necessary condition for inviscid instability. Examples 1. U (y) = tanh y, (either bounded, i.e. −1 < y < 1 or unbounded) U 0 = sech2 y,
U 00 = −2sech2 y tanh y.
Hence U 00 (ys ) = 0 =⇒ ys = 0 so U (ys ) = 0, So satisfies Rayleigh’s criterion. U 00 (U − U (ys )) = −2sech2 y tanh y tanh y ≤ 0 (equality only when y = 0). So satisfies Fjortoft’s criterion, =⇒ likely to be unstable for some disturbances. (Not guaranteed because conditions are not sufficient).
22
2. y − 1 < y < 1, 4 − y2 2y(y 2 + 12) = . (4 − y 2 )3
U (y) = U 00
U 00 (0) = 0 =⇒ satisfies Rayleigh criterion. U 00 (y)(U (y) − U (0)) =
2y 2 (y 2 + 12) ≥ 0, (4 − y 2 )4
Does not satisfy Fjortoft’s criterion =⇒ Flow is stable. (stable to inviscid disturbances). 3. Blasius Boundary layer (see later in course) f 000 + f f 00 = 0, with boundary conditions f (0) = f 0 (0) = 0,
f 0 (y) → 1 as y → ∞,
U (y) = f 0 (y). Solve numerically U 00 (y) < 0 for 0 < y < ∞, Does not satisfy Rayleigh criterion =⇒ Stable to inviscid disturbances. This can be seen in worksheet A.2.4. 3.3.2
Numerical calculation of eigenvalues c(K) for continuous U (y)
Example U (y)
− 1 < y < 1,
Rayleigh’s equation (U − c)(D2 − K 2 )˜ v − U 00 v˜ = 0
v˜(±1) = 0.
1. Fix K and choose a value of of c. 2. Fix v˜(−1) = 0 and d˜ v /dy(−1) = 1 (Non-zero value, doesn’t matter on exact value, linear equation). Hence have 2 initial conditions for 2nd order ODE. Solve numerically, i.e. advance y from −1 to 1. 3. Hence find v˜(1). If v˜(1) = 0, then we have chosen the correct value of c. If v˜(1) 6= 0, need to change value of c and re–solve equation. Tricky part is to choose new guess fore c, particularly since c is complex. For numerical solution for U = tanh y, see worksheet A.2.5. 23
3.3.3
Critical layers in inviscid flows
Rayleigh’s eqn (U − c)(D2 − K 2 )˜ v − U 00 v˜ = 0. Singular point when coefficient of highest derivative is zero. If U1 < U < U2 then for some c real and U1 < c < U2 we have some point y = yc where U (yc ) = c. y = yc is known as the critical layer.
y
yc
U1
c
U2
U
Figure 22: Critical layer. At the critical layer y = yc , U (yc ) − c = 0. Approximate Rayleigh’s equation for y ≈ yc . Let Y = y − yc . 1 U ∼ U (yc ) + Y U 0 (yc ) + Y 2 U 00 (yc ) + O(Y 3 ), 2 1 2 00 0 ∼ c + Uc Y + Y Uc + O(Y 3 ) where Uc0 = U 0 (yc ), 2 U 00 ∼ Uc00 + Y Uc000 + O(Y 2 ).
2 1 2 00 d 2 (Y + Y Uc ) − K v˜ − (Uc00 + Y Uc000 )˜ v 2 dY 2 2 d Uc000 1 Uc00 0 2 00 Y Uc 1 + Y 0 − K v˜ − Uc 1 + Y 00 v˜ 2 Uc dY 2 Uc 2 00 00 −1 d Uc 1 Uc Uc000 2 Y − K v˜ − 0 1 + Y 0 1 + Y 00 v˜ dY 2 Uc 2 Uc U !c! 2 d2 v˜ Uc00 Uc000 1 Uc00 2 Y − + K + − Y v˜ dY 2 Uc0 Uc0 2 Uc0 Uc0
using binomial expansion of 1 +
1 Uc00 Y U0 2 c
Y
−1
= 0, = 0, = 0, = 0,
.
d2 v˜ − (p + qY )˜ v = 0, dY 2 24
(3.3.20)
where p, q constants given by U . Equation valid for small Y . Can extend to higher powers of Y . Look for series solutions of (3.3.20) v˜ = Y
λ
∞ X
an Y n
a0 6= 0.
n=0
Y
∞ X
(n + λ)(n + λ − 1)an Y
n=0 ∞ X
λ(λ − 1)a0 Y −1 + Y
n+λ−2
−
∞ X
pan Y
n=0 ∞ X
(m + 1 + λ)(m + λ)am+1 Y m −
m=0
n+λ
−
∞ X
qan Y n+λ+1 = 0,
n=0 ∞ X
pam Y m −
m=0
qam−1 Y m = 0.
m=0
Equating coeffs of Y −1 λ(λ − 1)a0 = 0, =⇒ λ = 0 or λ = 1 (∵ a0 6= 0). Equating coeffs of Y 0 λ(λ + 1)a1 − pa0 = 0 a1 =
p a0 . λ(λ + 1)
Eq coeffs of Y 1 , obtain a2 in terms of a0 and a1 ect. Taking λ = 1 1 a1 = pa0 , 2 1 v˜ = Y (1 + pY ) + O(Y 3 ) 2 This is one solution of (3.3.20). Choosing λ = 0 =⇒ a1 = ∞. Can’t generate second solution in this way. To find second solution (Frobenius method). Let v˜ = v˜1 φ v˜1 is first solution. Subst into equation (3.3.20) v˜1 φ00 + 2˜ v10 φ0 = 0. Let f = φ0 , solve 1st order ODE f (Y ) = k˜ v1−2
(can set k = −1 to generate 2nd soln)
1 1 = − 2 2 (1 − pY + O(Y 2 )), 2 v˜ a0 Y 1 p φ = 2 + 2 ln Y + O(Y ). a0 Y a0
dφ dY
= −
Hence second solution 1 v˜2 = v˜1 φ = a0
1 pY 1 1 + pY + ln Y 1 + pY + · · · . 2 a0 2 25
NB: Appearance of ln term needs special treatment, see later. General solution is now v˜ = A˜ v1 + B˜ v2 . General Theory U (yc ) = c and assume U 0 (yc ) 6= 0. 2 linearly independent solutions. Tollmien inviscid solutions v˜1 (y) = (y − yc )P1 (y), U 00 v˜2 (y) = P2 (y) + c0 v˜1 (y) ln(y − yc ) Uc where P1 (y) and P2 (y) are analytic at y = yc and P1 (yc ) = P2 (yc ) = 1. Check this agrees with first two terms of solutions derived above. General solution v˜ = A˜ v1 (y) + B˜ v2 (y). If flow is bounded, i.e. y1 < y < y2 with v˜(y1 ) = v˜(y2 ) = 0, then v˜ = A˜ (˜ v1 (y)˜ v2 (y1 ) − v˜2 (y)˜ v1 (y1 )) , applying b.c at y = y1 . Applying b.c at y = y2 , for non–zero v˜ (i.e. A˜ = 6 0), need v˜1 (y2 ) v˜1 (y1 ) = . v˜2 (y1 ) v˜2 (y2 ) Each of v˜1 and v˜2 depends on c and K. So this gives us our dispersion relation. i.e. c(K). Branch cuts–brief revision ln Z is a multivalued function. Set Z = Rei(θ+2nπ) with |Z| = R, arg(Z) = θ, where n ∈ Z. ln Z = ln R + i(θ + 2nπ). Define a branch cut which originates from branch point at Z = 0, so that ln Z is uniquely defined and continuous everywhere except for a discontinuity across the branch cut. Define α − 2π < arg(Z) < α ln Z = ln R + iarg(Z). Clearly got discontinuity across the branch cut. [ln Z] = 2πi. Consider 2nd Tollmien inviscid solution. Branch point at y = yc , U (yc ) = c. If c is complex then yc is typically complex (e.g. U (y) = tanh y, c = λi =⇒ yc = tan−1 (λ)i). If c is real, may have branch point on the real axis. e.g. U (y) = 1 − y 2 26
− 1 < y < 1,
cut Z α Z=0 Figure 23: Branch cut.
yc yc
(b)
(a)
Figure 24: Branch points on real axis. If
c = cr
yc =
√
1 − cr .
Branch points on real axis for 0 < cr < 1. If yc is not on the real axis, define branch cut so that it doesn’t cross the real axis, i.e. parallel to the imaginary axis, away from the real axis. If branch point lies on real axis, have 2 choices (a) For y < yc ln(y − yc ) = ln |y − yc | − iπ. For y > yc ln(y − yc ) = ln |y − yc |. (b) For y < yc ln(y − yc ) = ln |y − yc | + iπ. For y > yc ln(y − yc ) = ln |y − yc |. Notes 1. Which branch cut to take is discussed in §4, effect of viscosity. 2. v2 (y) is continuous at y = yc =⇒ normal velocity v˜ is continuous.
27
3. As y → yc v2 ∼ 1 +
Uc00 (y − yc ) ln(y − yc ). Uc0 00
Discontinuity in v20 at y = yc , size= ±iπ UUc0 , c
=⇒ Discontinuity in
d˜ v , dy
=⇒ Jump in u˜(y), i.e. streamwise velocity across the critical layer, known as a phase jump.
28
4 4.1
Viscous stability Introduction
Non–dimensional Navier–Stokes 1 ∂uT + uT · ∇uT = −∇pT + ∇2 uT , ∂t R here R is the Reynolds number. 4.1.1
Boundary conditions for viscous flow
If solid boundary has velocity Us , then fluid velocity on surface is uT = U s , cf inviscid theory uT · n = Us · n, i.e. condition on normal velocity at boundary – no condition on tangential velocity. e.g. Fixed boundary at y = 0: viscous conditions uT = vT = wT = 0 on y = 0, inviscid conditions vT = 0 on y = 0. 4.1.2
Strictly parallel base flows uT = U (y)i + u(x, t), pT = P + p.
Subst. into N–S,
∂u ∂u 1 + U + u · ∇(U i) = −∇P − ∇p + (∇2 U )i + ∇2 u + O(2 ). ∂t ∂x R R
Equate coefficients of 0 : −∇P +
1 d2 U i = 0. R dy 2
Hence from j and k components ∂P ∂P = 0 and = 0, ∂y ∂z =⇒ P = P (x), i component dP 1 d2 U = , dx R dy 2 Func of x only Func of y only 29
(4.1.21)
Hence both sides must be constant. =⇒
∂P = constant, ∂x
(constant pressure gradient) and U (y) must be quadratic U (y) = ay 2 + by + c, (compare with inviscid theory where any base flow could be considered). Special cases 1. U (y) = 1 − y 2
−1 ≤ y ≤ 1 Poiseuille flow dP = constant 6= 0, dx
y=1
y=−1
Figure 25: Poiseuille flow 2. U (y) = y
−1 ≤ y ≤ 1 Plane Couette flow dP = 0, dx
i.e. P = constant.
i.e. gives viscous flow between 2 plates moving in opposite directions. (NB In general Couette flow describes flow between 2 cylinders rotating in opposite directions).
4.2
Orr–Sommerfeld Equation
Equation governing unsteady disturbances. From (4.1.21), equating O() coefficients, ∂u ∂u dU 1 +U +v i = −∇P + ∇2 u. ∂t ∂x dy R Now focus only on normal modes i(kx+lz−ckt) ˜ ei(kx+lz−ckt) = (˜ i.e. u = u u(y), v˜(y), w(y))e ˜ p = p˜(y)ei(kx+lz−ckt) ,
30
y=1
y=−1
Figure 26: Plane Couette flow 1 dU d˜ p u˜yy − K 2 u˜, v˜yy − K 2 v˜, w˜yy − K 2 w˜ . ik(U − c)˜ u + v˜ i = − ik p˜, , ilp˜ + dy dy R In component form 1 u˜yy − K 2 u˜ , R 1 U0 (U − c)˜ u− u˜yy − K 2 u˜ = −˜ p − v˜, ikR ik
ik(U − c)˜ u + U 0 v˜ = −ik p˜ +
or 2 D − K 2 − (ikR)(U − c) u˜ = ikR˜ p + RU 0 v˜, 2 D − K 2 − (ikR)(U − c) w˜ = ilR˜ p, 2 2 D − K − (ikR)(U − c) v˜ = RDp˜.
(4.2.22) (4.2.23) (4.2.24)
Also have continuity equation ∇ · u = 0, ik˜ u + D˜ v + ilw˜ = 0. Also have b. condns. e.g. u˜ = v˜ = w˜ = 0, on fixed boundaries e.g. y = y1 and y2 . 4.2.1
Squire’s Transformation
Using this transformation can simplify the system K uˆ = k˜ u + lw, ˜ K pˆ = p˜, k vˆ = v˜, ˆ = k R. R K ˆ New set of variables uˆ, vˆ, pˆ, R. 31
(4.2.25)
Notation : L = D2 − K 2 − (ikR)(U − c), ˆ = D2 − K 2 − (iK R)(U − c). k × (4.2.22) + l × (4.2.23) KLˆ u = i(k 2 + l2 )R˜ p + kRU 0 v˜, k = iK 2 R pˆ + kRU 0 vˆ, K ˆ p + K RU ˆ 0 vˆ = iK 2 Rˆ ˆ p + RU ˆ 0 vˆ. Lˆ u = (iK R)ˆ From (4.2.24) ˆ pˆ, Lˆ v = RD and from (4.2.25) iK uˆ + Dˆ v = 0. Eliminate pˆ 0 ˆ pˆ + RD(U ˆ D(Lˆ u) = iK RD vˆ), ˆ R 0 ˆ = iK Lˆ v + RD(U vˆ). ˆ R
Next eliminate uˆ 0 ˆ D(L(iK uˆ)) = −K 2 Lˆ v + iK RD(U vˆ), 2 0 ˆ −D(L(Dˆ v )) = −K Lˆ v + (iK R)D(U vˆ).
LHS: ˆ −D(L(Dˆ v )) = −D(D2 − K 2 − (iK R)(U − c))Dˆ v, 2 2 2 ˆ ˆ 0 Dˆ = −D (D − K )ˆ v + (iK R)(U − c)D2 vˆ + (iK R)U v. RHS: 00 ˆ ˆ −K 2 (D2 − K 2 )ˆ v + K 2 (iK R)(U − c)ˆ v + (iK R)(U vˆ + U 0 Dˆ v ).
LHS=RHS ˆ (U − c)D2 vˆ + U 0 Dˆ (iK R) v − K 2 (U − c)ˆ v − U 00 vˆ − U 0 Dˆ v = (D2 − K 2 )2 vˆ 1 (U − c)(D2 − K 2 )ˆ v − U 00 vˆ = (D2 − K 2 )2 vˆ ˆ iK R ORR–SOMMERFELD EQUATION. Notes ˆ → ∞ (i.e. viscosity vanishes) we obtain Rayleigh’s equation. 1. As R
32
2. (D2 − K 2 )2 vˆ = D4 vˆ − 2K 2 D2 vˆ + K 2 vˆ. 4th order differential. i.e. 4th order ODE for vˆ(y) – cf 2nd order equation (Rayleigh’s) for inviscid flows. so need 2 extra boundary conditions. For bounded flows (i.e. y1 < y < y2 ) inviscid conds: vˆ(y1 ) = vˆ(y2 ) = 0, viscous flow u(y1 ) = u(y2 ) = 0, i.e. vˆ(y1 ) = vˆ(y2 ) = uˆ(y1 ) = uˆ(y2 ) = 0. But iK uˆ = −Dˆ v, dˆ v dˆ v (y1 ) = (y2 ) = 0. dy dy
i.e. 3.
k R K ˆ ≤ R. R ˆ = R
K = (k 2 + l2 )1/2
To each 3 dimensional solution to Orr–Sommerfeld equation, there corresponds a 2D solution at a lower Reynolds number. (cf Results for inviscid – 2D modes (l = 0) grow fastest).
4.3
Squire’s Equation
As for inviscid case, need another equation involving u˜(y) and w(y) ˜ to fully determine velocity. As before define ∂u ∂w η= − ∂z ∂x ((∇ ∧ u)2 - i.e. vorticity in y direction). Normal modes η˜ = il˜ u − ik w. ˜ Recall (4.2.22) (4.2.23)
L(˜ u) = ikR˜ p + RU 0 v˜, L(w) ˜ = ilR˜ p,
il(4.2.22) − ik(4.2.23) L(˜ η ) = ilRU 0 v˜ Kˆ 0 = il RU v˜. k This is known as SQUIRE’S EQUATION. 33
(4.3.26)
Forced by v˜ – solution of Orr–Sommerfeld eqn (OSE), but may also have eigensolutions. i.e. Consider solution with v˜ = 0 and L(˜ η ) = 0, i.e.
ˆ − c) η˜ = 0. (D2 − K 2 ) − iK R(U
NB ˆ The eigensolutions to the Squire equation are damped, i.e. Im(c) < 0 for all K, R. Proof Multiply by η˜∗ and
4.4
R y2 y1
dy.
Formal theory of eigenvalues of OSE
Could try the same method as for Rayleigh eqn to find e-values c. Assume bounded flow y1 < y < y2 . 1. vˆ(y1 ) =
dˆ v (y1 ) = 0, dy
ˆ fix K and R. 2. For numerical solution, need 2 extra conds at y1 . Fix d2 vˆ (y1 ) = 1. dy 2 Choose a value for
d3 vˆ (y ) dy 3 1
and a value for c.
dˆ v 3. Integrate eqn (i.e. solve numerically starting from y = y1 ) If vˆ(y2 ) = dy (y2 ) = 0, then have correct c. In practice enormous numerical problems, no one uses this method.
Alternative method 4th order eqn =⇒ 4 lin indep solns vˆ1 , vˆ2 , vˆ3 , vˆ4 , vˆ(y) = A1 vˆ1 + A2 vˆ2 + A3 vˆ3 + A4 vˆ4 . 4 conditions vˆ(y1 ) = vˆ(y2 ) = and 4 unknowns A1 , A2 , A3 , A4 . Solve in matrix form vˆ1 (y1 ) vˆ1 (y2 ) 0 vˆ1 (y1 ) vˆ10 (y2 )
vˆ2 (y1 ) vˆ2 (y2 ) vˆ20 (y1 ) vˆ20 (y2 )
dˆ v dˆ v (y1 ) = (y2 ) = 0, dy dy
vˆ3 (y1 ) vˆ3 (y2 ) vˆ30 (y1 ) vˆ30 (y2 ) 34
vˆ2 (y1 ) A1 vˆ2 (y2 ) A2 vˆ20 (y1 ) A3 vˆ20 (y2 ) A4
= 04 .
Non–zero solutions if determinate of matrix is zero. Hence, formally ˆ c) = 0, F (K, R, (problem is finding indep solns). ˆ we can find a set of solutions for c Hence given K and R (n)
c(n) = c(n) r + ici . NB For bounded flows, there is an infinite set of discrete eigenvalues. 4.4.1
Finding eigenvalues numerically using Chebyshev polynomials
We can use Chebyshev polynomials to approximate solutions to the OSE (and Rayleigh equation), and acquire the ‘global’ set of eigenvalues for a particular base flow, U (y). Chebyshev polynomials are an orthogonal set of functions defined by Tn (x) = cos(n cos−1 x), where n is an integer which gives the order of the polynomial. The first five Chebyshev polynomials are T0 T1 T2 T3 T4
= = = = =
1, x, 2x2 − 1, 4x3 − 3x, 8x4 − 8x2 + 1.
We can approximate a function f (x) defined on [−1, 1] by f¯(x) =
N −1 X
an Tn (x),
n=0
where an are the coefficients of each polynomial that need to be found. One advantage of Chebyshev polynomials is |Tn (x)| ≤ 1 for all x ∈ [−1, 1], so each term in f¯(x) is bounded. Rather than work out the coefficients of each Chebyshev polynomials in our expansion, we can find the numerical value of each Chebyshev polynomial some point at x = X via the recurrence relation Tn+1 (X) = 2XTn (X) − Tn−1 (X), by noting that T0 (X) = 1 and T1 (X) = X. Then the numerical value of the k th derivative of each polynomial at x = X can be found via the recurrence relation (k)
T0 (X) = 0, (k)
(k−1)
T1 (X) = T0 (k) T2 (X)
=
Tn(k) (X) =
(X),
(k−1) 4T1 (X), (k−1) 2nTn−1 (X)
+
n (k) Tn−2 (X) for n ≥ 3. n−2
35
Therefore we have to map our solution domain to y ∈ [−1, 1], and substitute vˆ =
N −1 X
an Tn (y)
(4.4.27)
n=0
into the Orr-Sommerfeld equation. Evaluating this function at the special points jπ yj = cos N −1 (known as collocation points) for j = 2, · · · , N −3 will produce N −4 equations for the N unknowns an . Also, by using the 4 boundary conditions vˆ(±1) = vˆ0 (±1) = 0 and noting that Tm (±1) = (±1)m
Tm0 (±1) = (±1)m−1 m2 ,
(CHECK)
gives us a system of N equations and N unknowns. This system can be written as Aa = cBa where c is the eigenvalue and a is the vector of the unknown coefficients an . This equation can then be solved using a matrix eigenvalue solver (see books on matrix methods or the MATLAB help for more information). The resulting numerical solution will produce N eigenvalues and for each eigenvalue it will produce a vector a of coefficients which can then be used to produce the eigenfunction from (4.4.27). The more Chebyshev polynomials used in the approximation, the more eigenvalues will be found, however there are numerical problems in using this method. An example of this approach for solving the Poiseuille flow problem is given in worksheet A.3.1. 4.4.2
Summary of numerical results for OSE
ˆ sufficiently small, all c(n) < 0 so flow is stable. 1. For R i ˆ increases, one eigenvalue crosses the real axis, 2. As R (1)
ci > 0,
i.e. and hence we have an unstable mode.
ˆ 3. Hence we obtain a region in (K, R)–plane where flow is unstable.
4.5
Viscous stability analysis for Poiseuille flow U (y) = 1 − y 2 ,
− 1 ≤ y ≤ 1,
00
U (y) = −2. Does not satisfy Rayleigh’s criterion (i.e. U 00 = 0 somewhere in the flow). Hence it is NOT inviscidly unstable. ˆ We will see that the flow is unstable for ranges of K and R. OSE: 1 (D2 − K 2 )2 vˆ, (U − c)(D2 − K 2 )ˆ v − U 00 vˆ = ˆ iK R ˆ and can consider problem as 2D (use Squire’s transform to convert back to 3D). will dropˆoff R, 36
4.5.1
Numerical results
See handouts from books for calculations of eigenvalues for R = 10, 000. Figure 3.1a Schmid and Henningson – 3 classes of eigenvalues A, P, S. Split modes into odd vˆ(−y) = −ˆ v (y), even vˆ(−y) = vˆ(y). See also Drazin and Reid, fig 4.19. This just shows even modes. Notes 1. Infinite number of e-values on the S branch, finite on A and P. 2. Average mean flow between y = y1 and y2 Z y2 Z 1 1 1 2 U (y) dy = = 1 − y 2 dy = , y2 − y1 y1 2 −1 3 same as cr (phase speed) of S modes. 3. For R > Rc (K), one e-value crosses real axis but only one mode. This is an even mode, subsequent analysis focus on even modes. 4. The curve R = Rc (K) is known as the Neutral Stability Curve. See fig 3.8 – Schmid and Henningson. Minimum Rc (K) as K varies is known as the critical Reynolds number. Rcrit = min(Rc (K)) as K varies, = 5767 corresponding to K ≈ 1.02. Typical neutral curves
K Upper branch UNSTABLE lower branch R
R crit Figure 27: Typical neutral curve If R < Rcrit – system is stable to all disturbances.
37
4.5.2
Asymptotic theory for ODEs
Consider the ODE
dy + y = x, (4.5.28) dx with x ≥ 0 and 0 < 1. One b. cond y(0) = 1. Solve naively by assuming = 0, and we have y = x. But this contradicts the boundary condition!!!!!!! dy The paradox arises from dx not becoming small as → 0, this term is as important as the other two.
y y=x
1
x Figure 28: The approximate solution y = x is plausible, provided neither x nor y is small. Where x is small we expect a boundary layer to join the line y = x to the boundary condition x = 0, y = 1. Rescale x = Xα where X = O(1) and α =constant. (NB Because ODE is linear in y we don’t have to rescale w.r.t y). dy + y = x, dx dy + y = Xα , α dX dy + α−1 y = 2α−1 X. dX (A) (B) (C)
1. Does (A) balance (B)? dy + α−1 y = 0 =⇒ α = 1. dX 1 Therefore, term (C) is → 0 as → 0, so consistent. 2. Does (A) balance (C)?
38
dy 1 = 2α−1 X =⇒ α = . dX 2 −1/2 But (B) is now y → ∞ as → 0, so inconsistent. 3. Does (B) balance (C)? α−1 y = 2α−1 X, y = α X, Therefore α = 0 and y = X = x, solution we already had. Returning to (1) α = 1, x = X. dy + y = X, dX The solution of this which satisfies y(0) = 1 is called the inner solution. Let yinner (X) = y0 (X) + y1 (X) + O(2 ), then ODE gives (y00 + y10 ) + y0 + y1 + O(2 ) = X, (y00 + y0 ) + (y10 + y1 − X) + O(2 ) = 0, where prime denotes diff w.r.t X. Now compare different orders of . O(1) :
y00 + y0 Z dy0 y0 =⇒ y0 y0
= 0
satisying
y(0) = 1,
Z = −
dX,
= Ae−X = e−X .
A = 1 to satisfy b.c.,
At O(), y10 + y1 − X = 0, d y1 eX = XeX , dX Z −X XeX dX, y1 = e Z −X X X y1 (X) = e Xe − e dX , y1 (X) = X − 1 + Be−X . the boundary condition gives y(0) = 1, y0 (0) + y1 (0) + O(2 ) = 1, ∴ y0 (0) = 1 and y1 (0) = 0 (no 0 s on RHS). 39
This gives us y1 (0) = 0 − 1 + B = 0, B = 1, y1 (X) = X − 1 + e−X . Recap: yinner (X) = e−X + X − 1 + e−X + O(2 ), inner solution for which X = O(1), i.e. x = O(). The outer solution is y = x which you can check formally by substituting youter (x) = Y0 (x) + Y1 (x) + O(2 ), in the outer region x . The solution consists of two expressions, each valid in a different region. We hope the expressions match to one another in an overlap region, known as the matching region. In this region we have two limiting processes to consider, 1. x → 0 in outer solution, 2. X → ∞ with → 0 and x remaining fixed (X = x−1 ), for the inner solution. Limit process (2) applied to our inner solution yinner = e−X + e−X + X − 1 , x = e−x/ + e−x/ + − 1 . In letting → 0 with x fixed, both terms e−x/ are exponentially small w.r.t → 0. This leaves for small x − 1 = x − , yinner ∼ so yinner ∼ x = youter . This means that the 2 solutions in the 2 regions are all that’s needed, there is no third region. The above analysis can also be recreated by considering the exact solution to (4.5.28) which is x
y = x − + (1 + )e− . See worksheet A.3.2 for exact solution and inner and outer approximations. The above problem motivates thoughts on boundary layers in other fluid problems. 4.5.3
Asymptotic theory of Orr–Sommerfeld equation for large R
OSE with boundary conditions vˆ(−1) = Dˆ v (−1) = vˆ(1) = Dˆ v (1) = 0. Inviscid Solutions
40
If KR 1, we can then write 1 (1) vˆ . RK Substitute into OSE, =⇒ Rayleigh’s equation for vˆ(0) . vˆ = vˆ(0) +
(4.5.29)
(U − c)(D2 − K 2 )ˆ v (0) − U 00 vˆ(0) = 0, 2nd order ODE for vˆ(0) (y) – can’t satisfy 4 b. conds given. (typical of singular perturbations, where the small parameter is multiplying the highest derivative, cf §4.5.2) Solving this equation we obtain the 2 Tollmien inviscid solutions (see §3) vˆ1 = (y − yc )P1 (y), U 00 vˆ2 = P2 (y) + c0 vˆ1 (y) ln(y − yc ). Uc Writing vˆ(0) = A1 vˆ1 + A2 vˆ2 , the expansion (4.5.29) is going to become invalid when 1. y is close to yc because vˆ2 is not analytic in neighborhood of yc due to log term. Alternative approach is that U − c is small close to yc and so different terms balance in full OSE. 2. y is close to ±1 because then vˆ(0) can’t satisfy the 2 conditions at y = −1 and 2 conditions at y = 1. (cf inviscid analysis when only needed to satisfy vˆ = 0 at y = ±1) Hence need a boundary layer at y = −1 and 1 in which viscosity is important. NB Since unstable mode is even, concentrate only on even modes. Consider OSE for −1 < y < 0. Boundary conds vˆ(−1) = Dˆ v (−1) = 0, Dˆ v (0) = D3 vˆ(0) = 0. Hence no need to consider the boundary layer at y = 1. Viscous solutions Look for asymptotic forms of other 2 linearly independent solutions of OSE, vˆ3 and vˆ4 . Wall layer at y = −1 Let y = −1 + Y , 1. U (y) = 1 − y 2 = 2Y + O(2 ), U 00 (y) = −2. 41
Subs into OSE (2Y − c)
2 1 d2 1 d2 1 2 2 − K vˆ + 2ˆ v= −K vˆ. 2 dY 2 iKR 2 dY 2 dˆ v dY
To satisfy the b. conds vˆ(Y = 0) = terms from RHS. Largest term on LHS
(Y = 0) = 0, we need a dominant balance including viscous
−
c d2 vˆ . 2 dY 2
Largest on RHS 1 1 d4 vˆ . iKR 4 dY 4 Hence choose = (KR)−1/2 1 (∵ KR 1). Leading order equation becomes d2 vˆ d4 vˆ + ic = 0. dY 4 dY 2 Two solutions vˆ = 1 and vˆ = Y . These are just approximations to the inviscid solutions close to the wall. Other solutions are vˆ = exp(±ωY ) where ω 2 = −ic. Want solution which decays away from the boundary, so we can ‘match’ the solution with main part of the fluid. vˆBL = exp(−ωY ), ω = (−ic)1/2 = e−iπ/4 c1/2 for real c, = exp(−(y + 1)(KRc)1/2 e−iπ/4 ). Outer approximation to vˆ3 and vˆ4 Let
Z vˆ = exp
y
g(y) dy
for arbitrary constant σ.
σ
This is a standard technique for solving differential equations which involve rapid oscillations, known as the WKB approximation (or WKBJ). Substitute this expression into OSE to give a nonlinear equations for g(y) Z D(ˆ v ) = g exp g dy , Z 2 0 2 D (ˆ v ) = (g + g ) exp g dy . So 1 g 4 + 6g 2 g 0 + 4gg 00 + 3g 02 + g 000 − 2K 2 (g 2 + g 0 ) + K 4 = (U − c)(g 2 + g 0 − K 2 ) − U 00 . (4.5.30) iKR
42
As before balance terms on LHS (viscous) and RHS. Need g large for balance i.e. i.e.
g4 = O(g 2 ), KR g = O((KR)1/2 ).
Hence, set g = (iKR)1/2 g0 (y) + g1 (y) + (iKR)−1/2 g2 (y) + · · · , substitute into (4.5.30) and equate powers of (KR)1/2 . Obtain g02 = U − c, g0 (y) = ±(U (y) − c))1/2 , 5 U0 g1 (y) = − . 4U −c Hence Z
vˆ = exp
g dy , Z 1/2 = exp (iKR)
y
Z
g0 (y) dy exp
y
g1 (y) dy
1 + O (KR)1/2
.
Z 5 y U0 5 g1 (y) dy = − dy = − ln(U − c). 4 U −c 4 Z y −5/4 1/2 1/2 vˆ(y) = (U − c) exp ±(iKR) (U − c) dy 1 + O (KR)1/2 . Z
y
yc
Hence vˆ3 = (U − c)−5/4 exp −(KR)1/2 Q(y) , vˆ4 = (U − c)−5/4 exp +(KR)1/2 Q(y) , where iπ/4
Z
y
Q(y) = e
(U − c)1/2 dy.
yc
(NB: Defined so arg (U − c)1/2 = −π/2 for U − c < 0). If we let y → −1, then vˆ3 → vˆBL , the viscous wall solution we obtained before (CHECK). It’s clear that these solutions are not valid in the critical layer where U − c → 0. Inner approximation to vˆ3 and vˆ4 – valid in critical layer We are interested in the neutral curve, i.e. Rc (K), on which Im(c) = 0. On neutral curve, c is real and y = yc lies on the real y axis. Let y − yc = ˆZ, Z is new inner variable. U (y) = U (yc ) + (y − yc )U 0 (yc ) + · · · , = c + ˆZUc0 + O(ˆ2 ). 43
Substitute into OSE ˆZUc0
2 1 1 d2 1 d2 00 2 2 vˆ. − K vˆ − Uc vˆ = −K ˆ2 dZ 2 iKR ˆ2 dZ 2
Balance largest terms on LHS ZUc0 d2 , ˆ dZ 2 and RHS
1 1 d4 . iKR ˆ4 dZ 4
Choose ˆ = (KR)−1/3 . To leading order
d2 − iUc0 Z 2 dZ
d2 vˆ = 0, dZ 2
and this is a 4th order equation – 4 linearly independent solutions. vˆ = 1 vˆ = Z the Tollmien inviscid solutions. To find other 2 solutions, let Ω =
d2 vˆ . dZ 2
local approximation to v˜2 , local approximation to v˜1 , Find a transform from Z to Zˆ to convert equation to d2 Ω ˆ = 0, − ZΩ 2 ˆ dZ
known as Airy’s equation – commonly appears in stability analysis. ˆ and Bi(Z). ˆ See methods textbooks and MATLAB (or MAPLE) help. Two solutions Ai(Z) Summary of asymptotic structure Viscous wall layer at y = −1, width (KR)−1/2 . Critical layer – width (KR)−1/3 . NB wider than wall layer. Poiseuille flow U = 1 − y2, √ yc = ± 1 − c for c ≤ 1. So unless c is small, the critical layer is away from the wall. 4.5.4
Summary of solutions of the Orr–Sommerfeld equation
1. Two inviscid solutions vˆ1 valid as leading order approximation to viscous solution everywhere. vˆ2 valid as leading order approximation, except in critical layer. Can find correction to vˆ2 in critical layer, but doesn’t enter calculation of neutral stability. 44
2. Two viscous modes vˆ3 and vˆ4 vˆ3 exponentially large in −1 ≤ y < yc , vˆ4 exponentially small in −1 ≤ y < yc , from WKB approximations. Also, vˆ3 exponentially small in yc < y ≤ 0, vˆ4 exponentially large in yc < y ≤ 0. Expression for vˆ3 and vˆ4 from WKB valid everywhere except in critical layer. Separate expressions for vˆ3 and vˆ4 in critical layer (using matched asymptotic expansions). Close to boundary y = −1, expression for vˆ3 as a viscous boundary layer, comes from limit of WKB form as y → −1. 3. Structure of solution Recall neutral stability curve
K Upper branch UNSTABLE lower branch R
R crit Figure 29: Typical neutral curve
On neutral stab. curve KR 1, and c = cr + ici , ci = 0 (∵ Neutral stability) and cr 1 (not proved here, but follows from later). For small c, position of critical layer? U (y) = 1 − y 2 ,
− 1 ≤ y < 0.
Critical layer at yc where U (yc ) = c, 1 − yc2 = c, √ yc = ± 1 − c. So for cr 1, critical layer is close to boundary y = −1. 45
Size of c determines whether the critical layer overlaps the narrower viscous wall layer. Lower branch
c = O (KR)−1/3 , =⇒ yc + 1 = O (KR)−1/3 . y
y c =O(KR) (KR)
−1/3
−1/2
x
Figure 30: Viscous wall layer inside critical layer i.e. viscous wall layer within critical layer. Often called Triple Deck theory.
3
outer
2
critical
1
viscous+critical
Figure 31: Triple deck structure Upper branch
c = O (KR)−1/5 , yc + 1 = O (KR)−1/5 . viscous wall layer lies outside critical layer. Referred to as Five–Deck Theory.
46
y
yc =O(KR) (KR)
−1/3
−1/5
O(KR)
−1/2
x
Figure 32: Viscous wall layer outside critical layer
5 3/4 2 1
Figure 33: Five–deck structure 4.5.5
Calculating neutral stability curve
vˆ3 gives viscous boundary layer at y = −1, vˆ4 gives viscous boundary layer at y = 0. But for even modes there is no boundary layer at y = 0. vˆ4 can be dropped from the general solution for vˆ in −1 ≤ y < 0 (can formally be justified with error calculation). Write general solution as vˆ = Aˆ v1 + vˆ2 + C vˆ3 , To satisfy boundary condition at y = 0, Dˆ v = D3 vˆ = 0 (∵ even mode), Aˆ v10 (0) + vˆ20 (0) + C vˆ30 (0) = 0. vˆ3 is exponentially small at y = 0
Check that Dˆ v=0 Write
=⇒
vˆ20 (0) =⇒ A = − 0 . vˆ1 (0)
D3 vˆ = 0 from OSE using U (y) = 1 − y 2 .
vˆ = Φ(y) + C vˆ3 (y), vˆ0 (0) Φ(y) = vˆ2 (y) − 20 vˆ1 (y). vˆ1 (0)
47
Satisfy boundary condition at y = −1 vˆ(−1) = 0, Dˆ v (−1) = 0,
=⇒ =⇒
Φ(−1) = −C vˆ3 (−1), Φ0 (−1) = −C vˆ30 (−1).
Hence eigenvalue relation for c in terms of K and R is then Φ0 (−1) vˆ0 (−1) = 3 . Φ(−1) vˆ3 (−1) inviscid solution viscous mode Series expansion for vˆ1 and vˆ2 and hence Φ(y). To find RHS we can use the different approximations for vˆ3 (y) that were obtained. 1. Viscous boundary layer approximation for vˆ3 (y) See fig 4.8 of Drazin and Reid As R → ∞, results agree with numerics for upper branch only. No lower branch predicted. On upper branch, viscous boundary layer form is good approximation to vˆ3 at the wall. Not when the critical layer overlaps. 2. WKB approx for vˆ3 (y) See fig 4.8 of Drazin and Reid Agrees with numeric results as R → ∞, but does not give good results for moderately large R, so not good for finding Rcrit . 3. Turning point approx to vˆ3 (y) i.e. solution valid in critical layer involves Airy functions. Good agreement with numerics, somewhat surprizing. See fig 4.11 of Drazin and Reid.
48
5 5.1
Boundary layers Formulation
Non–dimensional N–S equation 1 ∂u + u · ∇u = −∇p + ∇2 u, ∂t R
∇ · u = 0.
For a 2D flow, can write velocity potential in terms of a streamfunction ψ(x, y) u = (ψy , −ψx ).
∂ ∂ ψyt + ψy − ψx ψy = −px + ∂x ∂y ∂ ∂ −ψxt − ψy ψx = −py − − ψx ∂x ∂y
1 2 ∇ ψy , R 1 2 ∇ ψx . R
(5.1.31) (5.1.32)
Eliminate p ∂ ∂ (5.1.31) − (5.1.32), ∂y ∂x ∂ ∂ ∂ 1 2 ∇ ψ + ψy − ψx ∇2 ψ = ∇4 ψ, ∂t ∂x ∂y R
(5.1.33)
(Alternative derivation is to take ∇∧ of N–S to eliminate p). Note For 2D flows ∇2 ψ = −ω, where vorticity ω = ωk. Hence we can write (5.1.33) as 1 ∂ω ∂(ψ, ω) − = ∇2 ω. ∂t ∂(x, y) R This is known as vorticity–streamfunction equation.
5.2
General theory for R 1
Outer expansion ψ = ψ0 +
1 ψ1 + O(R−2 ). R
Leading order inviscid ∂ ∂ ∂ 2 ∇ ψ0 + ψ0y − ψ0x ∇2 ψ0 = 0. ∂t ∂x ∂y Solve subject to u · n = 0 on solid surface (standard inviscid b.c), i.e. ψ =constant. In general u · s 6= 0. i.e. Flow does not satisfy the viscous b. condition that u = 0 on surface. Let s be a coordinate along the surface and u · s = U (s, t) = 49
∂ψ . ∂n
= u · s.) (Exercise: check ∂ψ ∂n Hence need a boundary layer at the wall, i.e. rescale n, n = N with 1. N is b–layer variable. Also equation is nonlinear for ψ, so also rescale ψ ˆ ψ = δ ψ. Since
∂ψ ∂ ψˆ = , ∂N δ ∂n ˆ
∂ψ → U (s), we choose δ = 1 to ensure ∂N is O(1). then as N → ∞ and n → 0, ∂ψ ∂n Next choose (as a function of R) so viscous terms enter at leading order. For simplicity consider a flat plate.
y = Y, ψ = ψˆ ∂ ∂ 1 ∂2 ∇2 = + = + O(1). ∂x2 ∂y 2 2 ∂Y 2 2
2
Subst. into (5.1.33) to give ∂ 1 ˆ ∂ 1 1 ˆ 1 ∂ ˆ ˆ ˆ ψY Y + ψY − ψx ψY Y = ψY Y Y Y + smaller terms. 2 2 ∂t ∂x ∂Y R 4 Choose R2 = 1 to balance terms. = R−1/2 , gives width of the b–layer. ψˆY Y t + ψˆY ψˆY Y x − ψˆx ψˆY Y Y = ψˆY Y Y Y . Integrate once wrt Y ψˆY t + ψˆY ψˆY x − ψˆx ψˆY Y = ψˆY Y Y + f (x, t), (check by differentiating). Let Y → ∞, need to match the outer/inviscid solution, i.e. ∂ ψˆ → U (x, t), ∂Y and so
∂ 2 ψˆ → 0 etc. ∂Y 2
Substituting ∂U ∂U +U = f (x, t). ∂t ∂x B–layer equation becomes ψˆY t + ψˆY ψˆY x − ψˆx ψˆY Y − ψˆY Y Y = Ut + U Ux . B. conds ˆ = 0) = 0 (∵ wall is streamline), ψ(Y ψˆY (Y = 0) = 0 (∵ tangential velocity = 0), ψˆY (Y → ∞) → U (x, t) (matching to outer). 50
5.3
Uniform flow past a flat plate U =1 also steady flow ψˆY ψˆY x − ψˆx ψˆY Y − ψˆY Y Y = 0.
try to separate variables, but first we need a transform for Y variable. Let ψˆ = h(X)φ(η) where η = g(x)Y, X = x. Chain rule ∂ ∂X ∂ ∂η ∂ = + , ∂x ∂x ∂X ∂x ∂η ∂ ∂ = + g0Y , ∂X ∂η ∂ g0η ∂ = + , ∂X g ∂η ∂ ∂η ∂ = , ∂y ∂Y ∂η ∂ = g . ∂η Substitute into (5.3.34). ψˆxY
! ∂ ψˆ g 0 ∂ ψˆ +η , ∂X g ∂η ∂ g0 0 = g h φ + η hφη , ∂η g 0 gh g0h 0 = g h φη + φη + η φηη . g g ∂ = g ∂η
Similarly for other terms (check these). 1 g0h g0h 1 g0h 0 0 φη h φη + φη + η φηη − φηη h φ + η φη − φηηη = 0 g g g g g (gh)0 2 h0 φ − φφηη − φηηη = 0. g2 η g B. conds φ(0) = 0, φη (0) = 0, ghφη → U as η → ∞, U φη → . gh Fn of η only Fn of X only 51
(5.3.34)
Hence choose gh = U . For flat plate – uniform flow U = 1, i.e. gh = 1. 0 Also choose hg = 1, we get ODE for φ(η)., i.e. φηηη + φφηη = 0, known as Blasius equation (see earlier handout for numerical solution). Since g = h0 and hg = 1, hh0 = 1, h2 = x + C, 2 p h = 2(x + C). For convenience, choose C = 0, g = (2x)−1/2 . Hence ψˆ = (2x)1/2 φ(η) where η = 5.3.1
Y . (2x)1/2
Notes on Blasius solution
The x component of velocity in the b–layer u = ψy = ψˆY
where η=
= ghφη , = φη (η),
Y R1/2 y = . (2x)1/2 (2x)1/2
If we return to dimensional variables ( use ∗ to denote dimensional variable). Let U ∗ be a typical velocity Let L∗ be a typical length. U* BL
Figure 34: Boundary layer flow on a flat plate U ∗ L∗ ν ∗ But there is no length scale in the problem, so L should not appear in expression for u∗ (x∗ , y ∗ ). R=
u=
u∗ , U∗
x=
x∗ , L∗
52
y=
y∗ , L∗
η=
U ∗ L∗ ν
1/2
y∗ L∗
2x∗ L∗
−1/2
=
U∗ ν
1/2
y∗ , (2x∗ )1/2
so no dependence on L∗ , as required. NB In dimensional variables ψ ∗ = (2νU ∗ x∗ )1/2 φ(η). Width of the boundary layer–δ ∗ Various measures 1. Simplest δ ∗ is point at which u is 99% of the free–stream value. u
δ∗
y*
Figure 35: Boundary layer thickness 99% measure. If φη (ηc ) = 0.99
∗
δ =
2ν U∗
1/2
x∗1/2 ηc .
Numerically ηc = 3.47,
∗
δ = 4.91
νx∗ U∗
1/2 .
2. The displacement thickness δ1∗ is defined by δ1∗
Z
0 δ0
= 0
∗
∗
∗
u∗ 1− ∗ U
dy ∗ .
Here u → U (slip velocity) as y → δ, where δ is an approximation to the boundary layer thickness. In terms of inner (or b–layer) variable
53
Matching region U*1
Outer inviscid solution boundary layer
δ
Figure 36: Matching region.
∞=R1/2 δ
Z
δ1∗
(1 − φη (η))
= 0
=
2νx∗ U∗
1/2 Z
2νx∗ U∗
1/2 dη,
∞
(1 − φη (η)) dη. 0
Integral exists because φη → 1 (in fact tends exponentially). Numerically
∞
Z
(1 − φη (η)) dη ≈ 1.2167, 0
δ1∗
≈ 1.72
νx∗ U∗
1/2 .
3. Definition of momentum thickness
θ
∗
0 δ0
u∗ u∗ = 1 − ∗ dy ∗ , ∗ U U 0 ∗ 1/2 Z ∞ 2νx φη (1 − φη (η)) dη, = U∗ 0 ∗ 1/2 νx ≈ 0.660 . U∗ Z
Each measure takes the form ∗
δ =k
νx∗ U∗
1/2 ,
where k is a dimensionless O(1) constant. Get similar results for more general boundary layers, but different k. Commonest measure to use is the displacement thickness. Note in all cases that the width of the boundary layer increases with x∗1/2 . 54
Sometimes define a Reynolds number based on boundary layer thickness ∗ 1/2 ∗ νx U U ∗δ∗ = k , Rδ ∗ = ∗ ν U ν ∗ ∗ 1/2 U x = k , ν 1/2
= kRx∗ . Rδ∗ is Reynolds number based on b–layer thickness. Rx∗ is Reynolds number based on b–layer thickness downstream.
5.4
Flow past flat plate with arbitrary slip velocity U (x)
e.g. Wind tunnel with curved walls.
Figure 37: Wind tunnel with curved walls. Same PDE in b–layer as before ψˆY → U (x) as Y → ∞. In general, can’t find transform η = g(x)Y so solution is seperable, but with suitable choice of g(x) can simplify the equation, ψˆY ψˆY x − ψˆx ψˆY Y − ψˆY Y Y = U Ux . Let η = g(X)Y, X = x, ψˆ = h(X)φ(X, η). Chain rule ηg 0 ˆ ψˆx = ψˆX + ψη , g ηg 0 = h φ + hφX + hφη , g extra term cf U = 1 case 0
55
Same analysis as before, but 2 extra terms φηηη +
h0 (gh)0 h U UX φφηη − 2 φ2η + (φX φηη − φηX φη ) = − 3 . g g g g h
Boundary condition U (X) . gh
φη → Choose g and h to simplify equation,
h0 = 1, g
gh = U so g = h0 h0 h = U, Z X h2 U dx, = 2 c
convenient to let c = 0,
so Z h = 2
1/2
X
,
U dx
0
Z g = U 2
−1/2
X
U dx
.
0
Hence, φηηη + φφηη + β 1 − φ2η + γ (φX φηη − φηX φη ) = 0, where RX 2 0 U dx dU (gh)0 UX β = = 2 = , g2 g U2 dX RX 2 0 U dx h U γ = = 2 = . g g U2 NB Often transform the streamwise coordinate, X, Z X ξ= U (x) dx, 0
to put the equation in standard form φηηη + φφηη + β 1 − φ2η + 2ξ (φX φηη − φηX φη ) = 0, where β is a function of ξ. β is known as the pressure gradient parameter (or Hartree parameter).
56
NB Pressure gradient is dp = −U UX , dX (from non–dimensional Navier–Stokes or Bernoulli). Assuming U > 0. dp < 0 known as favourable pressure gradient – flow is accelerated downstream. If UX > 0, dX p β>0 High
Low
x
Figure 38: Figure of pressure gradient. If UX < 0, adverse pressure gradient. 5.4.1
Falkner–Skan solutions
If β =constant, then solution with
∂φ ∂X
= 0 exists.
φηηη + φφηη + β 1 − φ2η = 0. True if U (X) = AX n (check) β=
2n . n+1
Aside Do any other U (X) give β =constant? Need R U dx UX = β0 UZ2 β0 U 2 U dx = , UX U 2 UXX U = β0 2U − , UX2 U 2 UXX 1 = 2− U, UX2 β0 U UXX = kUX2 . 57
Nonlinear 2nd order ODE. Try U = eµX , µ2 e2µX = kµ2 e2µX , i.e. solution for k = 1. Any more? 5.4.2
Numerical results for Falkner–Skan
Method 1. φ(0) = φη (0) = 0 choose φηη (0) = a. 2. Solve ODE by integration from η = 0. If φη (ηmax ) 6= 1 ηmax 1, then adjust a. Summary of Results Solutions do exist with Umax > 1, i.e. Disregard on physical grounds. φη 1
η
Figure 39: Falkner–Skan profile which is neglected. 1. If n > 0 =⇒ β > 0, Unique solution for φη (η). 2. −0.0904 < n < 0 =⇒ −0.1978 < β < 0, Two solutions exist, one with φηη (0) > 0, and one with φηη (0) < 0, corresponds with flow reversal. 3. β < −0.1978 No solutions exist. i.e. solutions exist for favourable pressure gradients and small adverse pressure gradients. See MATLAB handout A.4.1. 58
η
η
Flow reversal
φ
ηη
(0)>0
φη
φ
ηη
(0)<0
φη
Figure 40: Falkner–Skan profiles for −0.1978 < β < 0. 5.4.3
Summary of Falkner–Skan profiles
1. Adverse pressure gradients (β < 0) thicken boundary layer. 2. φηη (0) < 0 =⇒ flow reversal at wall (rarely occurs –see later). Such flows formally exist for small adverse pressure gradients. 3. As β decreases from positive values down to −0.1978 uy (0) ∝ φηη (0) decreases to zero. 4. Boundary layer separation
Figure 41: Boundary layer separation. At simplest level, flow reversal signals b–layer separation.
Figure 42: Boundary layer separation. 5. In general, solve φ(X, η) as a PDE marching downstream from X = 0. 59
If φηη (η = 0) approaches zero, we expect separation to occur. Example Rankine body – Source in a uniform flow ξ
Figure 43: Rankine body. Can find slip velocity on surface U (ξ)
β
ξ β min
Figure 44: Rankine body pressure gradient. Results suggest if βmin > −0.1978 b–layer stays attached. βmin can be slightly less than −0.1978 and stay attached (due to presence of φξ terms) but larger -ve values =⇒ separation.
5.5
Stability of boundary layers
In §2 and 3 we considered stability of steady flows of the form u = U (y)i, i.e. strictly parallel flows. For viscous flows we showed that only strictly parallel flows were U (y) = Ay 2 + By + C. Two flows in particular are commonly studied.
60
1. Poiseuille Flow U (y) = 1 − y 2 . Neutral stability curve,
K Upper branch UNSTABLE lower branch R
R crit
Figure 45: Poiseuille flow neutral stability curve. 2. Couette flow U (y) = y
− 1 ≤ y ≤ 1.
Corresponds to flow between 2 plates moving in opposite directions. Couette flow is stable, Now consider boundary layer on a flat plate – zero pressure gradient (i.e Blasius boundary layer). Using same notation as before ψ = R−1/2 ψˆ
y = R−1/2 Y, Y ψˆ = (2x)1/2 φ(η) η = . (2x)1/2 Let u = U i + V j, ∂ψ ∂ ψˆ = = φη (η), ∂y ∂Y ˆ ∂ψ Y −1/2 ∂ ψ −1/2 −1/2 1/2 = − = −R = −R (2x) φ − (2x) φη . ∂x ∂x (2x)3/2
U = V
Note that V = O(R−1/2 ). i.e. For R 1 =⇒ V U , i.e. the flow is almost parallel. Now consider a point dimensional distance x∗ from leading edge of plate. Recall there is no natural length scale in problem for non–dimensionalisation of Navier–Stokes. Choose L∗ = x∗ so U ∗ L∗ U ∗ x∗ R= = . ν ν 61
Hence at fixed point
x∗ x = ∗ = 1, L
and
Y η=√ . 2 Then repeat stability analysis for strictly parallel flows, ignoring the j component of velocity. (can be built into theory as a correction – non–parallel correction). Basic structure is the same, with viscous wall layer and critical layer inside the steady boundary layer. Obtain similar results U = φη (η)
K Upper branch UNSTABLE lower branch R
R crit
Figure 46: Blasius boundary layer, neutral curve. Main difference in analysis is that instead of having 4 b. conditions at the 2 walls as in Poiseuille flow, here extra conditions are that vˆ decays exponentially for large y. Results See fig 4.22 Drazin and Reid R1 = Reynolds number based on displacement thickness. R1 =
U ∗δ∗ 1/2 = 1.72Rx∗ . ν
where the notation is α1 ∝ K. Fig 4.22 in neutral stability curve for Blasius flow. Fig 4.24 – illustrates effect of pressure gradient (i.e. β) on critical Reynolds number. Adverse pressure gradients reduce critical R. i.e. Adverse pressure gradients =⇒ b–layer is more unstable.
5.6
Comparison of experiment and theory for neutral curve for flat plate
In developing theory for stability we considered disturbance of the form ei(kx−ωt) , 62
(ignoring z–dependency, i.e. 2D theory). Take k real and found a dispersion relation, formally D(k, ω) = 0, k real, solved D(k, ω) = 0 to give ω = ωr + iωi , ei(kx−ωr t) eωi t . So if disturbance is known at t = 0, we can find whether disturbance grows or decays with time. This is known as TEMPORAL STABILITY. In experiments it is hard to set up conditions so disturbance is known everywhere at t = 0. Easier to force a disturbance at a particular point. e.g. Oscillating strip on a plate, or a vibrating ribbon in the flow above the plate.
x=0
x=x
0
Figure 47: Oscillating strip or vibrating ribbon. In this case, we know disturbance at x = x0 for all time. e.g. oscillating at frequency ω, say u(x0 , t) = Aeiωt . Hence, suggests we take ω real and solve the dispersion relation D(k, ω) = 0 to give a complex wavenumber k. ei(kx−ωt) = ei(kr x−ωt) e−ki x . ki > 0, disturbance decays downstream
STABLE
x x0 Figure 48: Spatially stable. ki < 0, disturbance grows downstream This is known as SPATIAL STABILITY.
63
UNSTABLE
x x0
Figure 49: Spatially unstable. Theory Neutral curve for spatial stability is obtained in similar method as for temporal. ω
unstable
x
Figure 50: Neutral curve. Experiments Schubauer and Skramstad (1947) Wind tunnel test in general. Test section Summary of experiment 1. Incoming flow as uniform as possible (uses meshes etc to cut turbulence and lining on walls of tunnel to cut noise/vibration). 2. Vibrating ribbon – known frequency (ω) and amplitude. 3. Hot wire probes accurately measure flow speed close to plate. Results See figs 4.29 and 4.30 (Drazin and Reid)
64
meshes to reduce turbulence
Test section
Figure 51: Wind tunnel experiment.
hot wire probes vibrating ribbon
Figure 52: Test section. Parameters F =
ω Rδ∗
U ∗δ∗ , = ν 1/2 = 1.72Rx∗ ,
R
δ∗
where R
x∗
U ∗ x∗ = . ν
Parallel theory vs Non–parallel theory
δ ∗ = x∗1/2 . Parallel theory – Orr–Sommerfeld equation. non–parallel theory – extra terms in Orr–Sommerfeld.
65
δ∗ parallel theory not so good (significant growth in b−l thickness)
x* flow almost parallel expect parallel theory to be good
Figure 53: Non–parallel b-layer.
66
6
Laminar/Turbulent transition on a flat plate
6.1
Introduction
Flat plate in uniform flow
Laminar flow
Turbulent flow Transition point
Figure 54: Flat plate in uniform flow. Detection of transition point 1. Direct measurement of flow speed (using hot wire probe). 2. Visualisation – e.g. introduce smoke traces. 3. Surface measurements e.g. pressure/shear. From experiments position of transition is very sensitive to experimental setup. 1. Smoothness of plate. 2. Turbulent level of incoming flow. 3. Noise/acoustic disturbances in tunnel.
6.2
Stability theory
Consider U and ν fixed, then we can draw the neutral stability curve in dimensional x∗ , ω ∗ plane (∗ denotes dimensional quantities). (see spatial stability theory). ω∗ ω∗m
stable
ω∗0
unstable
x* x*1 0
x*
2
x*
Figure 55: Neutral stability curve for spatial stability. 67
Consider a vibrating ribbon frequency ω0∗ positioned at x∗0 . For x∗0 < x∗ < x∗1 disturbance decays exponentially. For x∗1 < x∗ < x∗2 disturbance grows exponentially. For x∗2 < x∗ < ∞ disturbance decays exponentially. Amplitude
initial amplitude
x*0
x*
x*
x*2
1
Figure 56: Schematic of disturbance amplitude. But in deriving Orr–Sommerfeld equation, we assumed the unsteady disturbance was small compared to steady flow. If amplitude of the unsteady disturbance grows large, linear equation is no longer valid and we need a nonlinear equation. With linear equation, solution ∝ eiωt continue to be ∝ eiωt . With nonlinear equation, these terms give rise to e2iωt , e3iωt etc (i.e. higher harmonics) and it is this that eventually leads to turbulence in the b–layer. Hence transition point loosely corresponds to onset of nonlinear behaviour. Expect transition point to lie downstream of neutral stability point, assuming that initial forcing is small. Amplitude
nonlinearity at this amp
A B C x*A x*
0
x*
1
x*
B
x*2
x*
Figure 57: Amplitude schematic. x∗A =transition point for disturbance A. Large initial forcing (by vibrating ribbon etc) =⇒ transition point closer to neutral stability point. Notes 1. Wavelength of disturbance is short at lower branch – important later. 2. Unstable waves in b–layers are called Tollmien–Schlichting waves (T–S waves). 68
3. In practical cases, we don’t have a forcing in the b–layer, but instead b–layer is forced by unsteady disturbances in the free–stream (e.g. noise). So to find transition point, need to determine how T–S waves are generated by the interaction of free–stream and the b–layer. known as boundary layer receptivity.
6.3
Receptivity theory
Consider
∗
u∗∞ = U ∗ (1 + δe−iω t )i δ 1, as velocity far from flat plate.
Figure 58: Velocity far from plate. Effectively, acoustic wave approaching the plate. How does the unsteady disturbance generate T–S waves. ∗ Length scale from disturbance l∗ = Uω∗ . Reynolds number based on this length Re =
U ∗ l∗ U ∗2 = . ν νω ∗
Assume Re 1. From earlier, the wavelength of T–S waves is short, so generation of T–S waves occurs when steady flow solution varies on short streamwise scale. 1. e.g. leading edge of plate – rapid b–layer growth
Figure 59: Leading edge receptivity. 2. Roughness on plate
69
Figure 60: Roughness receptivity.
6.4
Leading edge receptivity
See Goldstein - Journal of Fluid Mechanics (1983) 127, and my notes on this research paper.
6.5
Receptivity due to roughness element
Roughness centred at x∗ = x∗0 , y ∗ = H ∗ f (x∗ − x∗0 )
f (x) → 0 as |x| → ∞.
H*
x*
0
Figure 61: Roughness element. We know the wavelength of a disturbance near to the lower branch, x∗n (ω ∗ ), is given by λ∗ = O(3 x∗n ). This ‘suggests’ that unstable disturbances will be generated by interaction of unsteadiness in the free–stream (ω ∗ ) with the roughness on the surface of the plate, if the roughness has a similar length scale to wavelength of the unsteady disturbance (i.e. width O(3 x∗n )). Consider roughness close to neutral point (or lower branch) i.e. x∗0 = O(x∗n ). ∗ With this length scale we obtain the same triple deck structure as before. Here y = xy ∗ . 0 If we take H ∗ = 5 hx∗0 , then roughness has same scale as the inner or wall layer. If h 1, then roughness is small compared to the wall layer – means we can solve linearised equations in wall layer. Boundary conditions on the plate u(y ∗ = H ∗ f (X)) = 0, v(y ∗ = H ∗ f (X)) = 0. 70
u*=U*(1+ δ e
−iω∗ t*
)
x*
0
3
ε x* 0
Figure 62: Roughness element. ~
y= ε3y
upper or outer
^
y= ε4y
main
5
y= ε Y
wall layer
x*0 X=O(1)
X=x*−x*0 ε 3 x*0
Figure 63: Roughness element triple deck. In terms of non–dimensional parameters in lower deck u(Y = hf (X)) = v(Y = hf (X)) = 0, where h 1 and f (X) = O(1). We can expand u and v as a Taylor series u(Y = hf ) = u(0) + hf u0 (0) + O(h2 ). Boundary condition becomes u(0) = −hf (X), u0 (0) and similarly for v(Y ). Find that, if: 1. Roughness is upstream of neutral point then wave decays with distance downstream but then grows when lower branch point is reached. 2. Roughness is downstream of neutral point, then disturbance grows. Roughness elements close to lower branch point create unstable disturbances with highest magnitude =⇒ earliest transition.
71
2) 1)
x*
x*n
0
x*0
Figure 64: Schematic of disturbance amplitude.
Ω2 Ω1 Fluid
Figure 65: Flow between two cylinders.
6.6
Importance of stability properties
1. Mixing of fluids u(r) is velocity in θ direction say. Stability depends on Ω1 , Ω2 , r1 , r2 . No mixing of fluid if stable. 2. Drag on a body Drag force: Inviscid (known as form drag) F = b–layer is turbulent.
R
ρ ds. Viscous contribution from b–layer, larger if
Form drag is related to size of wake behind body (large wake =⇒ large drag). (a) Streamline body 72
u Ω 2 r2 Ω 1 r1
r2
r1
r
Figure 66: Velocity of fluid.
Figure 67: Streamline body. (b) Blunt body
Figure 68: Blunt body. In (b) if we want to minimise drag, we want boundary layer to be turbulent to delay separation and makes wake narrower. Decrease in form drag outweighs small increase in boundary layer drag. i.e. golf ball rough to make boundary layer turbulent and reduce total drag. For streamline bodies (i.e. (a)) the viscous boundary layer drag outweighs form drag, so try to reduce it by making boundary layer laminar for as long as possible. i.e. Delay transition point. Typical techniques – smooth surfaces. also suction through the surface stabilises the flow.
73
A A.1 A.1.1
MATLAB Worksheets Kelvin–Helmholtz instability Kelvin–Helmholtz general stability
Plotting the contours of constant growth rate for the Kelvin–Helmholtz instability. clear gamma=0.5; [k,l]=meshgrid(0:0.01:2.0,-2.0:0.01:2.0); f=k.^2-gamma*(k.^2.0+l.^2.0).^(0.5); contour(k,l,f,[0,0.5,1.0]) 2
1
l 0 −1
−2 0
0.5
1 k
1.5
2
Figure 69: Kelvin–Helmholtz growth rates γ = 0.5
gamma=1.0; For fixed k > γ, largest growth rate is for l = 0, i.e. 2D disturbances are most unstable.
74
2
1
l 0 −1
−2 0
0.5
1 k
1.5
2
Figure 70: Kelvin–Helmholtz growth rates γ = 1.0
75
A.2 A.2.1
Inviscid Stability Inviscid stability of mixing layer
Plotting the value of c2 as a function of X = Kb. clear x=0:0.01:5; zero=0*x; f=(1-0.5./x).^2-0.25*exp(-4*x)./x.^2; plot(x,f,x,zero,’k’) 1
0.5
f(X) 0 −0.5
−1 0
1
2
3
4
X
Figure 71: Mixing layer stability.
76
5
A.2.2
Unbounded piecewise linear mixing layer
A MATLAB code to solve the piecewise linear mixing layer problem. clear clc syms y b P Q R S k c bb=1.0; v1=P*exp(k*y); dv1=diff(v1,y); y1=-b; u1=-1; du1=diff(u1,y); v2=Q*exp(k*y)+R*exp(-k*y); dv2=diff(v2,y); y2=b; u2=y/b; du2=diff(u2,y); v3=S*exp(-k*y); dv3=diff(v3,y); u3=1; du3=diff(u3,y); eqn1=subs(v1,y,y1)-subs(v2,y,y1); eqn2=subs(v2,y,y2)-subs(v3,y,y2); eqn3=subs(dv2,y,y1)-subs(dv1,y,y1)-subs(v1,y,y1)*(du2-du1)/(subs(u1,y,y1)-c); eqn4=subs(dv3,y,y2)-subs(dv2,y,y2)-subs(v2,y,y2)*(du3-du2)/(subs(u2,y,y2)-c); [soln_P,soln_Q,soln_R]=solve(eqn1,eqn2,eqn3,P,Q,R); disp=expand(subs(eqn4,{P,Q,R},{soln_P,soln_Q,soln_R})/S); disp2=disp(1); c1=solve(disp2,c); for j=2:1001 k1(j)=(j-1.)/1000.; c11r(j)=double(subs(real(c1(1)),{b,k},{bb,k1(j)})); c11i(j)=double(subs(imag(c1(1)),{b,k},{bb,k1(j)})); c12r(j)=double(subs(real(c1(2)),{b,k},{bb,k1(j)})); c12i(j)=double(subs(imag(c1(2)),{b,k},{bb,k1(j)})); end figure(1) plot(k1,c11r,’k’,k1,c12r,’k’) 77
figure(2) plot(k1,c11i,’k’,k1,c12i,’k’) 0.5
1
0.5
cr
ci 0
0
−0.5
−0.5 0
0.2
0.4
0.6
0.8
1
K
−1 0
0.2
0.4
0.6
K
Figure 72: Unbounded mixing layer stability.
78
0.8
1
A.2.3
Two part piecewise linear boundary layer flow
A MATLAB code to solve the two part piecewise linear boundary layer problem. clear clc syms y b P Q R S k c bb=-0.2; v1=P*(exp(k*y)-exp(-k*y)); dv1=diff(v1,y); y1=0.5; u1=2.0*b*y; du1=diff(u1,y); v2=Q*exp(k*y)+R*exp(-k*y); dv2=diff(v2,y); y2=1.0; u2=2*(1-b)*y+(2*b-1); du2=diff(u2,y); v3=S*exp(-k*y); dv3=diff(v3,y); u3=1.0; du3=diff(u3,y); eqn1=subs(v1,y,y1)-subs(v2,y,y1); eqn2=subs(v2,y,y2)-subs(v3,y,y2); eqn3=subs(dv2,y,y1)-subs(dv1,y,y1)-subs(v1,y,y1)*(du2-du1)/(subs(u1,y,y1)-c); eqn4=subs(dv3,y,y2)-subs(dv2,y,y2)-subs(v2,y,y2)*(du3-du2)/(subs(u2,y,y2)-c); [soln_P,soln_Q,soln_R]=solve(eqn1,eqn2,eqn3,P,Q,R); disp=expand(subs(eqn4,{P,Q,R},{soln_P,soln_Q,soln_R})/S); disp2=disp(1); c1=solve(disp2,c); for j=2:3001 k1(j)=(j-1.)/1000.; c11r(j)=double(subs(real(c1(1)),{b,k},{bb,k1(j)})); c11i(j)=double(subs(imag(c1(1)),{b,k},{bb,k1(j)})); c12r(j)=double(subs(real(c1(2)),{b,k},{bb,k1(j)})); c12i(j)=double(subs(imag(c1(2)),{b,k},{bb,k1(j)})); end for j=1:2000 yy(j)=(j-1.0)/1000.0; 79
if (yy(j)
2
1
1
0.5
0
0
0.5
0 0 2
1
0.5 0 −0.5 0
2
4
1
2
1 2
0 0
4
Figure 73: b = 0.2.
80
2
1
1
0.5
0
0
0.5
0 0 2
1
0.05 0
2
4
1
2
1
−0.05 0
2
0 0
4
Figure 74: b = 0.49.
2
1
1
0.5
0 1
0 0 2
0
1
−1 0
0
0.5
2
1
0 0
4
Figure 75: b = 0.8.
81
2
4
1
2
2
1
1
0
0
0
0.5
−1 0 2
1
0.5 0 −0.5 0
2
4
0
5
1 2
0 −5
4
Figure 76: b = −0.2.
82
A.2.4
Blasius boundary layer solution
Define a function file ‘blasius.m’ containing the ODE information function ysol=blasius(x,y) ysol=[y(2); y(3); -y(1)*y(3)]; Then solve ODE using clear xspan=[0 6]; y0=[0; 0; 0.4696]; [x,y]=ode45(@blasius,xspan,y0); zero=0*y(:,1); plot(x,y(:,2),x,zero,’k’) 1.5
1
0.5
0 0
2
4
Figure 77: Plot of U (y)
plot(x,y(:,3),x,zero,’k’)
83
6
0.5 0.4 0.3 0.2 0.1 0 0
2
4
Figure 78: Plot of U 0 (y).
84
6
A.2.5
Stability analysis for U = tanh y - numerical results
For most U (y), Rayleigh’s equation (U − c)(D2 − K 2 )v − U 00 v = 0, can not be solved exactly. Here we consider U = tanh y in the unbounded region −∞ < y < ∞. The boundary conditions on v(y) are v → 0 as |y| → ∞. Since U 00 /(U − c) → 0 as |y| → ∞ the conditions of decay on v(y) become v → C− eKy v → C+ e−Ky which can be written as
v 0 (y) → v(y)
as
y → −∞, as y → ∞,
K y → −∞ . −K y → ∞
Thus for given K we must find c so that these boundary conditions are satisfied. This turns out to be a difficult numeric problem since c is typically complex. Here I summarise the method used and then present some results. General numerical approach
1. We solve the ODE in the interval y ∈ [−ym , ym ] starting from y = −ym where ym is large. I used ym = 5. The initial conditions then become v(−ym ) = ,
v 0 (−ym ) = K,
where the value of does not matter since the equation is linear. 2. We then choose a value of c, integrate (i.e. solve) the equation numerically and then find the value of v 0 (ym )/v(ym ). If the correct value of c has been chosen, this value should be close to −K, otherwise a new value of c is chosen. 3. The subtle part of this calculation is in choosing a new value of c using the previous results. For this particular problem there is some simplification possible. v(y) = sechy is a solution of Rayleigh’s equation if c = 0 and K = 1. For 0 < K < 1 it turns out that Re(c) = 0 with Im(c) a decreasing function of K. (NB It is possible to prove that if the unstable eigenvalue is unique, then the real part of c is independent of K). Thus to find c(K) we let c = λi and reduce λ from 1 until v 0 (ym )/v(ym ) ≈ −K. The MATLAB scritps are function vsol=rayleigh(y,v,k,c) vsol=[v(2); -2.0*tanh(y)*sech(y)^2.0*v(1)/(tanh(y)-c)+k^2.0*v(1)]; and
85
clear epsilon=0.1; lambda=1.0; options = odeset(’AbsTol’,1e-8,’RelTol’,1e-8); j=1; for k=0.05:0.05:0.95 a=1000.0; b=1.0; min=10000000.0; while (abs(a/b+k)
figure(1) plot(kk,lambdalambda) figure(2) plot(kk,growth) The results for c(K) = λi 0 < K < 1 are given below: 86
K 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00
λ(K) Kλ(K) 1.0000 0.0000 0.9157 0.0458 0.8367 0.0837 0.7644 0.1147 0.6975 0.1395 0.6352 0.1588 0.5769 0.1731 0.5221 0.1827 0.4705 0.1882 0.4215 0.1897 0.3750 0.1875 0.3307 0.1819 0.2883 0.1730 0.2477 0.1610 0.2086 0.1460 0.1710 0.1283 0.1347 0.1078 0.0995 0.0846 0.0654 0.0589 0.0323 0.0307 0.0000 0.0000
Table 1: Stability results for U = tanh y.
87
1
0.2
0.8
0.15
0.6 0.1 0.4 0.05
0.2
(a)
0 0
0.2
0.4
0.6
0.8
1
(b)
0 0
0.2
0.4
0.6
0.8
1
Figure 79: (a) Im(c)(K) and (b) growth rate as function of K. Thus we see that the system is unstable for 0 < K < 1. In fact the system is stable (Im(c) < 0) for K > 1 but finding the eigenvalues for K > 1 turns out to be quite complicated. Recall that K = (k 2 + l2 )1/2 and that the growth rate is given by ω(k, l) = kIm(c) = kλ(K). The numerical results show that λ is monotonically decreasing function of K and so for a given k, the maximum growth rate occurs when l = 0. Using the commands in the MATLAB worksheet below, you can verify that these eigenvalues are correct. Choose a value of K and compare the eigensolution v(y) when c = c(K) with v(y) when c take other values. We can plot the eigensolution using clear epsilon=0.1; k=0.6; c=0.2883*i options = odeset(’AbsTol’,1e-8,’RelTol’,1e-8); yspan=[-5 5]; v0=[epsilon; epsilon*k]; [y,v]=ode45(@rayleigh,yspan,v0,options,k,c); last=length(y); plot(y,v(:,1)) For the case K = 0.6 we have found the approximate value of c = 0.2883i we find an eigensolution which decays exponentially at large y, but if we choose an incorrect value of c, 0.4i say, then the eigensolution consists of both eKy and e−Ky parts of the solution, and does not tend to zero for large y.
88
4 1.2 1
3
0.8 2
0.6 0.4
1
0.2
(a)
0 −5
0
5
0
(b)−5
0
Figure 80: (a) c = 0.2883i and (b) c = 0.4i for K = 0.6.
89
5
A.3 A.3.1
Viscous Stability Eigenvalues for Poiseuille flow.
Finding the eigenvalues for Poiseuille flow U (y) = 1 − y 2 for y ∈ [−1, 1] with K = 1.0 and ˆ = 10000. R % Matlab script to find global eigenvalues for Poiseuille flow clear N=121; R=10000.0; alpha=1.0; % The collocation points are at ybar=cos(m*pi/(N-1)) for m=0..N-1 % Setup matrix t of Chebyshev polynomials and their derivatives at each collocation % point for m=N-2:-1:2; ybar=cos(m*pi/(N-1)); t=0.0; t(1,1)=1.0; t(1,2)=ybar; for ii=2:N-1; t(1,ii+1)=2.0*ybar*t(1,ii)-t(1,ii-1); end for j=2:5; t(j,1)=0.0; t(j,2)=t(j-1,1); t(j,3)=4.0*t(j-1,2); for k=4:N; t(j,k)=2.0*(k-1.0)*t(j-1,k-1)+(k-1.0)/(k-3.0)*t(j,k-2); end end % Evaluate the base flow at value of ybar U=1.0-ybar^2.0; dU=-2.0*ybar; d2U=-2.0; % Setting up matrices of Orr-Sommerfeld equation
for j=1:N; a(N-m,j)=U*(t(3,j)-alpha^2*t(1,j))-d2U*t(1,j)-1.0/(i*alpha*R)*(t(5,j)-2.0*alpha^2*t b(N-m,j)=t(3,j)-alpha^2*t(1,j); end end 90
% Boundary conditions for j=1:N; a(1,j)=1.0; a(2,j)=(j-1.0)^2.0; a(N-1,j)=(-1.0)^(j-2.0)*(j-1.0)^2.0; a(N,j)=(-1.0)^(j-1.0); b(1,j)=0.0; b(2,j)=0.0; b(N-1,j)=0.0; b(N,j)=0.0; end % find eigenvalues c [V,D]=eig(a,b); % put eigenvalues into plottable form for j=1:N; realc(j)=real(D(j,j)); imagc(j)=imag(D(j,j)); end plot(realc,imagc,’+’,[0 1],[0 0]) axis([0 1 -1 0.1])
0 −0.2 −0.4 −0.6 −0.8 −1 0
0.2
0.4
0.6
0.8
1
ˆ = 10000.. Figure 81: Temporal eigenvalues for plane Poiseuille flow with K = 1 and R
91
A.3.2
Asymptotic solution to ODE
The exact solution to
dy + y = x, dx with x ≥ 0 and 0 < 1 with one b. cond y(0) = 1, is
x
yexact = x − + (1 + )e− . The 1st and 2nd approximations to the outer solution are (1)
youter = x, (2)
youter = x − , while the 1st and 2nd approximations to the inner solution are (1)
yinner = e−x/ , (2) yinner
−x/
= e
+ e
−x/
x + −1 .
The 5 solutions can be plotted using the code below clear clc eps=0.1; x=0:0.025:1; y=x-eps+(1+eps)*exp(-x/eps); yo1=x; yo2=x-eps; yi1=exp(-x/eps); yi2=exp(-x/eps)+eps*(exp(-x/eps)+x/eps-1); plot(x,y,’k+’,x,yo1,’r’,x,yo2,’b’,x,yi1,’g’,x,yi2,’c’)
92
1 0.8 0.6 0.4 0.2 0 −0.2 0
0.2
0.4
0.6
0.8
1
Figure 82: Asymptotic approximations to the exact solution.
93
A.4 A.4.1
Boundary Layers Falkner–Skan velocity profiles
To plot Falkner–Skan profiles for pressure gradient β = 2n/(n + 1), we solve f 000 + f f 00 + β(1 − f 02 ) = 0. To do this, create the file f_skan.m , function vsol=f_skan(x,f,n) vsol=[f(2); f(3); -f(1)*f(3)-2.0*n/(n+1)*(1.0-f(2)^2.0)]; and use along with clear n=1.0; options = odeset(’AbsTol’,1e-8,’RelTol’,1e-8); fspan=[0 10]; f0=[0.0; 0.0; 1.23259]; [y1,f1]=ode45(@f_skan,fspan,f0,options,n); last=length(y1); %f1(last,2) n=0.3; fspan=[0 10]; f0=[0.0; 0.0; 0.9002]; [y2,f2]=ode45(@f_skan,fspan,f0,options,n); n=0.0; fspan=[0 10]; f0=[0.0; 0.0; 0.4696]; [y3,f3]=ode45(@f_skan,fspan,f0,options,n); n=-0.05; fspan=[0 10]; f0=[0.0; 0.0; 0.3098]; [y4,f4]=ode45(@f_skan,fspan,f0,options,n); 94
n=-0.05; fspan=[0 10]; f0=[0.0; 0.0; -0.1418]; [y5,f5]=ode45(@f_skan,fspan,f0,options,n); figure(1) plot(y1,f1(:,2),’k’,y2,f2(:,2),’k’,y3,f3(:,2),’k’,y4,f4(:,2),’k’) figure(2) zero=f5*0.0; plot(y4,f4(:,2),’k’,y5,f5(:,2),’k’,y5,zero,’k’) The value of f 00 (0) is found by trial and error by searching for when f 0 (10) ≈ 1.
1.2
1
1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0 0
2
4
6
−0.2 0
2
4
6
8
10
Figure 83: Falkner–Skan profiles for m = −0.05, 0, 0.3 and 1, and the two profiles for m = −0.05.
95
A.5
Transition
A.5.1
Summary of Goldstein (1983) JFM U ∗2 U∗ , Based on lengthscale , νω ∗ ω∗ = Re−1/6 ,
Re =
Note: Here we use ∗ to denote dimensional quantities. 2 regions – see Figure 1. 1.
∗
x = O(1) i.e x = O
U∗ ω∗
RECEPTIVITY REGION
Boundary layer is growing so near–parallel flow assumption cannot be made. Streamfunction for unsteady disturbance satisfies ∂ ∂ ∂ 00 0 3/2 (F ψ0η ) + 2x (F ψ0 − F ψ0η ) + 2ixψ0η = (2x) i− u∞ (x), ψ0ηηη + ∂η ∂x ∂x with ψ0 = ψ0η = 0 on η = 0 and ψ0η → (2x)1/2 u∞ (x) as η → ∞, where here we take u∞ = 1. These are eqns. (3.2)–(3.5). The solution takes the form ψ0 = ψST +
∞ X
(i)
Ci ψAP ,
i=1 (i)
where ψST is determined by the local flow (Stokes’s solution) and ψAP are eigensolutions which satisfy ψ0η → 0 as η → ∞. Ci are numerical coefficients determined by forcing near to the leading edge. Key thing about the eigensolutions is ψAP
λ(2x)3/2 = Cg0 (x, η) exp − 3U00
as x → ∞,
which is (3.6), where λ = =
e−7πi/4 3/2 ζn −πi/4
e
3/2
ζn = ρn eiπ ,
ρn > 0,
.
ρn
Hence Re(λ) > 0 and ψAP decreases in amplitude with distance downstream.
wavelength of disturbance
96
ψ , ψx ∼ x−1/2 ,
∼
(A.5.35)
(e.g. ψ = eikx , λ ∼ 1/k). Hence wavelength shortens with distance downstream. So terms in (2.19) ignored to obtain (3.2) e.g. 6 x∂ 2 /∂x2 , eventually become same size as terms retained. So (3.2) no longer valid when x = O(−2 ). 2. −2
∗
x = O( ) x = O
−2 U ∗ ω∗
ORR − SOMMERFELD or TRIPLE − DECK REGION. x1 = 2 x
(see 4.3 with r = 2), Z x i s κ(x1 , ) dx , ψ = G exp 0 wavelength (non − dim) λ ∼ , κ −1/2
(cf λ ∼ x−1/2 = x1
from (A.5.35)).
Variation of κ(x1 ) given in figure 3, from numerical solution of Orr–Sommerfeld. As → 0, stable for x˜1 < 3.03, and unstable for x˜1 > 3.03 (from figure 3(b)). Note: Neutral point x∗N
= 0.344
U∗ ω∗
U∗ νω ∗
1/3 (A.5.36)
From figure 3(a) see wavelength shortening up to x˜1 ≈ 0.9 and matching to first eigensolution as x˜1 → 0 (overlap region between receptivity and Orr–Sommerfeld). Hence first eigensolution matches onto T–S wave which starts to grow at neutral point (lower branch of neutral curve). So large C1 =⇒ Unstable wave larger at neutral point, =⇒ Nonlinear effects important soon after lower branch neutral point, =⇒ Earlier transition. C1 depends on nature of free–stream disturbance–determined numerically, known as Receptivity coefficient. Equation (A.5.36) gives the equation of lower branch of neutral curve in dimensional (ω ∗ , x∗ ) − −plane (valid for large x∗ ). Lower branch in α, R plane (i.e. Notation of Orr–Sommerfeld equation (4.10)) given by (4.47)– large R–limit. Triple–Deck Asymptotic solution of OSE as R → ∞ can be analysed as 3 different asymptotic forms in 3 layers:
97
Main: η =
y−3 = O(1), (2x)1/2 1/2
=⇒ y = O(3 x1/2 ) = O(2 x1 ), so y = O(2 ) at lower branch, y∗ = O(4 ). or x∗ Outer: η = O(−1 )
y = O()
y∗ = O(3 ), x∗
y = O(3 )
y∗ = O(5 ), ∗ x
Wall layer η = O() Wavelength
λ∗ =O = O(3 ). λ = O() ∗ x x These are conventional triple deck scales – non–dimensionalised by downstream position. (Varies from Goldstein only at higher orders).
98
