A tangential shear stress σθ*r* acts on a surface element with a radial normal, which is expressed by (A.5)

$$
\sigma\_{\theta r} = \mu \left( \frac{\partial v\_{\theta}}{\partial r} - \frac{v\_{\theta}}{r} + \frac{1}{r} \frac{\partial v\_{r}}{\partial \theta} \right) = \mu \left( \frac{\partial v\_{\theta}}{\partial r} - \frac{v\_{\theta}}{r} \right) = \mu \left. r \frac{\partial}{\partial r} \left( \frac{v\_{\theta}}{r} \right) \right|. \tag{3.43}
$$

Combining (3.41) and (3.43), we obtain

$$
\sigma\_{\theta r} = -\frac{2B\mu}{r^2} \,\, . \tag{3.44}
$$

Next we calculate the viscous moment, *M*, that acts on the interior cylinder per unit axial length. This moment is equal to the component σθ*r* evaluated at*r* = *R*1 and the area, 2π *R*1, on which this stress acts, multiplied by the lever arm, *R*1, the distance between the axis and the point where the force acts. We have

$$M = -2\pi R\_1^2 \frac{2B\mu}{R\_1^2} = 4\pi \mu \frac{(\omega\_2 - \omega\_1)R\_1^2 R\_2^2}{R\_2^2 - R\_1^2}.\tag{3.45}$$

This last relation indicates that we can measure the viscosity μ of a fluid in a Couette viscometer where the drive motor imposes a torque on one of the cylinders and we measure the resulting rotation speed of the other one.

#### *3.2.2 Circular Poiseuille Flow in a Cylindrical Pipe*

Poiseuille flow in a circular pipe with radius *R* is subject to the action of an imposed pressure gradient in direction *z* (Fig. 3.5). The flow is stationary. From the Navier– Stokes equations in cylindrical coordinates, we show first that the only non-zero component of the velocity is v*z*.

**Fig. 3.5** Circular Poiseuille flow in a cylindrical pipe

#### 3.2 Axisymmetric Stationary Flows 61

Given the hypotheses of axial symmetry and stationary flow, vθ = 0 and the only two components of velocity, v*r* and v*z*, are functions only of *r*. The continuity Eq. (A.20) is then

$$\frac{1}{r}\frac{\partial}{\partial r}(r\,\upsilon\_r) = 0\,. \tag{3.46}$$

Integration yields

$$r v\_r = f(z) \, .$$

But, since v*r* = 0 at the wall, *r* = *R*, we conclude that *f* (*z*) = 0 and thus that v*r* is zero everywhere in the flow. The Navier–Stokes equation for the radial component of velocity (A.21) reduces to ∂*p*/∂*r* = 0. The pressure depends only on *z* and not on *r*. The equation for the velocity component v*z* (A.23) yields

$$-\frac{dp}{dz} + \mu \left(\frac{\partial^2 v\_z}{\partial r^2} + \frac{1}{r} \frac{\partial v\_z}{\partial r}\right) = 0$$

or

$$\frac{dp}{dz} = \frac{\mu}{r} \frac{d}{dr} \left( r \frac{dv\_z}{dr} \right) \cdot$$

The left-hand side term only depends on *z*; on the right-hand side there is only dependence on *r*. Thus the two terms must be equal to a constant. Integrating, we obtain

$$v\_z = \left(\frac{dp}{dz}\right)\frac{1}{\mu}\left(\frac{r^2}{4} + A\ln r + B\right).$$

The velocity must be finite on the axis*r* = 0. This leads to *A* ≡ 0. Taking into account the condition v*z*(*R*) = 0, we have

$$v\_{\natural} = -\left(\frac{dp}{dz}\right)\frac{R^2}{4\mu}\left(1 - \left(\frac{r}{R}\right)^2\right) = \frac{GR^2}{4\mu}\left(1 - \frac{r^2}{R^2}\right),\tag{3.47}$$

with the definition *dp*/*dz* = −*G*. In Poiseuille flow, the velocity profile is parabolic. The maximum velocity at the center is

$$v\_{\text{max}} = -\left(\frac{dp}{dz}\right)\frac{R^2}{4\mu} \,. \tag{3.48}$$

Therefore the parabolic profile may be written

62 3 Exact Solutions of the Navier–Stokes Equations

$$
v\_z = v\_{\max} \left( 1 - \frac{r^2}{R^2} \right) \,. \tag{3.49}$$

The flow rate is obtained by integration over the section of the pipe. We have

$$\mathcal{Q} = 2\pi \int\_0^R v\_\varepsilon(r) \, r \, dr = -\left(\frac{dp}{dz}\right) \frac{\pi R^4}{8\mu} = \frac{\pi R^2 v\_{\max}}{2} \,. \tag{3.50}$$

The average, or flux, velocity obtained from the flux divided by the area of the section *S* is

$$
v\_{avg} = \frac{\mathcal{Q}}{\mathcal{S}} = \frac{v\_{\max}}{\mathcal{Q}}\,.\tag{3.51}$$

The maximum velocity is thus equal to twice the average velocity. The shear stress at the cylinder wall, which we denote τw, is given by the component σ*zr* evaluated at *r* = *R*

$$\sigma\_w = -\mu \frac{dv\_\varepsilon}{dr}|\_{r=R} = -\left(\frac{dp}{dz}\right)\frac{R}{2} = \frac{2\mu v\_{\max}}{R} = \frac{4\mu v\_{\arg}}{R} \,. \tag{3.52}$$

The sign change between τw and σ*zr* comes from the fact that τw represents the shear force exercised on the wall by the fluid. The friction coefficient is defined by the ratio of the stress at the wall to the average dynamic pressure

$$C\_f = \frac{\tau\_w}{\frac{\rho v\_{avg}^2}{2}} = \frac{8\mu}{\rho Rv\_{avg}} = \frac{8\nu}{Rv\_{avg}} = \frac{16}{Re\_D} \,, \tag{3.53}$$

with *ReD* being the Reynolds number based on the average velocity and the diameter of the section. It is common to define the head loss coefficient λ by the relation

$$-\left(\frac{dp}{dz}\right) = \frac{\rho v\_{avg}^2}{2}\frac{\lambda}{D}\,. \tag{3.54}$$

Thus it follows that

$$
\lambda = 4C\_f = \frac{64}{Re\_D} \,\mathrm{.}\tag{3.55}
$$

#### *3.2.3 Helical Flow Between Two Circular Cylinders in Relative Motion*

From the geometry point of view, this flow occurs in a configuration similar to the circular Couette flow in Fig. 3.4, with the same notation *R*1 and *R*2 for the radii and ω1, ω2 for the angular rates of rotation. Furthermore the viscous fluid between the cylinders is subjected to an axial pressure gradient. As the flow is steady-state, ∂/∂*t* = 0, and axisymmetric, ∂/∂θ = 0, the velocity profile is a function of *r* only. One has

$$v\_r = v\_r(r), \ v\_\theta = v\_\theta(r), \ v\_z = v\_z(r), \ p = p(r, z) \ . \tag{3.56}$$

With the fluid adherence to the wall, the boundary conditions are

$$\begin{aligned} v\_r(R\_1) = v\_r(R\_2) = 0, \quad &v\_\theta(R\_1) = \omega\_1 R\_1, \ v\_\theta(R\_2) = \omega\_2 R\_2, \\ v\_\varepsilon(R\_1) = v\_\varepsilon(R\_2) = 0 \end{aligned} \tag{3.57}$$

Similarly to the Couette flow solution, it is easy to show that the v*r* component is zero everywhere. The Navier–Stokes equations (A.21)–(A.23) become

$$\frac{1}{\rho} \frac{\partial p}{\partial r} = \frac{v\_{\theta}^{2}}{r} \,, \tag{3.58}$$

$$\frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial v\_{\theta}}{\partial r}\right) - \frac{v\_{\theta}^{2}}{r} = 0\ ,\tag{3.59}$$

$$-\frac{\partial p}{\partial z} + \mu \left( \frac{\partial^2 v\_z}{\partial r^2} + \frac{1}{r} \frac{\partial v\_z}{\partial r} \right) = 0. \tag{3.60}$$

The Couette solution (3.41) remains valid. Equation (3.58) yields

$$p = \rho \int\_{R\_1}^{r} \frac{v\_{\theta}^2}{r'} dr' + f(z) \,, \tag{3.61}$$

where vθ is the Couette profile and *f* (*z*)is an indeterminate function of*z*. Introducing (3.61) in (3.60), one gets

$$-\frac{df}{dz} + \mu \frac{1}{r} \frac{d}{dr} \left( r \frac{dv\_\varepsilon}{dr} \right) = 0 \,\,\,\,\,\tag{3.62}$$

The solutions are obtained taking the boundary conditions (3.57) into account

*f* (*z*) = −*Az* + *B* , (3.63)

$$v\_z(r) = \frac{A}{4\mu} \left[ -r^2 + \frac{R\_2^2 - R\_1^2}{\ln(R\_2/R\_1)} \ln r + \frac{R\_1^2 \ln R\_2 - R\_2^2 \ln R\_1}{\ln(R\_2/R\_1)} \right].\tag{3.64}$$

The pressure field is given by

$$p(r,z) = \rho \int\_{R\_1}^{r} \frac{v\_{\theta}^2}{r'} dr' - Az + B. \tag{3.65}$$

The pressure is known up to a constant *B*, which will give the reference pressure; the pressure gradient −*A* acts in the direction of the axis and finally, the first term of the right hand side member of (3.65) equilibrates the centrifugal forces of the rotating fluid. Note that the axial velocity is independent of the rotation speed of the cylinders, while the azimuthal velocity vθ is independent of the pressure gradient.

#### **3.3 Plane Transient Flows**

Let us now turn our attention to plane flows that depend on time. This situation leads to partial differential equations with independent variables of space and time. In order to arrive at an analytic solution of the problem, we use a change of variables to obtain an ordinary differential equation that is easier to solve.

#### *3.3.1 Transient Flow in a Semi-infinite Space*

An incompressible Newtonian viscous fluid occupies a half space (*x*2 ≥ 0), and is at rest for *t* < 0 (Fig. 3.6). At time *t* = 0, the rigid plane which limits the half space is instantaneously set into motion at the constant velocity *U* in the positive direction of axis *x*1. The motion is two-dimensional such that v3 = 0.

The boundary and initial conditions are given by

$$t < 0, \quad v\_1 = v\_2 = 0, \quad \forall \mathbf{x}\_1, \mathbf{x}\_2 \tag{3.66}$$

$$
\hbar \ge 0, \quad v\_1 = U, v\_2 = 0, \text{ for } \ge 0,\tag{3.67}
$$

$$
v\_1 = v\_2 = 0,\text{ for }\mathbf{x}\_2 \to \infty.\tag{3.68}$$

We assume that v1 and v2 are functions of *x*2 and *t*

$$v\_1 = v\_1(\mathbf{x}\_2, t), \ v\_2 = v\_2(\mathbf{x}\_2, t) \; , \tag{3.69}$$

and that the pressure *p* is a function only of *x*2 (there is no horizontal pressure gradient; the flow is generated entirely by the motion of the moving wall). The conservation of mass becomes

$$\frac{\partial v\_2(\mathbf{x}\_2, t)}{\partial \mathbf{x}\_2} = \mathbf{0} \,. \tag{3.70}$$

**Fig. 3.6** Unsteady flow in an infinite half space

The component v2 only depends on time, and with conditions (3.67) and (3.68), it is identically zero for all *t*. The Navier–Stokes equations become

$$
\rho \frac{\partial v\_1}{\partial t} = \mu \frac{\partial^2 v\_1}{\partial x\_2^2} \,, \tag{3.71}
$$

∂*p* ∂*x*2 = 0 . (3.72)

The pressure *p* is constant.

We can, if we wish, include the effect of gravity in the pressure calculation, by writing

$$\frac{\partial p}{\partial \mathbf{x}\_2} = -\rho \mathbf{g} \mathbf{x}\_2 \,. \tag{3.73}$$

Integration of this relation leads to the calculation of the hydrostatic pressure, where the pressure at a point is equal to the weight of the column of fluid located above that position. The hydrostatic pressure, as its name suggests, does not participate in the dynamics of the flow.

The motion equation (3.71) is a diffusion equation, of the same type as the "heat equation". We can transform this partial differential equation into an ordinary differential equation with a variable change that we obtain from dimensional analysis (cf. Sect. 2.5). Since the problem has no spatial scale other than the variable *x*2 nor time scale other than that of *t* itself, we combine them to form the non-dimensional group (compare with Eq. (2.36))

$$
\eta = \frac{\varkappa\_2}{2\sqrt{\nu t}}.\tag{3.74}
$$

This permits us to obtain an ordinary differential equation for which the solution is a function of η. It is called a self-similar solution because the velocity profile with respect to the variable *x*2 is similar for all times *t*.

Setting

$$v\_1 = U\left|f(\eta)\right|,\tag{3.75}$$

relation (3.71) becomes

$$f'' + 2\eta f' = 0\,,\tag{3.76}$$

with conditions

$$
\eta = 0, \; f = 1; \quad \eta \to \infty, \; f = 0 \; . \tag{3.77}
$$

Integrating (3.76), we obtain

$$f = A \int\_0^\eta e^{-\eta'^2} d\eta' + B \,. \tag{3.78}$$

Taking into account conditions (3.77), we have for η = 0, *B* = 1 and for η = ∞, *A* = −2/ √π where we introduced the error function erf(*x*) defined by Abramowitz and Stegun [1]

$$\text{erf}\left(\mathbf{x}\right) = \frac{2}{\sqrt{\pi}} \int\_0^{\mathbf{x}} e^{-\mathbf{r}^2} d\mathbf{r}\;,\tag{3.79}$$

such that erf(∞) = 1. Then

$$f = 1 - \text{erf}\,\eta\,,\tag{3.80}$$

and from (3.75) the velocity of the fluid for *t* > 0 is

$$v\_1 = U[1 - \text{erf}(\frac{\chi\_2}{2\sqrt{\nu t}})] \,. \tag{3.81}$$

The velocity profile v1/*U* as a function of η is shown in Fig. 3.7. For a fixed value of *t*, the variable η is proportional to *x*2. Then, we can deduce the velocity profile at every instant as a function of the distance from the wall. An interesting question is to know the depth of penetration of the wall motion into the semi-infinite space. More precisely, for a given *t*, what is the distance at which the velocity attains, for example, one per cent of the value of *U*? Examining the function erf, 1 − erf has the value 0.01 for η ∼ 2. Defined as such, the penetration depth δ is given by

$$
\eta\_{\delta} = \frac{\delta}{2\sqrt{\nu t}} \cong 2, \quad \delta \cong 4\sqrt{\nu t} \,, \tag{3.82}
$$

which is proportional to the square root of the kinematic viscosity and time. If the viscosity is very small, the fluid "sticks" less to the wall and it has a weak effect. If *t* tends to infinity, the velocity at every point in the half space goes to *U*.

#### *3.3.2 Flow on an Oscillating Plane*

Consider the flow produced by the periodic horizontal oscillation of a plate in its own plane (Fig. 3.8).

Equation (3.71) still applies, and we must solve it with the boundary conditions

$$
v\_1 = U \cos \omega t \quad \text{for} \quad \mathbf{x}\_2 = \mathbf{0} \;. \tag{3.83}$$

After the initial transient phenomena, the fluid velocity gradually becomes a periodic function of time at the same frequency as the plate oscillation. Here we examine this periodic regime. Assume that solution v1 is of the form

$$v\_1 = \mathfrak{N}\left(e^{i\alpha t}f(\mathbf{x}\_2)\right),\tag{3.84}$$

**Fig. 3.7** Transient flow in an infinite half space

**Fig. 3.8** Unsteady flow on an oscillating plane

where  means the real part of a complex expression. The combination of (3.71) and (3.84) yields

$$i a f = \nu \frac{d^2 f}{d x\_2^2} \dots$$

Recall that *i* 1/2 = *ei*π/4; then the only solution that remains finite as *x*2 → ∞ is

$$f = A \exp\left(-(1+i)(a/2\nu)^{1/2}x\_2\right) \dots$$

The imposition of the boundary condition (3.83) leads to *A* = *U* and the solution becomes

$$\upsilon\_1 = U \exp\left(-(\alpha/2\upsilon)^{1/2}\mathbf{x}\_2\right) \cos\left(\alpha t - (\alpha/2\upsilon)^{1/2}\mathbf{x}\_2\right) \,. \tag{3.85}$$

The velocity profile is a damped harmonic oscillation of amplitude *U e*−*x*2 √ω/2ν in a fluid where a layer at distance *x*2 has a phase lag of *x*2 √ω/2ν with respect to the motion at the wall. Two layers of fluid separated by the distance 2π(2ν/ω)1/2

oscillate in phase. This distance constitutes an estimation of the length of the motion and is called the viscous wave penetration depth. That it increases with viscosity and decreases with frequency is not surprising: if we slowly oscillate a flat plate in a sticky fluid, we expect to drag large masses of fluid along with the plate; on the other hand, if we move the plate rapidly in a fluid of low viscosity, we expect the fluid essentially to ignore the plate, except in a thin boundary layer.

#### *3.3.3 Channel Flow with a Pulsatile Pressure Gradient*

Blood flow in the vascular tree is driven by the pulsating pressure gradient produced by the heart that is acting as a pump. In order to avoid (temporarily) the geometrical complexity of cylindrical coordinates of blood flow in the arteries, we will tackle a simplified version of the problem, namely the plane channel flow under an oscillating pressure gradient.

Recall that the standard plane Poiseuille flow as shown in Fig. 3.2 with a steady constant pressure gradient denoted by *G* gives rise to a parabolic velocity profile. Let us add now an oscillating component characterized by the pulsation ω such that

$$-\frac{1}{\rho} \frac{\partial p}{\partial x\_1} = -G - C \cos \alpha t \,, \tag{3.86}$$

with *C* a constant obtained from experimental data, for example. For the sake of simplicity in the analytical treatment, it is customary to resort to Fourier representation and use the relation

$$-\frac{1}{\rho} \frac{\partial p}{\partial x\_1} = -G - \Re(Ce^{i\alpha t})\,. \tag{3.87}$$

As the steady state oscillating solution is sought for the velocity field, the solution is written as a complex function

$$\upsilon\_1 = \upsilon\_P + \Re\left(\mu(\boldsymbol{\omega}, \boldsymbol{x}\_2)e^{i\boldsymbol{at}}\right). \tag{3.88}$$

where the Poiseuille solution v*P* given by Eq. (3.27) corresponds to the constant pressure gradient.

The Navier–Stokes equations lead to the relation

$$\frac{\partial v\_1}{\partial t} = -\frac{1}{\rho} \frac{\partial p}{\partial x\_1} + \nu \frac{\partial^2 v\_1}{\partial x\_2^2} \,. \tag{3.89}$$

With Eqs. (3.87) and (3.88), Eq. (3.89) gives

$$i\omega u = -C + \nu \frac{\partial^2 u}{\partial x\_2^2} \,. \tag{3.90}$$

#### 3.3 Plane Transient Flows 69

The boundary conditions are

$$
\mu(h) = 0, \quad \frac{\partial \mu}{\partial x\_2}(0) = 0\,. \tag{3.91}
$$

The solution of (3.90) is

$$\mu = \Re \left[ \frac{iC}{\omega} \left( 1 - \frac{\cosh \sqrt{\frac{i\omega}{\nu}} x\_2}{\cosh \sqrt{\frac{i\omega}{\nu}} h} \right) \right]. \tag{3.92}$$

Taking the relation *i* 1/2 = (1 + *i*)/√2 into account, the real part of (3.92) yields the velocity field

$$v\_1 = v\_P - \frac{C}{\omega} \left[ \left( 1 - \frac{f\_1(\omega, \chi\_2)}{f\_3(kh)} \right) \sin \omega t - \frac{f\_2(\omega, \chi\_2)}{f\_3(kh)} \cos \omega t \right],\tag{3.93}$$

where the various notations are defined as follows

$$\begin{aligned} k &= \sqrt{\frac{\omega}{2\nu}}, \\ cc(\mathbf{x}) &= \cos(\mathbf{x})\cosh(\mathbf{x}) \ , \\ ss(\mathbf{x}) &= \sin(\mathbf{x})\sinh(\mathbf{x}) \ , \\ f\_1(\boldsymbol{\omega}, \mathbf{x}\_2) &= cc(k\boldsymbol{\omega}\_2)cc(kh) + ss(k\boldsymbol{\omega}\_2)ss(kh) \ , \\ f\_2(\boldsymbol{\omega}, \mathbf{x}\_2) &= cc(k\boldsymbol{\omega}\_2)ss(kh) - ss(k\boldsymbol{\omega}\_2)cc(kh) \ , \\ f\_3(\boldsymbol{\omega}) &= cc^2(\boldsymbol{\upchi}) + ss^2(\boldsymbol{\upchi}) \ . \end{aligned} \tag{3.94}$$

With ω = 1, in Fig. 3.9, the left part represents the flow for a low frequency case or when the viscous forces are important, i.e. *hk* 1, whereas the right part corresponds to high frequency forcing or to a fluid with very low viscosity. The low frequency solution may be obtained by taking the limit of Eq. (3.93) when *k* → 0. In that limit, *cc*(*x*) → 1, *ss*(*x*) → *x* 2 and therefore, one has

$$v\_1 = v\_P + \frac{Ch^2}{2\nu} \cos\omega t \left(1 - (\frac{\chi\_2}{h})^2\right),\tag{3.95}$$

such that the pulsating term is still a parabola with a modified amplitude. The high frequency case or the equivalent inviscid fluid may be treated with the limit *hk* 1. Then *cc*(*x*) → 1 2 *ex* cos *x* and *ss*(*x*) → 1 2 *ex* sin *x*. The limit solution reads

$$v\_1 = v\_P - \frac{C}{\alpha} \left( \sin \alpha t - \sin(\alpha t - \eta)e^{-\eta} \right) \,, \tag{3.96}$$

where the new variable η is defined as

**Fig. 3.9** Pulsating velocity field with ω = 1; *left*: *k* = 1/ √2; *right*: *k* = 5/ √2

$$
\eta = k(h - x\_2) = \frac{h - x\_2}{\sqrt{2v/\alpha}} \,. \tag{3.97}
$$

Note that the first term of the oscillating part is the response of the inviscid fluid (ν = 0) to the pressure gradient.

#### **3.4 Axisymmetric Transient Flows**

This section treats time dependent flows. The first case is the starting process from rest in a circular Poiseuille flow when the pressure gradient is applied at the initial time. The second example is the pulsating flow in a circular rigid pipe which is somehow connected to blood flow simulation.

#### *3.4.1 Starting Transient Poiseuille Flow*

We examine the transient flow in a circular pipe where the fluid starts from rest to reach the Poiseuille steady parabolic profile (3.47). The only non vanishing velocity component is v*z* and the pressure gradient goes instantaneously at *t* = 0 from zero to the value −*G* everywhere. The dynamic equation is from (A.23)

$$G + \mu \left(\frac{\partial^2 v\_z}{\partial r^2} + \frac{1}{r} \frac{\partial v\_z}{\partial r}\right) = \rho \frac{\partial v\_z}{\partial t} \,, \tag{3.98}$$

with the initial condition

$$v\_z(r,0) = 0, \quad 0 \le r \le R\ ,\tag{3.99}$$

and the boundary condition

$$v\_{\varepsilon}(\mathcal{R}, t) = 0, \forall t \;. \tag{3.100}$$

In order to render (3.98) homogeneous, let us change variables

$$w(r,t) = \frac{G}{4\mu} \left(R^2 - r^2\right) - v\_z(r,t) \,. \tag{3.101}$$

The new variable will be solution of the equation

$$\frac{\partial^2 w}{\partial r^2} + \frac{1}{r} \frac{\partial w}{\partial r} = \frac{1}{\nu} \frac{\partial w}{\partial t} \,, \tag{3.102}$$

with the initial condition

$$w(r,0) = \frac{G}{4\mu} \left(R^2 - r^2\right),\tag{3.103}$$

and the boundary condition

$$w(\mathcal{R}, t) = 0, \quad \forall t \;. \tag{3.104}$$

Through the transient phase, the velocity v*z* will increase till the steady state (3.47) is reached, whereas the transient perturbation w(*r*, *t*) will decay to zero. To solve (3.102), we proceed by separation of variables

$$w(r,t) = f(r)\mathbf{g}(t)\,. \tag{3.105}$$

Substituting in (3.102), one gets

$$\frac{d\mathbf{g}(t)}{dt} + \mathcal{C}\boldsymbol{\nu}\,\mathbf{g}(t) = \mathbf{0}\,\,,\tag{3.106}$$

$$\frac{d^2f}{dr^2} + \frac{1}{r}\frac{df}{dr} + C\;f=0\;\,,\tag{3.107}$$

where *C* is an arbitrary constant. The solution of (3.106) reads

$$\log(t) = B \exp(-C\nu t) \,. \tag{3.108}$$

As w(*r*, *t*) decreases with respect to time, we assume that *C* will involve only positive values so that *C* can be written λ2/*R*2. This notation will ease the next computations, as we will observe. Equation (3.107) then becomes

$$\frac{d^2f}{dr^2} + \frac{1}{r}\frac{df}{dr} + \frac{\lambda^2}{R^2}f = 0\,\,. \tag{3.109}$$

The change of variable λ*r*/*R* = *z* leads (3.109) to the canonical form of the Bessel equation (C.1) whose general solution is given by Eq. (C.2) with the Bessel functions *Jk* and *Yk* . Consequently, the solution of (3.109) is

72 3 Exact Solutions of the Navier–Stokes Equations

$$f = C\_1 J\_0(\frac{\lambda r}{R}) + C\_2 Y\_0(\frac{\lambda r}{R})\,. \tag{3.110}$$

As *Y*0 goes to −∞ when *r* → 0, one concludes that *C*2 ≡ 0 for w to remain finite on the axis. The general solution of (3.102) becomes

$$w(r,t) = C\_3 J\_0(\frac{\lambda r}{R}) \exp(-\frac{\lambda^2}{R^2} \nu t) \,. \tag{3.111}$$

The solution (3.111) verifies condition (3.104) for λ values, denoted λ*n*, given by the zeroes of the Bessel function *J*0

$$J\_0(\lambda\_n) = 0 \,. \tag{3.112}$$

The solution is obtained as

$$w(r,t) = \frac{G}{4\mu} \sum\_{n=1}^{\infty} c\_n J\_0(\frac{\lambda\_n r}{R}) \exp(-\frac{\lambda\_n^2}{R^2} \nu t),\tag{3.113}$$

and the coefficients *cn* are determined by (3.103):

$$R^2 - r^2 = \sum\_{n=1}^{\infty} c\_n J\_0(\frac{\lambda\_n r}{R})\,. \tag{3.114}$$

To solve Eq. (3.114), let us recall the orthogonality properties of Bessel functions as expressed by Lommel integrals

$$\int\_{0}^{1} z \, J\_{n}(\lambda\_{i} z) J\_{n}(\lambda\_{j} z) dz = 0, \quad \lambda\_{i} \neq \lambda\_{j} \,, \tag{3.115}$$

$$\int\_0^1 z \left. J\_n^2(\lambda\_i z) dz \right| = \frac{1}{2} \left[ J\_n'(\lambda\_i) \right]^2. \tag{3.116}$$

Solution of (3.114) is obtained with *z* = *r*/*R* as

$$c\_n = \frac{2\mathcal{R}^2}{[J\_0'(\lambda\_n)]^2} \int\_0^1 (1 - z^2) z \, J\_0(\lambda\_n z) dz \, . \tag{3.117}$$

The evaluation of the two integrals in (3.117) is carried out using successively the recurrence relationships (C.4) and then (C.3) with = 2, *m* = 0. This yields

pipe

$$\int\_{0}^{1} z \, J\_{0}(\lambda\_{n} z) dz = \frac{1}{\lambda\_{n}} J\_{1}(\lambda\_{n}) \,, \tag{3.118}$$

$$\int\_0^1 z^3 \ J\_0(\lambda\_n z) dz = \frac{1}{\lambda\_n^4} \left[ \lambda\_n^3 J\_1(\lambda\_n) + 2\lambda\_n^2 J\_0(\lambda\_n) - 4\lambda\_n J\_1(\lambda\_n) \right], \quad (3.119)$$

$$= \frac{1}{\lambda\_n} J\_1(\lambda\_n) - \frac{4}{\lambda\_n^3} J\_1(\lambda\_n) \ .$$

One finds with the help of the relation [*J* 0(λ*n*)] 2 = [*J*1(λ*n*)] 2: (cf. (C.5) for *m* = 0)

$$c\_n = \frac{8\mathcal{R}^2}{\lambda\_n^3 J\_1(\lambda\_n)}\,. \tag{3.120}$$

Taking Eqs. (3.101), (3.113) and (3.120) into account, the velocity profile is

$$v\_{\varepsilon}(r,t) = \frac{G}{4\mu} \left(R^2 - r^2\right) - \frac{2GR^2}{\mu} \sum\_{n=1}^{\infty} \frac{J\_0(\frac{\lambda\_n r}{R})}{\lambda\_n^3 J\_1(\lambda\_n)} \exp(-\frac{\lambda\_n^2}{R^2} \nu t) \,. \tag{3.121}$$

Figure 3.10 shows the velocity variation with respect to time.

#### *3.4.2 Pulsating Flow in a Circular Pipe*

Let us consider the Poiseuille flow in a circular pipe subject to the action of an oscillating pressure gradient. This problem has been analyzed by Uchida [112] and Womersley [122] and constitutes a modeling of the blood flow in a rigid artery, an assumption far from physiological phenomena as arterial walls deform and move under pressure waves [6, 127]. As the cardiac cycle is time-periodic, the pressure gradient can be represented by a Fourier series (a dozen of modes are sufficient)

$$\frac{\partial p}{\partial z} = \sum\_{k=0}^{\infty} c\_k e^{ik\alpha t} \,. \tag{3.122}$$

The continuous component *c*0 corresponds to the time average of the pressure gradient and produces the Poiseuille profile (3.47). In (3.122), ω is the signal frequency such that ω = 2π/*T* , with *T* the period of the phenomenon. The pressure gradient corresponds to the real part of the complex representation. The flow governing equation is obtained from (A.23)

$$-\frac{\partial p}{\partial z} + \mu \left( \frac{\partial^2 v\_z}{\partial r^2} + \frac{1}{r} \frac{\partial v\_z}{\partial r} \right) = \rho \frac{\partial v\_z}{\partial t} \tag{3.123}$$

and the solution is sought in terms of the Fourier series

$$w\_{\varepsilon}(r,t) = \sum\_{k=0}^{\infty} \hat{w}\_k(r) \, e^{ikat} \,. \tag{3.124}$$

Inserting (3.124) in (3.123), one gets

$$\frac{d^2\hat{w}\_k}{dr^2} + \frac{1}{r}\frac{d\hat{w}\_k}{dr} - \frac{i\alpha k}{\nu}\hat{w}\_k = \frac{c\_k}{\mu} \,. \tag{3.125}$$

The boundary conditions are

$$\frac{d\hat{w}\_k}{dr} \mid\_{r=0} = 0, \quad \hat{w}\_k \mid\_{r=R} = 0 \; . \tag{3.126}$$

The particular solution of (3.125) is quickly found as −*i*ω*k* ν wˆ*k* = *ck* μ .

Introducing the dimensionless variable *z* = *r*/*R*, the homogeneous equation becomes

$$\frac{d^2\hat{w}\_k}{dz^2} + \frac{1}{z}\frac{d\hat{w}\_k}{dz} - \frac{i\mathcal{R}^2\alpha k}{\nu}\hat{w}\_k = 0\,,\tag{3.127}$$

whose solution is given by

$$
\hat{w}\_{kh} = C\_1 J\_0(\alpha \sqrt{k} z \,\text{i}^{3/2}) + C\_2 Y\_0(\alpha \sqrt{k} z \,\text{i}^{3/2}) \,, \tag{3.128}
$$

where the quantity α is the dimensionless Womersley number

$$
\alpha = R \sqrt{\frac{\alpha}{\nu}} \,. \tag{3.129}
$$

The Womersley number is the ratio of the radius to the penetration depth and is a characteristic feature of pulsatile blood flow. Typical values of α in the aorta range from 20 for a human in good health to 8 for a cat. Another way of interpreting the Womersley number consists in estimating the distance from the rigid wall, say δ, where the viscous forces and the inertia are of equal magnitude. Near the wall, viscosity is dominant and a rough estimate of the viscous forces is μ*U*/δ2. Near the symmetry axis, inertia dominates and yields the estimate ρω*U*. Equating the two forces leads to the definition

$$
\delta^2 = \frac{\nu}{\omega} \,. \tag{3.130}
$$

Therefore

$$
\alpha = \frac{R}{\delta} \,. \tag{3.131}
$$

If α is large, the viscous effects are confined to a region very close to the wall. In the centre of the flow, the dynamics will be driven by the equilibrium of inertia and pressure forces, resulting in a velocity profile that will be more blunt than the parabolic profile that comes from the balance of viscous and pressure forces.

The functions *J*0 and *K*0 with the complex argument are Kelvin functions of order zero [1]. As the boundary conditions impose a finite velocity on the axis, it follows that *C*2 ≡ 0 because *Y*0 and *Y* 0 → ∞ when *r* → 0. The velocity profile is

$$\begin{split} v\_z(r,t) &= -\frac{c\_0}{4\mu} R^2 \left( 1 - (\frac{r}{R})^2 \right) \\ &+ \frac{R^2}{\mu} \sum\_{k=1}^{\infty} \mathfrak{N} \left( \frac{ic\_k}{\alpha\_k^2} \left[ 1 - \frac{J\_0(\alpha\_k \frac{r}{R} i^{3/2})}{J\_0(\alpha\_k i^{3/2})} \right] e^{ik\alpha t} \right) . \end{split} \tag{3.132}$$

In this last relation, we have set α*k* = α √*k*. The complex function *J*0(*z i* 3/2) with *z* real ≥ 0 is decomposed in a real part *Ber*(*z*) (real Bessel) and an imaginary part *Bei*(*z*) (imaginary Bessel) such that

$$J\_0(z\,i^{3/2}) = Ber(z) + i\,\,Bei(z)\,, \tag{3.133}$$

and the final velocity profile reads

$$\begin{split} v\_{\varepsilon}(r,t) &= -\frac{c\_{0}}{4\mu}R^{2}\left(1-\left(\frac{r}{R}\right)^{2}\right) \\ &+\frac{R^{2}}{\mu}\sum\_{k=1}^{\infty}\mathfrak{N}\left(\frac{ic\_{k}}{a\_{k}^{2}}\left[1-\frac{Ber(\boldsymbol{\alpha}\_{k}\frac{r}{R})+i\operatorname{Bei}(\boldsymbol{\alpha}\_{k}\frac{r}{R})}{Ber(\boldsymbol{\alpha}\_{k})+i\operatorname{Bei}(\boldsymbol{\alpha}\_{k})}\right]e^{ik\boldsymbol{\alpha}t}\right). \end{split} \tag{3.134}$$

The velocity profile for mode *k* = 1 is shown in Fig. 3.11 for Womersley numbers between 3.34 and 6.67.

**Fig. 3.11** Pulsatile velocity profile for various values of the Womersley number

The flow rate across the pipe section is given by

$$\begin{split} \mathcal{Q} &= 2\pi \int\_0^R v\_z(r, t) \, r \, dr \\ &= -\frac{c\_0}{8\mu} \pi \, R^4 + \frac{\pi R^4}{\mu} \sum\_{n=1}^\infty \Re \left( \frac{ic\_k}{a\_k^2} F(\alpha\_k) e^{ikat} \right), \end{split} \tag{3.135}$$

$$F(\boldsymbol{\alpha}\_k) = 1 - 2J\_1(\boldsymbol{\alpha}\_k \, \mathrm{i}^{3/2}) / \left[ \boldsymbol{\alpha}\_k \, \mathrm{i}^{3/2} J\_0(\boldsymbol{\alpha}\_k \, \mathrm{i}^{3/2}) \right]. \tag{3.136}$$

The cases of slow and rapid pulsation are interesting limit cases.

#### **<IP_ADDRESS> Slow Pulsation**

In this case, α is supposed to be small. The Taylor series of (3.133) (cf. [1]) is

$$\begin{split} J\_0(z\,\,^{3/2}\_{\cdot}) &= 1 - \frac{(z/2)^4}{(2!)^2} + \frac{(z/2)^8}{(4!)^2} - \frac{(z/2)^{12}}{(6!)^2} \cdots \\ &\quad + i\left(\frac{(z/2)^2}{(1!)^2} - \frac{(z/2)^6}{(3!)^2} - \frac{(z/2)^{10}}{(5!)^2} - \cdots\right) .\end{split} \tag{3.137}$$

Taking the dominant terms into account, one obtains successively

$$\frac{J\_0(\alpha\_k \frac{r}{R} i^{3/2})}{J\_0(\alpha\_k i^{3/2})} = \frac{1 + i \frac{(\alpha\_k \frac{r}{R})^2}{4} - \frac{(\alpha\_k \frac{r}{R})^4}{64} + \dotsb}{1 + i \frac{\alpha\_k^2}{4} - \frac{\alpha\_k^4}{64}}.\tag{3.138}$$

Developing the inverse of the denominator, one has

$$(1 + i\frac{\alpha\_k^2}{4} - \frac{\alpha\_k^4}{64})^{-1} = 1 - i\frac{\alpha\_k^2}{4} - \frac{\alpha\_k^4}{16} + \frac{\alpha\_k^4}{64} + \dotsb$$

$$\simeq 1 - \frac{i\alpha\_k^2}{4} - \frac{3\alpha\_k^4}{64} \,. \tag{3.139}$$

Eventually we write

$$\frac{J\_0(\alpha\_k \frac{r}{R} i^{3/2})}{J\_0(\alpha\_k i^{3/2})} = 1 - i \frac{\alpha\_k^2}{4} \left( 1 - (\frac{r}{R})^2 \right) - \frac{\alpha\_k^4}{64} \left( (\frac{r}{R})^4 - 4(\frac{r}{R})^2 + 3 \right) + O(\alpha\_k^6) \,. \tag{3.140}$$

The velocity is now

$$v\_{\varepsilon}(r,t) = -\frac{c\_0}{4\mu}R^2\left(1 - (\frac{r}{R})^2\right) - \frac{R^2}{\mu}\sum\_{k=1}^{\infty} \frac{c\_k}{4}\left(1 - (\frac{r}{R})^2\right)\cos(k\omega t)$$

$$+ \frac{c\_k a\_k^2}{64}\left((\frac{r}{R})^4 - 4(\frac{r}{R})^2 + 3\right)\sin(k\omega t) \,. \tag{3.141}$$

If the continuous Poiseuille component of the velocity is neglected and if we set v*z*,*c* = −*c*1*R*2/(4μ) and α1 = α, the first mode *k* = 1 becomes

$$\frac{v\_z}{v\_{z,c}} = \left(1 - (\frac{r}{R})^2\right) \cos(\omega t) + \frac{a^2}{16} \left((\frac{r}{R})^4 - 4(\frac{r}{R})^2 + 3\right) \sin(\omega t) \,\,. \tag{3.142}$$

The flow comprises a Poiseuille part in phase with the imposed pressure gradient and this is perfectly logical as the pulsation is slow and therefore, the flow has time to adjust itself to that variation. The additional term presents a phase shift of π/2. Figure 3.12 shows the velocity profile for values of α = 1/ √2 and ω*t* = 0,π/4,π/2, 3π/4. Starting from the Poiseuille flow at time ω*t* = 0, at mid period (ω = π/2), the profile is still positive while the corresponding pressure gradient vanishes. The phase difference disappears at mid cycle (ω*t* = π) and the Poiseuille flow is recovered.

**Fig. 3.12** Slow pulsatile flow with α2 = 1/2

#### **<IP_ADDRESS> Rapid Pulsation**

The parameter α*r*/*R* takes large values and the axis (*r* = 0) is excluded from the analysis. For high values of its argument, the asymptotic development of *J*0(*z*) (cf. Eq. (14.144) of Arfken et al. [3]) is such that

$$J\_0(z) = \sqrt{\frac{2}{\pi z}} \cos \left( z - \frac{\pi}{4} \right) + O(|z|^{-1}), \quad \text{with} \quad |\arg(z)| < 2\pi \text{ .} \tag{3.143}$$

Using the relation *i* 3/2 = *ei*3π/4 and *s* = α*kr*/*R*, we perform the next algebraic calculation

$$\begin{split} J\_0(e^{i3\pi/4}s) &= e^{-i3\pi/8} \sqrt{\frac{2}{\pi s}} \cos \left( e^{i3\pi/4}s - \frac{\pi}{4} \right) \\ &= e^{-i3\pi/8} \sqrt{\frac{2}{\pi s}} \left( \cos \left( i\frac{s}{\sqrt{2}} \right) - \left( \frac{s}{\sqrt{2}} + \frac{\pi}{4} \right) \right) \\ &= e^{-i3\pi/8} \sqrt{\frac{2}{\pi s}} \left( \cos \left( i\frac{s}{\sqrt{2}} \right) \cos \left( \frac{s}{\sqrt{2}} + \frac{\pi}{4} \right) + \sin \left( i\frac{s}{\sqrt{2}} \right) \sin \left( \frac{s}{\sqrt{2}} + \frac{\pi}{4} \right) \right) \\ &= e^{-i3\pi/8} \sqrt{\frac{2}{\pi s}} \left( \cosh \left( \frac{s}{\sqrt{2}} \right) \cos \left( \frac{s}{\sqrt{2}} + \frac{\pi}{4} \right) + i \sinh \left( \frac{s}{\sqrt{2}} \right) \sin \left( \frac{s}{\sqrt{2}} + \frac{\pi}{4} \right) \right) \\ &= e^{-i3\pi/8} \sqrt{\frac{2}{\pi s}} \cosh \left( \frac{s}{\sqrt{2}} + i \left( \frac{s}{\sqrt{2}} + \frac{\pi}{4} \right) \right) . \end{split}$$

Neglecting the decaying exponential in cosh as we deal with large values of the argument, we obtain

$$J\_0(e^{i3\pi/4}s) = e^{-i3\pi/8} \sqrt{\frac{1}{2\pi s}} e^{\frac{s}{\sqrt{2}}} e^{i\left(\frac{s}{\sqrt{2}} + \frac{\pi}{4}\right)}\,. \tag{3.144}$$

Consequently one finds

$$\frac{J\_0(\alpha\_k \frac{r}{R} i^{3/2})}{J\_0(\alpha\_k i^{3/2})} \approx \frac{1}{\sqrt{r/R}} e^{-(1+i)\frac{\alpha\_k}{\sqrt{2}}(1-\frac{r}{R})} \,. \tag{3.145}$$

The first mode of the velocity profile in absence of the continuous Poiseuille component yields

$$\frac{v\_{\varepsilon}}{v\_{\varepsilon,c}} = \frac{4}{\alpha^2} \left[ \sin(\omega t) - \frac{1}{\sqrt{r/R}} e^{-\frac{\alpha}{\sqrt{2}} (1 - \frac{r}{R})} \sin\left(\omega t - \frac{\alpha}{\sqrt{2}} (1 - \frac{r}{R})\right) \right] + O\left(\frac{1}{\alpha^4}\right). \tag{3.146}$$

We observe that the resulting flow is in complete phase shift by a factor of π/2 with respect to the pressure gradient. We must again draw the attention of the reader to the assumptions of this approximation which bans the presence of the axis.

#### **3.5 Plane Periodic Solutions**

Many exact solutions of the Navier–Stokes equations are obtained for spatial periodic conditions. In this section we consider a two-dimensional (2D) solution due to Walsh [117].

Let us first proof the following theorem.

**Theorem 3.1** (Walsh) *Let us consider a vector field u in the domain that satisfies*

$$
\nabla^2 \mathfrak{u} = \lambda \,\mathfrak{u},\tag{3.147}
$$

$$\text{div } \mathfrak{u} = 0 \,. \tag{3.148}$$

*Then the velocity v* = *e*νλ*t u satisfies the Navier–Stokes equations (1.73)–(1.74) with a pressure p such that*

$$
\nabla p = -\mathbf{v} \cdot \nabla \mathbf{v} \;. \tag{3.149}
$$

The vector *v* is divergence free as is also *u*. Furthermore,

$$\frac{\partial \boldsymbol{\upsilon}}{\partial t} = \boldsymbol{\upsilon} \boldsymbol{\lambda} \boldsymbol{\upsilon} = \boldsymbol{\upsilon} \boldsymbol{\Delta} \boldsymbol{\upsilon} \,. \tag{3.150}$$

It remains to prove that the nonlinear term is a gradient, i.e. that **∇** × (*v* · **∇***v*) = 0. This amounts to showing that

$$\frac{\partial}{\partial \mathbf{x}\_2} \left( \upsilon\_1 \frac{\partial \upsilon\_1}{\partial \mathbf{x}\_1} + \upsilon\_2 \frac{\partial \upsilon\_1}{\partial \mathbf{x}\_2} \right) = \frac{\partial}{\partial \mathbf{x}\_1} \left( \upsilon\_1 \frac{\partial \upsilon\_2}{\partial \mathbf{x}\_1} + \upsilon\_2 \frac{\partial \upsilon\_2}{\partial \mathbf{x}\_2} \right), \tag{3.151}$$

as **curl ∇** = **0**. This is evident by incompressibility and relation (3.150).

In the 2D case, we resort to the streamfunction ψ, assuming that it is an eigenfunction of the Laplacian with eigenvalue λ. Consequently, *u* = (∂ψ/∂*x*2, −∂ψ/∂*x*1) satisfies (3.147)–(3.148) with the same λ. Therefore, *e*νλ*t* ψ is the streamfunction of the associated Navier–Stokes flow. If we have a periodic domain of size 2π, then the eigenvalues λ are of the form λ = −(*k*2 *x*1 + *k*2 *x*2 ), with *kx*1 and *kx*2 positive integers. For given *kx*1 , *kx*2 , the linearly independent eigenfunctions are

$$\begin{aligned} & \cos(k\_{x\_1} x\_1) \cos(k\_{x\_2} x\_2), \sin(k\_{x\_1} x\_1) \sin(k\_{x\_2} x\_2) \ & \dots \\ & \cos(k\_{x\_1} x\_1) \sin(k\_{x\_2} x\_2), \sin(k\_{x\_1} x\_1) \cos(k\_{x\_2} x\_2) \ & \dots \end{aligned}$$

It is possible to build up complicated geometrical patterns by combination of the eigenfunctions named *n*, *m* eigenfunction by Walsh, with λ = −(*n*2 + *m*2). A theorem in number theory shows that integers of the form *p*2*i* and *p*2*i*+1, where *p* is an integer prime number such that *p* ≡ 1 (mod 4), may be written as sums of squares in exactly *i* + 1 manners. For example, 625 = 54 = 252 = 242 + 72 = 202 + 152. Figure 3.13 displays the streamlines corresponding to ψ = sin(25*x*1) + cos(25*x*2) − sin(24*x*1) cos(7*x*2) + cos(15*x*1) cos(20*x*2) − cos(7*x*1)sin(24*x*2).

**Fig. 3.13** ψ isocontours in the square (0,π/4)2

#### **3.6 Pipe Flow**

Engineering and industrial applications require very often to move various fluids flowing through pipes: oil, water, grain, etc. To model that situation, we will consider a steady state laminar longitudinal flow in a pipe of arbitrary section in order to obtain afterwards the optimal shape of the pipe. The material presented in this section is much influenced by the book by Langlois and Deville [49].

Let *x*3 be the axis of the Cartesian coordinate system parallel to the pipe generatrix. Let us assume that a mechanical device like a pump generates a constant pressure gradient *G* such that −∂*p*/∂*x*3 = *G*. This pressure gradient gives rise to a rectilinear laminar flow along the pipe. Therefore, the velocity profile is

$$v\_1 = v\_2 = 0, \quad v\_3 = v\_3(\mathbf{x}\_1, \mathbf{x}\_2) \,. \tag{3.152}$$

It is easily verified that the continuity Eq. (1.73) is satisfied and also the Navier– Stokes equations (1.74) in directions *x*1 and *x*2. If the body force is constant, the third Navier–Stokes equation yields

$$
\Delta v\_3 + \frac{G}{\mu} = 0 \,, \tag{3.153}
$$

where *G* = *G* + ρ*b*3 with the assumption that *p* = *C* − *G x*3,*C* = *cnst* and is the two-dimensional Laplacian, i.e.

$$
\Delta = \frac{\partial^2}{\partial \mathbf{x}\_1^2} + \frac{\partial^2}{\partial \mathbf{x}\_2^2} \,. \tag{3.154}
$$

Equation (3.153) is a special case of Poisson's equation, a linear, second-order partial differential equation of the elliptic type. We shall look for solutions which satisfy the boundary condition

$$
v\_3(\mathbf{x}\_1, \mathbf{x}\_2) = 0 \quad \text{on } \Gamma \,, \tag{3.155}$$

where **Γ** is the periphery of the pipe cross section.

The boundary value problem represented by (3.153) and (3.155) is equivalent to the Dirichlet problem

$$
\nabla^2 \mu = 0 \,, \tag{3.156}
$$

$$
\mu = f(\mathbf{x}\_1, \mathbf{x}\_2) \quad \text{on } \Gamma. \tag{3.157}
$$

To show the equivalence, we need only set

$$v\_3(\mathbf{x}\_1, \mathbf{x}\_2) = \frac{G}{2\mu} \left[ \mu(\mathbf{x}\_1, \mathbf{x}\_2) - f(\mathbf{x}\_1, \mathbf{x}\_2) \right],\tag{3.158}$$

where *f* is *any* function satisfying

$$
\Delta f = \text{2.}\tag{3.159}
$$

#### *3.6.1 Polynomial solutions*

Many polynomials have constant Laplacian; some of these, when equated to zero, yield the equation of a closed contour. Thus setting v3(*x*1, *x*2) equal to such a polynomial, multiplied by an appropriate constant, will satisfy (3.153) and (3.155), with **Γ** the contour on which the polynomial vanishes. We reject immediately all linear polynomials for two reasons: they vanish only on straight lines, not on closed contours; their Laplacians are not only constant, they vanish. Consideration of quadratic polynomials, however, bears fruit.

#### **<IP_ADDRESS> The Elliptical Pipe**

We note first that the Laplacian of any quadratic polynomial is constant. We next recall that equating any quadratic polynomial to zero yields the equation of a conic section. Only the ellipses are of interest to us, for circles were treated in the Sect. 3.2.2 and none of the other conic sections are closed contours. The algebra is simplified by placing the centroid of the ellipse on the *x*3-axis, and orienting the *x*1- and *x*2-axes along the axes of the ellipse as shown in Fig. 3.14. Thus the most general ellipse can be represented by

$$\left(\frac{\chi\_1}{a}\right)^2 + \left(\frac{\chi\_2}{b}\right)^2 = 1\,\,.\tag{3.160}$$

**Fig. 3.14** Elliptic pipe

#### 3.6 Pipe Flow 83

We find by inspection that a solution to (3.153) which vanishes on the ellipse (3.160) is provided by

$$v\_3 = \frac{G}{2\mu} \frac{a^2 b^2}{a^2 + b^2} \left( 1 - \frac{x\_1^2}{a^2} - \frac{x\_2^2}{b^2} \right) \,. \tag{3.161}$$

The volume flow rate is easily calculated. We have

$$Q = \int \int v\_3 \, d\mathbf{x}\_1 d\mathbf{x}\_2 \,, \tag{3.162}$$

where the integration is carried out over the ellipse. With (3.161),

$$Q = \frac{\pi G}{4\mu} \frac{a^3 b^3}{a^2 + b^2} \,. \tag{3.163}$$

This result may be rewritten

$$\mathcal{Q} = \frac{GA^2}{4\pi\mu} \frac{R}{1+R^2} \,, \tag{3.164}$$

where *A* = π*ab* is the area of the ellipse and *R* = *a*/*b* is the ratio of the semi-axes. We then obtain

$$\frac{\partial \mathcal{Q}}{\partial R} = \frac{GA^2}{4\pi\mu} \frac{1 - R^2}{(1 + R^2)^2} \,. \tag{3.165}$$

Hence, for μ, *G*, *A* all fixed, *Q* has an extremum at *R* = 1, which is easily shown to be a maximum. We thus find that the circular pipe is more efficient than any elliptical pipe, in the sense that the circular pipe produces a greater volume flow for a given pressure gradient than does an elliptical pipe of the same cross-sectional area.

#### **<IP_ADDRESS> The Triangular Pipe**

It seems surprising that we can express in closed form the flow through so unlikely a cross section as an equilateral triangle, but such is the way of polynomial methods. If we place the origin at the intersection of the medians and let the negative *x*1-axis pass through one vertex as in Fig. 3.15, the equation of an equilateral triangle with height 3*a* becomes

$$(\mathbf{x}\_1 - a)(\mathbf{x}\_1 - \sqrt{3}\,\mathbf{x}\_2 + 2a)(\mathbf{x}\_1 + \sqrt{3}\,\mathbf{x}\_2 + 2a) = \mathbf{0} \,\, . \tag{3.166}$$

Note that the two last parentheses in (3.166) are the equations for the two slanted sides of the triangle. Since

#### **Fig. 3.15** Triangular pipe

$$
\Delta \left[ (\mathbf{x}\_1 - a)(\mathbf{x}\_1 - \sqrt{3}\,\mathbf{x}\_2 + 2a)(\mathbf{x}\_1 + \sqrt{3}\,\mathbf{x}\_2 + 2a) \right] = 12a\,,\tag{3.167}
$$

we have

$$v\_3 = \left(\frac{G}{12a\mu}\right)(a - x\_1)(x\_1 - \sqrt{3}x\_2 + 2a)(x\_1 + \sqrt{3}x\_2 + 2a)\,. \tag{3.168}$$

It is readily verified that the maximum velocity occurs at the origin, and that its value there is *Ga*2/3μ.

#### *3.6.2 The Rectangular Pipe*

The flow in a rectangular pipe cannot be solved using the previous approach of polynomial solutions as the equation of a rectangle does not have constant Laplacian. Here we resort to the method of separation of variables.

The rectangular section of the pipe has sides 2*a* and 2*b* and is depicted in Fig. 3.16.

The governing equation (3.153) is subjected to the boundary condition

$$
v\_3 = 0 \quad \text{on} \quad \mathbf{x}\_1 = \pm a \quad \text{and} \quad \mathbf{x}\_2 = \pm b \; . \tag{3.169}$$

The solution of (3.153) is decomposed in a particular solution v3*p* taking the constant term into account and a homogeneous periodic solution v3*h* of the equation

$$
\Delta v\_{3h} = 0 \,. \tag{3.170}
$$

#### **Fig. 3.16** Rectangular pipe

The particular solution is provided by the plane Poiseuille problem (3.27) written as

$$
v\_{3p} = \frac{G}{2\mu} \left(b^2 - \chi\_2^2\right) \,. \tag{3.171}$$

Consequently the homogeneous solution must satisfy the boundary conditions

$$v\_{3h} = \frac{G}{2\mu}(\mathbf{x}\_2^2 - b^2) \quad \text{on } \mathbf{x}\_1 = \pm a \text{ },\tag{3.172}$$

$$v\_{3h} = 0 \quad \text{on } \mathfrak{x}\_2 = \pm b \; . \tag{3.173}$$

The homogeneous solution of (3.170) is obtained by separation of variables such that

$$w\_{3h} = \sum\_{n=0}^{\infty} X\_n(\mathbf{x}\_1) Y\_n(\mathbf{x}\_2) \,. \tag{3.174}$$

One has

$$\frac{X\_{n}^{\prime\prime}}{X\_{n}} = -\frac{Y\_{n}^{\prime\prime}}{Y\_{n}} = k\_{n}^{2} \,. \tag{3.175}$$

The *Yn* solution reads

$$Y\_n(\mathbf{x}\_2) = \cos k\_n \mathbf{x}\_2 \,. \tag{3.176}$$

With the boundary condition (3.173), the coefficient *kn* is obtained as *kn* = (2*n* + 1)π/(2*b*). The corresponding *Xn*(*x*1) is determined from

$$\frac{X\_{n}^{\prime\prime}}{X\_{n}} = k\_{n}^{2}\,,\tag{3.177}$$

i.e.,

$$X\_n = A\_n \cosh k\_n \mathbf{x}\_1 + B\_n \sinh k\_n \mathbf{x}\_1 \ . \tag{3.178}$$

Because of the symmetry of the problem, we have *Bn* = 0. Thus we let

$$w\_{\mathfrak{M}}(\mathbf{x}\_1, \mathbf{x}\_2) = \sum\_{n=0}^{\infty} A\_n \cosh k\_n \mathbf{x}\_1 \cos k\_n \mathbf{x}\_2 \,, \tag{3.179}$$

and set about determining the coefficients *An* so as to satisfy the boundary conditions (3.172). With the symmetry of the hyperbolic cosine, we need

$$\sum\_{n=0}^{\infty} A\_n \cosh[(n+\frac{1}{2})\frac{\pi a}{b}] \cos k\_n x\_2 = \frac{G}{2\mu} (\mathbf{x}\_2^2 - b^2) \,. \tag{3.180}$$

The cosine functions are orthogonal in Fourier analysis. They satisfy the relation

$$\frac{1}{b} \int\_{-b}^{b} \cos k\_n x\_2 \, \cos k\_m x\_2 \, d\mathbf{x}\_2 = \delta\_{mn} \,. \tag{3.181}$$

Multiplying the two sides of (3.180) by cos *km x*2 and integrating, we have

$$A\_m \cosh[(m+\frac{1}{2})\frac{\pi a}{b}] = \frac{G}{2\mu b} \int\_{-b}^{b} (\mathbf{x}\_2^2 - b^2) \cos k\_m \mathbf{x}\_2 \, d\mathbf{x}\_2 \,. \tag{3.182}$$

Utilizing the relation

$$\int x\_2^2 \cos x\_2 \, dx\_2 = 2x\_2 \cos x\_2 + (x\_2^2 - 2) \sin x\_2 \, \, , \tag{3.183}$$

one gets

$$A\_m \cosh[(m+\frac{1}{2})\frac{\pi a}{b}] = \frac{G}{2\mu} \frac{4(-1)^{m+1}}{b k\_m^3} = \frac{G}{2\mu} \frac{32(-1)^{m+1} b^2}{(2m+1)^3 \pi^3} \,. \tag{3.184}$$

The final solution reads

$$v\_3 = \frac{G}{2\mu} \left[ b^2 - x\_2^2 + \frac{32b^2}{\pi^3} \sum\_{n=0}^{\infty} \frac{(-1)^{n+1} \cosh(2n+1) \frac{\pi x\_1}{2b} \cos(2n+1) \frac{\pi x\_2}{2b}}{(2n+1)^3 \cosh(2n+1) \frac{\pi a}{2b}} \right]. \tag{3.185}$$

Figure 3.17 shows isocontour lines for v3/v*max* = *cst* in a channel of rectangular section with *a* = 2*b* = 4.

**Fig. 3.17** v3/v*max* = 0.2, 0.4, 0.6, 0.8, 1.0 isocontours for a rectangular section with *a* = 2*b* = 4

#### **Exercises**

**3.1** Consider the two-dimensional Couette–Poiseuille flow obtained by superimposing Couette flow induced by the constant velocity motion *U* of the upper wall and Poiseuille flow resulting from a pressure gradient in direction *x*1. Calculate the velocity profile, the shear stress, and the flow rate.

**3.2** A solid sphere of radius *R* is immersed in an incompressible viscous Newtonian fluid that fills the space and is at rest at infinity. The sphere rotates about its diameter at a constant angular velocity . Assume that the Reynolds number is less than one and neglect the volume forces. The streamlines are circles centered on the rotation axis in planes perpendicular to this axis. Working in spherical coordinates, calculate the velocity profile.

**3.3** With the same hypotheses as in the preceding exercise, examine the flow of a fluid between two spheres of radii *R*1 and *R*2 such that *R*1 < *R*2, which rotate at the angular velocities 1 and 2 about a common, fixed axis. Calculate the velocity profile. This solution is called the spherical Couette flow.

**3.4** A cylinder of radius *R*1 moves parallel to its axis with a constant velocity *U* inside a fixed, coaxial cylinder of radius *R*2.

Calculate the velocity field of a viscous fluid which fills the space between the two cylinders. Find the friction force per unit length that acts on the moving cylinder.

**Fig. 3.18** Couette flow with free surface

#### **3.5 Couette flow with a free surface**

Consider the circular Couette flow between two cylinders. The viscous fluid fills the annular space under a free surface that is in contact with the ambient air assumed to behave as an inviscid fluid. Figure 3.18 shows the Couette device.

The outer cylinder rotates with the velocity vθ (*R*2) = ω2*R*2, while the inner cylinder rotates at velocity vθ (*R*1) = ω2*R*1.

#### **3.6 Plane Couette flow with two layers**

We consider the plane Couette flow of two incompressible fluids of different viscosities and densities as shown in Fig. 3.19. The steady flow is oriented in direction *x*1 and is generated by the upper wall moving with constant velocity*U*. The two layers of height *h*1 and *h*2, respectively, do not mix and the interface is positioned at *x*2 = 0. The bottom wall is fixed. It is assumed that no pressure gradient influences the flow.

Evaluate the velocity field and the resulting vorticity.

#### **3.7 Bubble dynamics**

A spherical bubble of inviscid gas is contained in an infinite volume of viscous fluid. We suppose that the pressure *pg* of the gas inside the bubble varies with time. As a consequence the radius *R* of the bubble will also vary with time. The varying **Fig. 3.19** Plane Couette flow with two immiscible fluids

bubble generates a velocity field within the fluid which in turn generates a stress field.

Compute the velocity v*r*, the pressure distribution in the fluid, the governing equation for the bubble dynamics and investigate the treatment of the interface.

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Chapter 4 Vorticity and Vortex Kinematics**

In Sect. 1.8, the incompressible Navier–Stokes equations were derived for a viscous Newtonian fluid in terms of primitive variables: velocity and pressure. The phenomena observed in fluid flows have been interpreted by an equilibrium between the inertial forces, the pressure gradient, the volume forces such as gravity, and the viscous forces. In this chapter, we take a different point of view based on the concept of vorticity.

The presence of vorticity in a flow is an indication of the importance of the viscous effects, given that they are generated by viscous stresses. However rotational flows are computed by solving Euler equation, which show the presence of vortices. For example, the flow over a delta wing may be obtained by solving Euler equation. Vorticity can also be generated by a baroclinic mechanism for compressible flow as in Rayleigh–Taylor instability.

Therefore, under certain assumptions, vorticity possesses the following properties:

Thus the vorticity produced on a solid wall introduces the notion of a boundary layer for which we are led to modify certain conclusions coming from the theory of irrotational perfect fluids. In turbulence, flow dynamics is mostly the result of the stretching or shortening of vortex filaments and their deformation.

#### **4.1 Kinematic Considerations**

The velocity gradient tensor *L* can be decomposed into the sum of a symmetric strain rate tensor *d* and an antisymmetric rotation rate tensor *ω*˙ according to Eq. (1.36). The tensor *d* is given by (1.33) and *ω*˙ by (1.35). Recall that the dual vector ˙ , corresponding to the rotation rate tensor, is the rotation rate vector introduced by (1.37).

In fluid mechanics, we classically introduce the vorticity vector *ω*, defined as the curl of the velocity (1.40). To acquire familiarity with the concept of vorticity, we study the flow near a stagnation point at the origin. The velocity components are such that we have, with a constant *C*,

$$v\_1 = C\boldsymbol{x}\_1,\ v\_2 = -C\boldsymbol{x}\_2,\ v\_3 = 0\,. \tag{4.1}$$

We easily calculate that for this flow *ω* = **0**. A flow with zero vorticity is called *irrotational*.

Now consider the plane Poiseuille flow in a channel of height *h*. If the coordinate system has its origin on the lower wall, the velocity profile, (3.19) with definition (3.24), is given by the relation

$$v\_1 = 4v\_{\max} \frac{\chi\_2}{h} (1 - \frac{\chi\_2}{h}) \,, \tag{4.2}$$

with v*max* being the maximum velocity on the centerline of the channel at *x*2 = *h*/2. The only component of the vorticity is ω3. It is perpendicular to the plane of the flow and its value is

$$
\omega\_3 = \varepsilon\_{321} \frac{\partial v\_1}{\partial x\_2} = -\frac{4v\_{\text{max}}}{h} (1 - \frac{2x\_2}{h}) \,. \tag{4.3}
$$

In this case, the absolute value of the vorticity attains a maximum at the two walls and goes to zero on the centerline of the channel.

From these examples we can conclude that the concept of vorticity has no relation to the curvature of the streamlines. In the first case, the streamlines are curved (ψ = *C x*1*x*2), but the vorticity is zero; while in the second example, the streamlines are straight lines and there is finite vorticity.

The Stokes theorem or the curl theorem transforms the flux of the curl of a vector through a surface *S* into the line integral of that vector around the curvy boundary *C* of the surface. The theorem reads

#### 4.1 Kinematic Considerations 93

#### **Fig. 4.1** Vortex tube

**Theorem 4.1** (Stokes theorem)

$$\int\_{S} \boldsymbol{\omega} \cdot \boldsymbol{\mathfrak{n}} \, dS = \oint\_{C} \boldsymbol{\mathfrak{v}} \cdot \boldsymbol{\mathfrak{r}} \, dl \, . \tag{4.4}$$

From the definition of vorticity, (1.40), and the Stokes theorem (4.4), we obtain the identity

$$I(S) = \int\_{S} \boldsymbol{\mathfrak{o}} \cdot \boldsymbol{\mathfrak{n}} \, dS = \int\_{S} \boldsymbol{\mathfrak{curl}} \, \boldsymbol{\mathfrak{v}} \cdot \boldsymbol{\mathfrak{n}} \, dS = \oint \boldsymbol{\mathfrak{v}} \cdot \boldsymbol{\mathfrak{r}} \, dl = \boldsymbol{\Gamma} \,, \tag{4.5}$$

where *I*(*S*) is the intensity of the vortex tube. The curvilinear integral in (4.5) defines the velocity circulation, , along the closed curve *C*, of the unit tangent vector *τ* . It is thus equal to the vorticity vector flux through an arbitrary surface bounded by the curve. In the following, this property will permit us to systematically link the concept of circulation to an interpretation in terms of vorticity.

Recall that a vortex line (Fig. 4.1) is a line tangent at all its points to the vorticity vector, and that a vortex tube is a family of vortex lines circumscribed by a closed curve. The intensity of a vortex tube, for a surface *S* defined by a closed line enclosing the vortex tube, is the flux *I*(*S*) of vorticity through the surface.

**Theorem 4.2** (Helmholtz) (Vorticity properties) *Helmholtz main theorems about vorticity are as follows:*

• *a vortex tube can only end connected to itself or extend to infinity unless it is cut by a wall.*

The proof of these theorems can be found in Panton's book [71].

#### **4.2 Dynamic Vorticity Equation**

#### *4.2.1 General Equation*

The formulation of the equation that governs vorticity dynamics requires the establishment of certain preliminary relations.

The acceleration term *a* can be written as follows:

$$
\mathfrak{a} = \frac{\partial \mathfrak{v}}{\partial t} + \mathfrak{a} \times \mathfrak{v} + \mathfrak{g} \operatorname{rad} \left( \frac{\mathfrak{v} \cdot \mathfrak{v}}{2} \right), \tag{4.6}
$$

that can be proved by induction

$$\begin{split} a\_{i} &= \frac{\partial v\_{i}}{\partial t} + \varepsilon\_{ijk} \boldsymbol{\omega}\_{j} v\_{k} + \frac{\partial}{\partial \boldsymbol{x}\_{i}} \left( \frac{\boldsymbol{v}\_{j} \boldsymbol{v}\_{j}}{2} \right), \\ &= \frac{\partial \boldsymbol{v}\_{i}}{\partial t} + \varepsilon\_{ijk} \boldsymbol{\varepsilon}\_{jlm} \left( \frac{\partial \boldsymbol{v}\_{m}}{\partial \boldsymbol{x}\_{l}} \right) \boldsymbol{v}\_{k} + \boldsymbol{v}\_{j} \frac{\partial \boldsymbol{v}\_{j}}{\partial \boldsymbol{x}\_{i}} \,, \\ &= \frac{\partial \boldsymbol{v}\_{i}}{\partial t} + (\delta\_{kl}\delta\_{im} - \delta\_{km}\delta\_{il}) \left( \frac{\partial \boldsymbol{v}\_{m}}{\partial \boldsymbol{x}\_{l}} \boldsymbol{v}\_{k} \right) + \boldsymbol{v}\_{j} \frac{\partial \boldsymbol{v}\_{j}}{\partial \boldsymbol{x}\_{l}} \,, \end{split}$$

or

$$a\_i = \frac{\partial v\_i}{\partial t} + v\_k \frac{\partial v\_i}{\partial x\_k} \dots$$

This last expression is none other than the definition of acceleration (1.23). Thus, relation

$$\operatorname{curl} a = \frac{D\omega}{Dt} - \omega \cdot \operatorname{grad} v \tag{4.7}$$

is an identity. As can be seen, applying the curl operator to relation (4.6) leads to

$$
\underline{\text{curl}\,a} = \frac{\partial}{\partial t} \underline{\text{curl}\,v} + \underline{\text{curl}(\omega \times \mathbf{v})} + \underline{\text{curl}\,g} \,\underline{\text{grad}\,a} \left(\frac{\mathbf{v} \cdot \mathbf{v}}{2}\right)
$$

$$
\underline{\text{curl}\,a} = \frac{\partial \omega}{\partial t} + \underline{\text{curl}(\omega \times \mathbf{v})} \tag{4.8}
$$

or

$$\operatorname{curl} \boldsymbol{a} = \frac{\partial \boldsymbol{\omega}}{\partial t} + \operatorname{curl} (\boldsymbol{\omega} \times \boldsymbol{\nu}) \,. \tag{4.8}$$

The term **curl**(*ω* × *v*) can be developed as follows:

$$\mathbf{curl}(\boldsymbol{\omega} \times \boldsymbol{\mathfrak{v}}) = \boldsymbol{\mathfrak{v}} \cdot \mathbf{grad}\,\boldsymbol{\omega} - (\nabla \boldsymbol{\mathfrak{v}})\,\boldsymbol{\omega} + \boldsymbol{\omega} \,\mathrm{div}\,\boldsymbol{\mathfrak{v}} - \boldsymbol{\mathfrak{v}} \,\mathrm{div}\,\boldsymbol{\mathfrak{w}}\,\boldsymbol{\mathfrak{z}}\,\tag{4.9}$$

The last term of (4.9) is zero as div **curl** *v* = 0. From (4.8) and (4.9), it follows that

$$\operatorname{curl} a = \frac{D\omega}{Dt} - (\nabla v)\,\omega + \omega \operatorname{div} \,\mathfrak{v}\,\,\omega$$

From the mass conservation equation (1.50), we obtain the relation

$$\operatorname{curl} a = \frac{D\omega}{Dt} - (\nabla v) \text{ } \omega \text{ },$$

which is equivalent to Eq. (4.7) that can be written in the form

$$\frac{D\boldsymbol{\omega}}{Dt} = (\boldsymbol{\omega} \cdot \boldsymbol{\sigma} \boldsymbol{\omega} \boldsymbol{d}) \, \boldsymbol{\upsilon} + \boldsymbol{\text{curl}} \, \boldsymbol{a} \, . \tag{4.10}$$

From the conservation of momentum (1.58), we write

$$\frac{D\boldsymbol{\omega}}{Dt} = (\boldsymbol{\omega} \cdot \boldsymbol{\sigma} \boldsymbol{\sigma} \boldsymbol{\omega}) \, \boldsymbol{\nu} + \boldsymbol{\text{curl}} \, (\boldsymbol{b} + \frac{1}{\rho} \, \text{div} \, \boldsymbol{\sigma}) \, . \tag{4.11}$$

Using the constitutive equation (1.67) in (4.11), we have

$$\frac{D\boldsymbol{\omega}}{Dt} = (\boldsymbol{\omega} \cdot \boldsymbol{\sigma} \boldsymbol{\sigma} \boldsymbol{\omega} d) \, \boldsymbol{v} + \boldsymbol{\epsilon} \boldsymbol{\text{curl}} \, \boldsymbol{b} - \frac{1}{\rho} \boldsymbol{\text{curl}} (\nabla p) + 2 \boldsymbol{\nu} \, \boldsymbol{\text{curl}} (\boldsymbol{\text{div}} \, \boldsymbol{d}) \, \boldsymbol{\dots} \qquad (4.12)$$

If the body force is conservative, it can be derived from a potential χ, as is the case for gravity. Then we write

$$b = -\nabla \chi \,. \tag{4.13}$$

Consequently, **curl** *b* = **0**, and this term disappears from (4.12). We adopt this hypothesis for the rest of the discussion. Furthermore as for the scalar field *p* one has **curl**(∇ *p*) = 0, Eq. (4.12) is simplified as

$$\frac{D\boldsymbol{\omega}}{Dt} = (\boldsymbol{\omega} \cdot \mathbf{grad})\,\boldsymbol{v} + 2\boldsymbol{\nu}\,\mathbf{curl}(\mathrm{div}\,\boldsymbol{d})\,\boldsymbol{\omega}.\tag{4.14}$$

The left-hand side of relation (4.14) contains the material derivative of the vorticity. On the right-hand side, we find two terms that describe the deformation (stretching and shrinking) and the curvature (bending-tilting) of the vortex lines and the viscous diffusion of the vorticity. We notice that in the incompressible case, the vorticity equation does not contain any pressure contribution unlike the compressible case where a baroclinicity term appears (cf. Botsis-Deville [16]).

#### *4.2.2 Physical Interpretation of Vorticity Dynamics for the Incompressible Perfect Fluid*

For an incompressible fluid (ρ = *constant*), that is inviscid (ν = 0), Eq. (4.14) yields

$$\frac{D\boldsymbol{\omega}}{Dt} = (\boldsymbol{\nabla}\boldsymbol{\nu})\,\boldsymbol{\omega} = (\boldsymbol{\omega}\cdot\boldsymbol{\nabla})\boldsymbol{\nu}\,. \tag{4.15}$$

The term

$$(\boldsymbol{\omega} \cdot \nabla)\boldsymbol{v}$$

does not correspond to any term in the Navier–Stokes equations written with the primitive variables, velocity and pressure. Let us examine what that term means from the physical point of view.

In Fig. 4.2, consider two neighboring points Pr and Q on a vortex line. The points P and Q also define a material line of length *dx* = *dx* , and we can show the equality

$$\frac{D(d\mathbf{x}\_i)}{Dt} = dv\_i = \frac{\partial v\_i}{\partial x\_j} d\mathbf{x}\_j \text{ or } \frac{D(d\mathbf{x})}{Dt} = d\mathbf{x} \cdot \mathbf{g} \, \mathbf{rad} \, \mathbf{v} \, . \tag{4.16}$$

This last equation simultaneously expresses the changes in length and direction of a material line element. Comparison of (4.15) and (4.16) shows that the vorticity vector *ω* plays a role analogous to that of the vector *dx*. Thus, *ω* behaves *as if* it were a material line element instantaneously coinciding with a portion of the vortex line. Let δ*v* be the relative velocity of the fluid at Q with respect to P. In relation (4.15), we can make the substitution

**Fig. 4.2** Portion of a vortex

$$(\nabla v)\,\,\omega = \|\,\,\omega\|\,\,\lim\_{\stackrel{P\,\,Q}{\geq 0}\to 0} \frac{\delta v}{P\,\,Q}\,\,.$$

One part of the change of *ω* measured by (4.15) comes from the rigid body rotation of the material line element (from the component of δ*v* normal to *ω*), and the other part is generated by the shrinking or stretching of the elementary line (from the component δ*v* parallel to *ω*). Finally, Eq. (4.15) can be interpreted as follows: the vorticity is transported by the fluid particles, while being oriented and deformed *as if* it were an elementary material vector.

#### *4.2.3 The Vorticity Number*

Truesdell [110] proposed to measure the vorticity in a flow by means of a dimensionless invariant named the vorticity number. This number evaluates the relative strengths of the rotation and stretching and is defined by the relation

$$\mathcal{W} = \sqrt{\frac{\dot{\boldsymbol{\omega}} : \dot{\boldsymbol{\omega}}}{\mathbf{d} : \mathbf{d}}} \,, \tag{4.17}$$

where *ω*˙ and *d* are the vorticity tensor (1.35) and rate of deformation tensor (1.33), respectively. The symbol : defines the scalar product of two tensors. For example, *d* : *d* = *di jdi j* .

The two limit cases are the irrotational motion where *ω*˙ = 0 and W = 0, and the rigid body motion where *d* = 0 with *ω*˙ = 0 and thus W = ∞. All flows with non zero vorticity will be measured through the vorticity number the range of which is in between 0 and infinity.

As an example, let us consider a generalized Poiseuille flow with the following velocity profile

$$v\_1 = v\_2 = 0, \quad v\_3 = v\_3(x\_1, x\_2) \,. \tag{4.18}$$

It is easily computed that W = 1.

#### **4.3 Vorticity Equation for a Viscous Newtonian Fluid**

We assume now that the viscosity μ is invariable. With (1.68), we write

$$\mathbf{div}\,d\mathbf{v}\,d=\frac{1}{2}\left(\mathbf{grad}\,(\mathrm{div}\,\mathbf{v})+\Delta\,\mathbf{v}\right)\tag{4.19}$$

or using the identity for an arbitrary vector *a*

$$
\Delta \mathfrak{a} = \nabla \cdot \nabla \mathfrak{a} = \nabla(\nabla \cdot \mathfrak{a}) - \text{curl } \text{curl } \mathfrak{a} \tag{4.20}
$$

the relation (4.19) becomes

$$\operatorname{div} d = \operatorname{grad} \left( \operatorname{div} \boldsymbol{\nu} \right) - \frac{1}{2} \operatorname{curl} \operatorname{curl} \boldsymbol{\nu} \,. \tag{4.21}$$

Taking the curl of (4.21) leads to

$$\operatorname{curl}(\operatorname{div} d) = -\frac{1}{2} \operatorname{curl}(\operatorname{curl} \omega) \,. \tag{4.22}$$

The vorticity dynamics equation is obtained by combining (4.14) and (4.22):

$$\frac{D\boldsymbol{\omega}}{Dt} = (\boldsymbol{\omega} \cdot \boldsymbol{\sigma} \boldsymbol{\omega} \boldsymbol{d}) \,\, \boldsymbol{v} - \boldsymbol{\nu} \, \boldsymbol{\epsilon} \, \text{curl} \, \boldsymbol{\epsilon} \, \boldsymbol{\omega} \, \, \text{s} \,. \tag{4.23}$$

Equation (4.23) may be rewritten with the help of (4.20) as

$$\frac{D\boldsymbol{\omega}}{Dt} = \left(\boldsymbol{\omega} \cdot \boldsymbol{\sigma} \boldsymbol{\sigma} \boldsymbol{\omega}\right) \boldsymbol{\upnu} + \boldsymbol{\nu} \,\Delta \boldsymbol{\omega} \,. \tag{4.24}$$

#### **Special Case for Two-Dimensional Flow**

For an incompressible two-dimensional flow, Eq. (4.23) becomes, with the notation ω3 = ω, *D*ω

$$\frac{D\omega}{Dt} = \nu \,\,\Delta\omega \,\,\, , \tag{4.25}$$

because in this special case the term (*ω* · *grad v*) is zero since *ω* is orthogonal to the flow plane and thus to *grad v*. We note that Eq. (4.25) is analogous to that for heat conduction, with the kinematic viscosity replacing the thermal diffusivity. We also notice that Eq. (4.25) is satisfied for ω = 0, that is, for an irrotational flow. However, that solution is inadequate. To understand why, we reason by analogy with the heat equation, which also allows an identically zero solution. We know from the study of heat flow, that any non-uniform distribution of temperature at the wall or non-zero heat flux will generate a variable temperature field in the material. Thus the analogy leads us to conclude that, in the case of a viscous fluid, the vorticity that is generated at the walls will diffuse out by shear and then be carried away by the flow. The creation of vorticity at the wall is the result of the shear stress on the wall. To obtain the value of vorticity at the wall, we can resort to the classical method of Green's functions [101].

#### **4.4 Circulation Equation**

In the context of the hypotheses introduced in the previous section, we prove that for a material curve *C*(*t*), along which the circulation of the velocity vector is (*t*), we can write the following:

$$\frac{d\Gamma}{dt} = -\oint\_{C(t)} \nu(\mathbf{curl}\,\mathbf{curl}\,\mathbf{v}) \cdot d\mathbf{x}\,\,. \tag{4.26}$$

This relation expresses the fact that the variation of the circulation along the material curve is due to the viscosity which dampens the motion.

To obtain (4.26), we must first prove that for a material curve *C*(*t*), we have the following identity:

$$\frac{d\Gamma}{dt} = \oint\_{C(t)} \mathbf{a} \cdot d\mathbf{x} \,. \tag{4.27}$$

For that purpose, we can write the equation

$$\frac{d\Gamma}{dt} = \frac{d}{dt} \oint\_{C(t)} \nu\_i d\mathbf{x}\_i = \frac{d}{dt} \oint\_{C\_0} V\_i \frac{\partial \mathbf{x}\_i}{\partial X\_j} dX\_j \ , \ .$$

in which*C*0 denotes the material curve*C*(*t*) at the instant *t* = *t*0 and *Xi* are the associated Lagrangian coordinates. Denoting by *Ai* and *Vi* the Lagrangian representations of acceleration (1.22) and velocity (1.9), we have:

$$\begin{aligned} \frac{d}{dt} \oint\_{C\_0} V\_i \frac{\partial \mathbf{x}\_i}{\partial X\_j} dX\_j &= \oint\_{C\_0} \left( A\_i \frac{\partial \mathbf{x}\_i}{\partial X\_j} + V\_i \frac{\partial V\_i}{\partial X\_j} \right) dX\_j \\ &= \oint\_{C(t)} a\_i d\mathbf{x}\_i + \oint\_{C\_0} \frac{\partial}{\partial X\_j} \left( \frac{V\_i V\_i}{2} \right) dX\_j \dots \end{aligned}$$

The last term of the right-hand side of this equality is zero on a closed curve.

With relation (1.67), which we use in the motion equation (4.21), taking into account the vector identity (4.20) and Eq. (1.74), we can write

$$\mathbf{a} = -\mathbf{grad}\,\chi - \frac{1}{\rho}\mathbf{grad}\,\,p + 2\boldsymbol{\upsilon}\,\mathbf{grad}\,\,(\mathrm{div}\,\boldsymbol{\upsilon}) - \boldsymbol{\upsilon}\,\mathbf{curl}\,\,\mathbf{curl}\,\,\boldsymbol{\upsilon}\,\,.$$

By the conservation of mass, it leads to

$$a = -\operatorname{grad}\left(\frac{p}{\rho} + \chi\right) - \nu \operatorname{curl} \operatorname{curl} \nu \ . \tag{4.28}$$

Inserting (4.28) in (4.27), we obtain (4.26).

#### **4.5 Vorticity Equation for a Perfect Fluid**

For an incompressible, perfect (ν = 0) fluid, the vorticity dynamics equation (4.24) becomes

$$\frac{D\boldsymbol{\omega}}{Dt} = \boldsymbol{\omega} \cdot \boldsymbol{\sigma} \, \text{grad } \boldsymbol{\mathfrak{v}}\,. \tag{4.29}$$

In the two-dimensional case, *ω* is orthogonal to *grad v* and this relation reduces to

$$\frac{D\alpha}{Dt} = 0\,\,. \tag{4.30}$$

From Eq. (4.29), we deduce that, for a perfect incompressible fluid, if the flow is irrotational at an instant, it remains so. In particular, an initially uniform flow will remain irrotational afterwards.

In the case of a perfect fluid, Eq. (4.26) yields Kelvin's theorem (cf. [71], Sect. 13.10)

**Theorem 4.3** (Kelvin theorem) *The circulation of the velocity along a closed material line does not change, for an incompressible perfect fluid*

$$\frac{d\Gamma}{dt} = 0\,. \tag{4.31}$$

#### **4.6 Bernoulli's Equation**

Bernoulli's equation is obtained from the Navier–Stokes equation (1.74) written for perfect fluids (μ = 0). Assume that the volume forces can be derived from a potential (4.13), then

$$\frac{D\mathbf{v}}{Dt} = -\frac{1}{\rho}\nabla\ p - \nabla\chi\ . \tag{4.32}$$

Using the vector identity

$$
\boldsymbol{v} \cdot \nabla \boldsymbol{v} = \boldsymbol{\omega} \times \boldsymbol{v} + \nabla \left( \frac{\boldsymbol{v} \cdot \boldsymbol{v}}{2} \right) \,, \tag{4.33}
$$

in the material derivative of the velocity, we obtain

$$\frac{\partial \mathfrak{v}}{\partial t} = -\mathfrak{w} \times \mathfrak{v} - \frac{1}{\rho} \nabla \ p - \nabla \left( \frac{v^2}{2} + \chi \right) \,. \tag{4.34}$$

We also assume that the flow is irrotational,*ω* = **0**. This assumption is strong, because real fluids produce rotational flows, such as those produced, for example, by the viscous effects near a wall. Thus, Eq. (4.34) can now be written as

#### 4.7 Vorticity Production on a Solid Wall 101

$$\frac{\partial \mathbf{v}}{\partial t} = -\frac{1}{\rho} \nabla \ p - \nabla \left( \frac{v^2}{2} + \chi \right) \,. \tag{4.35}$$

Since the flow is irrotational, the velocity field can be derived from a potential, , such that

$$
\mathfrak{v} = \nabla \Phi \,. \tag{4.36}
$$

Applying the divergence to (4.36) shows that the potential satisfies a Laplace equation

 = 0 . (4.37)

The Euler equation (4.35) yields

$$\nabla \left( \frac{\partial \Phi}{\partial t} + \frac{v^2}{2} + \chi \right) = -\frac{1}{\rho} \nabla \ p. \tag{4.38}$$

As the left-hand side of (4.38) corresponds to the gradient of a scalar function, the same must be the case for the right-hand side. Consequently, Eq. (4.38) becomes with the assumption ρ = cnst

$$\nabla \left( \frac{\partial \Phi}{\partial t} + \frac{v^2}{2} + \chi + \frac{p}{\rho} \right) = 0 \,\,\,\,\,\,\tag{4.39} $$

We integrate this equation to obtain the general form of *Bernoulli's equation:*

$$\frac{\partial \Phi}{\partial t} + \frac{p}{\rho} + \frac{v^2}{2} + \chi = C(t) \,. \tag{4.40}$$

If the flow is stationary, then (4.40) yields the steady state form of Bernoulli's equation

$$\frac{p}{\rho} + \frac{v^2}{2} + \chi = \mathcal{C},\tag{4.41}$$

which, as is suggested by the second term, is an integral of the energy. Therefore, Bernoulli's equation is a first integral of the Euler equation for the case of a stationary, irrotational, perfect fluid.

#### **4.7 Vorticity Production on a Solid Wall**

The presence of a solid wall in a flow generates vorticity. An important fact to report is the direct link between the viscous wall shear stress and the produced vorticity. It is proposed in this section to establish certain properties of the effort exercised by a viscous fluid on a fixed wall. To this end, we can show that for a viscous

**Fig. 4.3** Vorticity generated on a plane wall

incompressible fluid, considered in a point of the fixed wall, we have the following properties:

Let *S* be a plane wall, *n* the unit normal vector to *S* pointing to the outside of the fluid. The system of Cartesian rectangular axes is chosen such that the axes *x*1, *x*2 are located in *S* and the normal *n* is oriented in the negative direction of the axis *x*3, cf. Fig. 4.3.

The continuity equation expressed at the origin of the axes reads

$$
\frac{
\partial \upsilon\_1
}{
\partial \upsilon\_1
} + \frac{
\partial \upsilon\_2
}{
\partial \upsilon\_2
} + \frac{
\partial \upsilon\_3
}{
\partial \upsilon\_3
} = 0
\ .
\tag{4.42}
$$

As the components v1, v2 vanish on the wall, we obtain

$$\frac{\partial v\_1}{\partial x\_1} = \frac{\partial v\_2}{\partial x\_2} = 0 \; . \tag{4.43}$$

The relations (4.42) and (4.43) imply consequently

#### 4.7 Vorticity Production on a Solid Wall 103

∂v3 ∂*x*3 = 0 . (4.44)

Taking these relations into account, we develop the velocity components in Taylor series with respect to the normal direction to the wall. The dominating terms are

$$v\_1 = 0 + \frac{\partial v\_1}{\partial x\_3} \, |\_0 \, x\_3 + \dotsb \, , \tag{4.45}$$

$$v\_2 = 0 + \frac{\partial v\_2}{\partial x\_3} \, |\_0 \, x\_3 + \dotsb \,\tag{4.46}$$

$$v\_3 = 0 + 0 + \frac{\partial^2 v\_3}{\partial x\_3^2} \, |\_0 \, \frac{x\_3^2}{2} + \dotsb \, . \tag{4.47}$$

It is then possible to evaluate the slope of the wall streamline (w.s.), i.e. the streamline obtained when the distance to the wall goes to zero

$$\frac{d\mathbf{x}\_2}{d\mathbf{x}\_1}|\_{w.s.} = \tan\theta = \lim\_{\mathbf{x}\to 0} \frac{\upsilon\_2}{\upsilon\_1} = \frac{\frac{\partial \upsilon\_2}{\partial x\_3}|\_0}{\frac{\partial \upsilon\_1}{\partial x\_3}|\_0}.\tag{4.48}$$

In this relation θ is the angle at the origin of the axes between the wall streamline and the *x*2 axis. This particular streamline is indeed in the wall as the angle that it makes with the planes (*x*3, *x*1) and (*x*3, *x*2), that is the limits of the respective ratios v3/v1 and v3/v2, vanish altogether.

Let us consider now the vorticity lines. On the wall at the origin of the axes, the vorticity components are

$$\begin{aligned} \rho\_1 \mid\_0 &= \frac{\partial v\_3}{\partial x\_2} - \frac{\partial v\_2}{\partial x\_3} = -\frac{\partial v\_2}{\partial x\_3} \mid\_0 \\ \rho\_2 \mid\_0 &= \frac{\partial v\_1}{\partial x\_3} - \frac{\partial v\_3}{\partial x\_1} = \frac{\partial v\_1}{\partial x\_3} \mid\_0 \\ \rho\_3 \mid\_0 &= \frac{\partial v\_2}{\partial x\_1} - \frac{\partial v\_1}{\partial x\_2} = 0 \end{aligned} \tag{4.49}$$

The normal component of the wall vorticity is zero. It results that the vorticity vector is entirely located in the plane of the wall, where it is generated. We can also show that the wall vorticity lines (v.l.) are always orthogonal to the wall streamlines. Indeed, we have the relationship

$$\frac{d\mathbf{x}\_2}{d\mathbf{x}\_1}|\_{v.l.} = \frac{\alpha\_2}{\alpha\_1} = \frac{\frac{\partial v\_1}{\partial x\_3}|\_0}{-\frac{\partial v\_2}{\partial x\_3}|\_0} = -\frac{1}{\frac{dx\_2}{dx\_1}|\_{w.s.}}\tag{4.50}$$

When moving away from the wall, vorticity lines and streamlines do not necessarily remain orthogonal, especially in a three-dimensional flow.

Let us rotate the system of axes around the *x*3 axis in such a way that the angle θ be zero. By (4.48), one has

$$\frac{\partial v\_2}{\partial x\_3} \mid\_{0} = 0 \;. \tag{4.51}$$

The viscous component of the contact force exerted on the wall by the fluid will be denoted *t* v. We evaluate its value by the relations (4.43), (4.44), (4.51) and the Cauchy theorem (1.53)

$$
\mu\_1^v = n\_3 \sigma\_{13} = -\mu \,\frac{\partial v\_1}{\partial x\_3} = -\mu \,\omega\_2 \,. \tag{4.52}
$$

The wall contact viscous force is thus tangential and directly proportional to the wall vorticity. It is possible to generalize (4.52) for an arbitrary system of axes for a wall that is not a plane, cf. Berker [14]. One obtains the relation

$$t = p\,\mathfrak{n} + \mu(\mathfrak{n} \times \mathfrak{o})\,. \tag{4.53}$$

In the same way that the temperature of a wall gives no indication of the amount of energy leaving it, the wall vorticity does not give information on the vorticity intensity which goes in or out the flow. For the sake of comparison, the theory of heat flow establishes that the heat flux through a plane of normal *n* is given by the product *n* **·** *q* where *q* denotes the heat flux vector. By analogy, we define the diffusive vorticity flux for a viscous incompressible fluid by the equation

$$\begin{aligned} \xi &= -\left(\nabla \omega\right) \mathbf{n} \\ \xi\_i &= -n\_j \frac{\partial \omega\_i}{\partial x\_j} \,. \end{aligned} \tag{4.54}$$

The wall value of ζ*i* is computed from Eq. (4.6) and the momentum equation (1.74) with no body force term

$$\frac{\partial \boldsymbol{\upsilon}}{\partial t} + \nabla(\frac{\boldsymbol{\upsilon} \cdot \boldsymbol{\upsilon}}{2} + \frac{p}{\rho}) = -\boldsymbol{\omega} \times \boldsymbol{\upsilon} - \boldsymbol{\upsilon} \, \mathbf{curl} \, \boldsymbol{\omega} \, . \tag{4.55}$$

On a fixed wall the velocity *v* is zero and (4.55) yields

$$
\nabla p = -\mu \,\,\mathbf{curl}\,\boldsymbol{\omega}\,\,.\tag{4.56}
$$

By (4.49), (4.54) and (4.56), one obtains

$$\frac{\partial p}{\partial \mathbf{x}\_1} = \mu \frac{\partial \omega\_2}{\partial \mathbf{x}\_3} = \mu \,\xi\_2,\tag{4.57}$$

$$\frac{\partial p}{\partial \mathbf{x}\_2} = -\mu \frac{\partial \omega\_1}{\partial \mathbf{x}\_3} = -\mu \,\xi\_1. \tag{4.58}$$

We notice that it is necessary to have a pressure gradient along the wall to maintain vorticity production inside the fluid.

The third component of the diffusive vorticity flux is computed by the property **∇** · *ω* = 0. One has

$$\mathfrak{F}\_3 = -n\_3 \frac{\partial a\_3}{\partial x\_3} = -(\frac{\partial a\_1}{\partial x\_1} + \frac{\partial a\_2}{\partial x\_2})\,. \tag{4.59}$$

The necessary pressure gradient for the wall vorticity generation is produced at the flow start-up as the term ρ ∂*v*/∂*t* is the only one during the first instants that can be compensated by −**∇** *p*, because the viscous part of the Navier–Stokes operator comes into play on longer time scales, especially when the Reynolds number is large, cf. Morton [61].

#### **4.8 Flow Behind a Grid**

Kovasznay [45] examines the steady state two-dimensional exact solution of the Navier–Stokes equation for the laminar flow behind a periodic array of cylinders or rods. The velocity field is assumed to be such that v1 = *U* + *u*1, v2 = *u*2, where *U* is the mean velocity in the *x*1 direction. The vorticity equation (4.25) yields

$$\frac{\partial \omega\_3}{\partial t} + (U + \mu\_1) \frac{\partial \omega\_3}{\partial x\_1} + \mu\_2 \frac{\partial \omega\_3}{\partial x\_2} = \nu \Delta \omega\_3 \,. \tag{4.60}$$

Denoting the grid spacing by δ, we define the Reynolds number as *Re* = δ*U*/ν. The dimensionless vorticity becomes ω = ω3 δ/*U*. The other dimensionless variables are *x* = *x*1/δ, *y* = *x*2/δ, τ = *tU*/δ, 1 + *u* = v1/*U*, v = v2/*U*. The governing equation (4.60) becomes

$$\frac{\partial \boldsymbol{\omega}}{\partial \boldsymbol{\tau}} + (\mathbf{l} + \boldsymbol{\mu}) \frac{\partial \boldsymbol{\omega}}{\partial \boldsymbol{\omega}} + \boldsymbol{v} \frac{\partial \boldsymbol{\omega}}{\partial \mathbf{y}} = \frac{1}{Re} \Delta \boldsymbol{\omega} \,. \tag{4.61}$$

As steady state solutions are sought, the term ∂ω/∂τ vanishes. We are left with

$$
\Delta\omega - Re\frac{\partial\omega}{\partial x} - Re\left(\mu \frac{\partial\omega}{\partial x} + \upsilon \frac{\partial\omega}{\partial y}\right) = 0 \,. \tag{4.62}
$$

To build up the analytical solution, the trick consists in finding an expression that cancels the nonlinear term. The streamfunction is introduced to satisfy the continuity equation

$$u = \frac{\partial \Psi}{\partial \mathbf{y}}, \quad v = -\frac{\partial \Psi}{\partial \mathbf{x}}\,, \tag{4.63}$$

and therefore the vorticity is

$$
\omega = -\Delta \psi \,. \tag{4.64}
$$

Taking the *y* periodicity into account, the streamfunction is set up such that

$$
\psi = f(\mathbf{x}) \sin 2\pi \mathbf{y} \,. \tag{4.65}
$$

With (4.65), the nonlinear term of (4.62) gives

$$f'f'' - ff'''' = 0 \,. \tag{4.66}$$

Assuming that none of the derivatives vanish, we write

$$\frac{f'''}{f''} = \frac{f'}{f} \,. \tag{4.67}$$

Integrating (4.67) we obtain

$$f'' = k^2 \, f \,, \tag{4.68}$$

where *k* is a real or complex arbitrary constant. A further integration yields

$$f = Ce^{k\mathbf{x}}\,. \tag{4.69}$$

With the stream function

$$
\psi = Ce^{k\mathbf{x}} \sin 2\pi \mathbf{y} \tag{4.70}
$$

canceling the nonlinear term in (4.62), we have to seek a solution of the equation

$$
\Delta\omega - Re\frac{\partial\omega}{\partial x} = 0 \,\,. \tag{4.71}
$$

Setting

$$
\omega = \mathbf{g}(\mathbf{x}) \sin 2\pi \mathbf{y} \; , \tag{4.72}
$$

we have

$$\text{g}'' - \text{Re}\,\text{g}' - 4\pi^2 \text{g} = 0\,,\tag{4.73}$$

the solution of which is

$$g(\mathbf{x}) = Ae^{\lambda\_{1}\mathbf{x}} + Be^{\lambda\_{2}\mathbf{x}},\tag{4.74}$$

where

$$
\lambda\_{1,2} = \frac{Re}{2} \pm \sqrt{\frac{Re^2}{4} + 4\pi^2} \,. \tag{4.75}
$$

Combining (4.72) and (4.74), the vorticity is

$$
\omega = \left( A e^{\lambda\_{1^\infty}} + B e^{\lambda\_{2^\infty}} \right) \sin 2\pi \text{y } \,, \tag{4.76}
$$

while Eqs. (4.64) and (4.70) give

$$
\omega = \mathcal{C}(4\pi^2 - k^2)e^{kx}\sin 2\pi \text{y} \,\,. \tag{4.77}
$$

Comparison of (4.76) and (4.77) shows that two solutions are possible

$$k = \lambda\_1, \quad A = -Re\lambda\_1 C, \quad B = 0,\tag{4.78}$$

$$k = \lambda\_2, \quad A = 0, \quad B = -Re\lambda\_2 C \,. \tag{4.79}$$

The constant *C* is obtained by fixing the stagnation point at *x* = 0 and therefore *C* = −1/2π. The corresponding streamfunction is

$$
\psi = \mathbf{y} - \frac{1}{2\pi} e^{\lambda\_{2^\text{\text\text\text\text\text\text\text\text\text\texthead}}} \sin 2\pi \mathbf{y} \,\, . \tag{4.80}
$$

With *Re* = 40 and the corresponding λ2, the streamlines are shown in Fig. 4.4, with pairs of eddies generated behind the cylinders. The flow recovers uniformity downstream through the exponential term of the solution.

As the Kovasznay flow incorporates the nonlinear term, it is a good benchmark to test the numerical accuracy and space convergence of computational methods integrating the Navier–Stokes equation.

#### **4.9 Taylor–Green Vortex**

Taylor and Green [104] proposed an idealized model of a three-dimensional vortex field in order to test the dynamics of turbulence.

Orszag [66, 68] modified this model by shifting the origin of axis *x*3 by a factor of 1 2π to obtain a two-dimensional initial velocity field that was simpler to handle by spectral Fourier methods in a cubic box of period 2π. The velocity is given as

$$\begin{aligned} v\_1 &= \cos x\_1 \sin x\_2 \cos x\_3 \\ v\_2 &= -\sin x\_1 \cos x\_2 \cos x\_3 \\ v\_3 &= 0 \end{aligned} \tag{4.81}$$

The corresponding vorticity components are

$$\begin{aligned} \alpha\_1 &= -\sin x\_1 \cos x\_2 \sin x\_3 \\ \alpha\_2 &= -\cos x\_1 \sin x\_2 \sin x\_3 \\ \alpha\_3 &= -2\cos x\_1 \cos x\_2 \cos x\_3 \ . \end{aligned} \tag{4.82}$$

The initial streamlines are planar curves given by cos *x*1 cos *x*2 = *const* in planes *x*3 = *const*. Nonetheless, the velocity that will develop at later times is fully threedimensional. The initial vortex lines are the curves

$$
\sin x\_1 / \sin x\_2 = \text{const}, \ \sin^2 x\_1 \cos x\_3 = \text{const}, \tag{4.83}
$$

so they are twisted and induce a velocity field able to stretch them. The Taylor– Green vortex is perhaps the simplest example of self-induced vortex stretching by a three-dimensional velocity field.

Orszag did not give any detail to obtain (4.83). Therefore, to compute the vorticity lines, we will rely on a paper by Nore et al. [63], especially the appendix "Taylor– Green Clebsch potentials".

The Clebsch potentials allow to decompose the velocity field as follows

$$
\mathfrak{v} = \nabla \varphi + \lambda \nabla \mu \; , \tag{4.84}
$$

i.e. in a potential part (first term of the r.h.s.) and a rotational part (second term of r.h.s.). The rotational part is chosen to be a complex-lamellar field, that is a flow where the velocity field is orthogonal to its own **curl**. Call it *v*(ω). The complexlamellar field may be considered as a potential if it is divided by an integration factor λ such that

$$
\boldsymbol{\sigma}^{(\boldsymbol{\alpha})} = \lambda \nabla \boldsymbol{\mu} \,. \tag{4.85}
$$

The potential ϕ(*x*), λ(*x*), μ(**x**) are named Clebsch variables. Note that the decomposition (4.84) is not unique.

Taking into account the identity**∇** × (λ**∇**μ) = λ(**∇** × **∇**μ) + **∇**λ × **∇**μ, the vorticity of the field (4.84) yields

$$
\boldsymbol{\omega} = \mathbf{curl} \,\boldsymbol{\nu}^{(\boldsymbol{\omega})} = \nabla \boldsymbol{\lambda} \times \nabla \boldsymbol{\mu} \,. \tag{4.86}
$$

Geometrically speaking, the relation (4.86) shows that the surfaces of constant λ and μ are material vorticity surfaces. As λ and μ both contain vortex lines, their intersection describes the vortex lines. The potentials λ and μ must be invariants under vorticity flow dynamics.

Looking for a general invariant I, one must satisfy the transport equation

$$
\boldsymbol{\omega} \cdot \nabla \boldsymbol{I} = 0 \quad \text{or} \quad \boldsymbol{\omega}\_j \frac{\partial \boldsymbol{I}}{\partial \boldsymbol{x}\_j} = \boldsymbol{0} \,. \tag{4.87}
$$

Carrying through the algebra and dividing by cos *x*1 cos *x*2 sin *x*3, one writes

$$
\tan x\_1 \frac{\partial I}{\partial x\_1} + \tan x\_2 \frac{\partial I}{\partial x\_2} + 2 \frac{\frac{\partial I}{\partial x\_3}}{\tan x\_3} = 0 \;. \tag{4.88}
$$

This equation is solved by separation of variables. We divide (4.88) by I = *u*(*x*1)v(*x*2)w(*x*3) to obtain the relation

$$
\tan x\_1 \frac{\partial [\ln u(\mathbf{x}\_1)]}{\partial \mathbf{x}\_1} + \tan x\_2 \frac{\partial [\ln \upsilon(\mathbf{x}\_2)]}{\partial \mathbf{x}\_2} + 2 \frac{\partial [\ln w(\mathbf{x}\_3)]/\partial \mathbf{x}\_3}{\tan x\_3} = \mathbf{0} \,. \tag{4.89}
$$

Each term in (4.89) must be equal to a constant *c*1, *c*2, *c*3 such that

$$c\_1 + c\_2 + 2c\_3 = 0\,\,. \tag{4.90}$$

The first two terms of (4.89) are of the form

$$\tan x\_1 \frac{\partial [\ln \mu(x\_1)]}{\partial x\_1} = c\_1 \ , \tag{4.91}$$

with a general solution given by *u*(*x*1) = const(sin *x*1)*c*1 . The last term yields the equation

$$\frac{\partial [\ln w(\chi\_3)] / \partial \chi\_3}{\tan \chi\_3} = c\_3 \ , \tag{4.92}$$

with the solution given by w(*x*3) = const(cos *x*3)−*c*3 .

Two independent solutions of (4.90) are *c*1 = 1, *c*2 = 0, *c*3 = −1/2 and *c*1 = 0, *c*2 = 1, *c*3 = −1/2. The Clebsch potentials are

$$
\lambda = \sin(\mathbf{x}\_1) \sqrt{\cos(\mathbf{x}\_3)} \tag{4.93}
$$

$$
\mu = \sin(\mathfrak{x}\_2)\sqrt{\cos(\mathfrak{x}\_3)}\,,\tag{4.94}
$$

that give the vorticity lines (4.83).

Unlike the three-dimensional case that has no closed form solution, the twodimensional Taylor–Green vortex is amenable to an analytical solution

$$\begin{aligned} v\_1 &= e^{-2\upsilon t} \sin x\_1 \cos x\_2 \ , \\ v\_2 &= -e^{-2\upsilon t} \cos x\_1 \sin x\_2 \ , \end{aligned} \tag{4.95}$$

$$p = -\frac{1}{4}e^{-4vt} \left(\cos(2\mathbf{x}\_1) + \cos(2\mathbf{x}\_2)\right) \,. \tag{4.96}$$

This solution is used in numerical fluid mechanics to check the accuracy and stability of Navier–Stokes solvers.

#### **Exercises**

**4.1** Compute the vorticity in the circular Couette flow and verify the Stokes theorem (4.5).

**4.2** By applying Bernoulli's theorem for perfect fluids (4.41), show that the velocity of a jet exiting an orifice in a wall at a distance *h* from the free surface of the fluid in a container is

$$
v = \sqrt{2gh}\,.\tag{4.97}$$

#### **4.3 Hill's Vortex**

Hill's vortex [38] in an incompressible fluid represents the limit case where the vorticity is distributed in the volume inside a sphere of radius *R*. Outside the sphere the flow is considered irrotational. For example this situation concerns the constant velocity fall of a water drop in oil. The other limit case corresponds to filament vorticity distribution.

The problem in spherical coordinates has the *ω* vorticity components

$$
\omega\_{\varphi} = \frac{\alpha r \sin \theta}{R}, \quad \omega\_r = \omega\_\theta = 0, \quad \forall r \le R \tag{4.98}
$$

$$
\omega = \mathbf{0} \quad \forall r > R \; , \tag{4.99}
$$

with *C* a constant.

We assume the problem is axisymmetric with vϕ = 0.

#### **4.4 Simplified Vortex Study of the Draining of a Container**

The physical situation concerns the draining of a container whose free surface is located at the altitude *z* = 0 in a cylindrical coordinate system, with the *z* axis oriented downwards. The fluid flow is axisymmetric, irrotational, incompressible with radial and axial velocity components given respectively by

$$v\_r = -\frac{ar}{2}, \quad v\_z = az\ ,\tag{4.100}$$

**Fig. 4.5** Initial velocity distribution of the vortex sheet

$$\frac{a}{2}\rho\_{\natural}r^{2} + \nu r \frac{d\omega\_{\natural}}{dr} = C \,, \tag{4.101}$$

−*U*

with *C* a constant.

#### **4.5 Diffusion of a vortex sheet**

At the initial time *t* = 0, a vortex sheet coinciding with the plane *x*2 = 0 of the Cartesian coordinate system is given. This vortex sheet is a singularity that we will ignore. On each side of the sheet, the velocity distribution is uniform and such that (cf. Fig. 4.5)

$$\begin{aligned} v\_1(\mathbf{x}\_2, \mathbf{0}) &= \quad U \quad \forall \mathbf{x}\_2 > \mathbf{0} \\ v\_1(\mathbf{x}\_2, \mathbf{0}) &= -U \quad \forall \mathbf{x}\_2 < \mathbf{0} \end{aligned} \tag{4.102}$$

The fluid is incompressible. The body forces are neglected. The pressure is uniform. The flow is one-dimensional in the plane *Ox*1*x*2.

*x*1

Compute the time evolution of the vorticity resulting from the effects of viscous diffusion. To this end, follow the next steps.

• Show that the dynamic vorticity equation reduces to

$$\frac{\partial \alpha\_3}{\partial t} = \nu \frac{\partial^2 \alpha\_3}{\partial x\_2^2} \,. \tag{4.103}$$

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Chapter 5 Stokes Flow**

In this chapter, let us consider first the Stokes equation, valid for very slow flows that we conventionally call creeping flows or Stokes flows. These flows are dominated by viscous forces which are much larger than inertial forces. Examples come from technologies in such diverse domains as convection currents in high temperature glass melting furnaces, lubricants in bearings, and the flow of oils and mud (although the latter may have pronounced non-Newtonian behavior). In nature (another source of interesting cases), we find convection in terrestrial magma, the flow of lava, the swimming of fish, the propulsion of microorganisms, and the squirming of spermatozoon.

We assume that the Reynolds number *Re* - 1 and therefore the Navier–Stokes equations reduce to the Stokes equation. As the latter is linear, a complete analytic treatment is possible.

Taking the divergence of the Stokes equation (2.54) and taking into account the solenoidal character of the velocity field, we get

$$
\Delta p = 0 \,. \tag{5.1}
$$

The pressure is thus a harmonic function for a Stokes flow.

Taking the curl of the Stokes equation (2.54), we obtain

$$\frac{\partial \boldsymbol{\omega}}{\partial t} = \boldsymbol{\nu} \Delta \boldsymbol{\omega} \,, \tag{5.2}$$

where we have introduced the vorticity, *ω* = **curl** *v* using Eq. (1.40). If the flow is stationary, the components of the vorticity are also harmonic functions.

#### **5.1 Plane Creeping Flows**

Consider a plane flow for which we have

$$\mathfrak{v} = (v\_1(\mathbf{x}\_1, \mathbf{x}\_2, t), v\_2(\mathbf{x}\_1, \mathbf{x}\_2, t), \, \mathbf{0}); \quad p = p(\mathbf{x}\_1, \mathbf{x}\_2, t) \,. \tag{5.3}$$

In such a two-dimensional problem, incompressibility (1.73) is automatically satisfied by the introduction of a *stream function* ψ so that

$$
v\_1 = \frac{\partial \psi}{\partial x\_2}, \quad v\_2 = -\frac{\partial \psi}{\partial x\_1}.\tag{5.4}$$

As the vorticity reduces to a single component *ω* = (0, 0, ω), it follows that

$$
\omega = \frac{\partial \upsilon\_2}{\partial \mathbf{x}\_1} - \frac{\partial \upsilon\_1}{\partial \mathbf{x}\_2} = -\Delta \psi \; , \tag{5.5}
$$

and relation (5.2) becomes

$$\frac{\partial \Delta \psi}{\partial t} = \nu \Delta \Delta \psi \,. \tag{5.6}$$

For a stationary problem, we will have

$$
\Delta\Delta\psi = 0\,,\tag{5.7}
$$

showing that in this case the stream function is a biharmonic function.

In polar coordinates (*r*,θ), the conservation of mass becomes

$$\frac{1}{r}\frac{\partial}{\partial r}(rv\_r) + \frac{1}{r}\frac{\partial v\_\theta}{\partial \theta} = 0 \; . \tag{5.8}$$

A stream function ψ also exists such that

$$v\_r = \frac{1}{r} \frac{\partial \Psi}{\partial \theta}, \quad v\_\theta = -\frac{\partial \Psi}{\partial r} \,. \tag{5.9}$$

#### *5.1.1 Flow in a Corner*

Let us consider the flow in a corner as presented in Fig. 5.1. The lower wall is fixed while the wall inclined at an angle θ0 is in uniform translational motion at the constant velocity *U* in the direction *x*1. Near the origin, the velocity gradients are large; nonetheless, we expect the viscous forces to dominate in the neighborhood of the origin. To formulate the problem in steady state, we choose a coordinate system with the origin at the intersection of the two walls, in motion with the inclined wall. In this case, the boundary conditions are written as

$$\frac{1}{r}\frac{\partial\Psi}{\partial\theta} = -U, \quad \frac{\partial\Psi}{\partial r} = 0 \quad \text{at} \quad \theta = 0 \tag{5.10}$$

$$\frac{1}{r}\frac{\partial\psi}{\partial\theta} = 0, \quad \frac{\partial\psi}{\partial r} = 0 \quad \text{at} \quad \theta = \theta\_0 \,. \tag{5.11}$$

**Fig. 5.1** Flow in a corner of angle θ0. The coordinate system moves at the upper wall velocity *U*

The form of the boundary conditions suggests that we can write ψ in the following form:

$$
\psi = r \, f(\theta) \,. \tag{5.12}
$$

Substituting (5.12) in the biharmonic equation

$$\nabla^4 \psi = \left(\frac{1}{r}\frac{\partial}{\partial r} + \frac{\partial^2}{\partial r^2} + \frac{1}{r^2}\frac{\partial^2}{\partial \theta^2}\right) \left(\frac{1}{r}\frac{\partial \psi}{\partial r} + \frac{\partial^2 \psi}{\partial r^2} + \frac{1}{r^2}\frac{\partial^2 \psi}{\partial \theta^2}\right) = 0\,,\tag{5.13}$$

we find the relation

$$\frac{1}{r^3} \left( \frac{d^4 f}{d\theta^4} + 2 \frac{d^2 f}{d\theta^2} + f \right) = 0 \,, \tag{5.14}$$

for which the solution is

$$f(\theta) = A\sin\theta + B\cos\theta + C\theta\sin\theta + D\theta\cos\theta \,. \tag{5.15}$$

(Recall that if *H* is a harmonic function, θ *H* is a biharmonic function). The imposition of the boundary conditions (5.10) and (5.11) allows us to evaluate the constants which are

$$A, B, C, D = \left( -\theta\_0^2, 0, \theta\_0 - \sin\theta\_0 \cos\theta\_0, \sin^2\theta\_0 \right) \frac{U}{\theta\_0^2 - \sin^2\theta\_0}. \tag{5.16}$$

For the special case of a right angle, we have

$$\psi = \frac{rU}{(\frac{\pi}{2})^2 - 1} \left( - (\frac{\pi}{2})^2 \sin \theta + \frac{\pi}{2} \theta \sin \theta + \theta \cos \theta \right) \,, \tag{5.17}$$

from which we can easily obtain the velocity components

$$v\_r = \frac{U}{(\frac{\pi}{2})^2 - 1} \left( (1 - \frac{\pi^2}{4}) \cos \theta + \frac{\pi}{2} (\sin \theta + \theta \cos \theta) - \theta \cos \theta \right) \quad (5.18)$$

$$v\_{\theta} = -\frac{U}{(\frac{\pi}{2})^2 - 1} \left( - (\frac{\pi}{2})^2 \sin \theta + \frac{\pi}{2} \theta \sin \theta + \theta \cos \theta \right) \,. \tag{5.19}$$

In retrospect we can examine the correctness of the creeping flow assumption. We see that the acceleration components (A.11) and (A.12) evaluated with the preceding solution are proportional to *U*2/*r* with a factor that depends on θ, that is, of the order of unity. As for the viscous effects, they are of the order of μ*U*/*r* 2. Thus the creeping flow assumption is met when ρ*rU*/μ - 1. This is the case in the region close to the origin such that *r* ν*U*. Further away, the solution will no longer be correct as the inertial forces rapidly become of the same order of magnitude as the viscous forces.

It is interesting to compute the pressure from Eq. (A.22):

$$\frac{\partial p}{\partial \theta} = r\mu \left( \Delta v\_{\theta} + \frac{2}{r^2} \frac{\partial v\_r}{\partial \theta} - \frac{v\_{\theta}}{r^2} \right) = r\mu \left( -\Delta \frac{\partial \Psi}{\partial r} + \frac{2}{r^3} \frac{\partial^2 \Psi}{\partial \theta^2} + \frac{1}{r^2} \frac{\partial \Psi}{\partial r} \right), \tag{5.20}$$

which yields

$$p = \frac{2\mu}{r(\theta\_0^2 - \sin^2\theta\_0)} \left( (\frac{1}{2}\sin 2\theta\_0 - \theta\_0)\sin\theta - \sin^2\theta\_0\cos\theta \right) \,. \tag{5.21}$$

Observe that the pressure varying like *r*−1 becomes unbounded when we approach the corner. This dismal performance comes from the fact that at the corner the boundary conditions are not consistent with the real problem, which has always a tiny gap so that the forces remain finite.

#### **5.2 Two-Dimensional Corner Moffatt Eddies**

Plane Stokes flows occur in engineering or physical problems in the neighborhood of slots or cracks in a wall. This situation is modeled by analyzing the creeping flow close to the vertex of a sharp wedge with an aperture angle 2α formed by the walls θ = ±α as shown in Fig. 5.2.

The forcing mechanism generating the corner flow is "far" from the wedge vertex. For example, in the lid-driven square cavity problem, corner vortices are generated by the influence of the main primary vortex.

For the flow in a wedge, two geometrical configurations are possible: the flow is asymmetric with respect to the symmetry axis θ = 0; the flow is symmetric in θ and the symmetry axis plays the role of a mirror.

Here we will concentrate on the asymmetric case shown in Fig. 5.2 with the assumption of 2D creeping flow described by the stream function. The governing equation is still (5.7) which in polar coordinates is given by

$$
\Delta\Delta\psi = \left(\frac{\partial^2}{\partial r^2} + \frac{1}{r}\frac{\partial}{\partial r} + \frac{1}{r^2}\frac{\partial^2}{\partial \theta^2}\right)^2 \psi\ . \tag{5.22}
$$

We look for separable solutions of (5.22) expressed by

$$\psi = r^{\lambda} f(\theta) \,, \tag{5.23}$$

**Fig. 5.2** Geometry of the corner flow

where λ is a real or complex number. Carrying through the algebra, the biharmonic equation yields

For λ = 0, 1, 2, the solution of (5.24) is given by

$$\psi = r^{\lambda} (A e^{i\lambda \theta} + B e^{i(\lambda - 2)\theta}) \,. \tag{5.25}$$

Following Moffatt [60], the two-dimensional streamfunction is expanded in a series of basic solutions ψ = *r* λ *f* (θ )

$$\psi = \mathfrak{N} \sum\_{n=1}^{\infty} A\_n r^{\lambda\_n} f\_{\lambda\_n}(\theta) \,, \tag{5.26}$$

where the *An* are complex numbers and the λ*n* satisfy the condition

$$1 < \Re \lambda\_1 < \Re \lambda\_2 < \cdots \tag{5.27}$$

The first inequality imposes that the flow vanish at the origin, which is located at the corner. The remaining inequalities indicate that the first term in (5.26) dominates over the others and then

$$
\psi \approx A\_1 r^{\lambda\_{\parallel}} f\_{\lambda\_{\parallel}}(\theta) = A r^{\lambda} f\_{\lambda}(\theta) \,. \tag{5.28}
$$

As for the asymmetric solution v*r*(*r*, −θ ) = −v*r*(*r*,θ) and vθ (*r*, −θ ) = vθ (*r*,θ), the function *f*λ has to be even in θ. The term in between parentheses in solution (5.25) is rewritten as

$$A\sin\lambda\theta + B\cos\lambda\theta + C\cos(\lambda - 2)\theta + D\sin(\lambda - 2)\theta \text{ , }\tag{5.29}$$

and is such that the constants *A* and *D* vanish:

$$f\_{\lambda} = B \cos \lambda \theta + C \cos(\lambda - 2)\theta \text{ .} \tag{5.30}$$

The no slip boundary conditions on the two walls *f*λ(±α) = *f* λ(±α) = 0 yield the system

$$B\cos\lambda\alpha + C\cos(\lambda - 2)\alpha = 0\ ,\tag{5.31}$$

$$B\lambda\sin\lambda\alpha + C(\lambda - 2)\sin(\lambda - 2)\alpha = 0\ . \tag{5.32}$$

For a nonzero solution the determinant of this system must vanish, i.e.

$$
\sin 2(\lambda - 1)\alpha + (\lambda - 1)\sin 2\alpha = 0 \,\,. \tag{5.33}
$$

#### *5.2.1 Real Solutions for λ (α >* **73***.***15◦***)*

The nonlinear equation (5.33) gives real solutions for an angle α > 73.15◦. Figure 5.3 shows the real solutions obtained as the intersections of the sine function sin 2(λ − 1)α in black and the straight line −(λ − 1)sin 2α in green with respect to the variable (λ − 1)α. The smallest value is the relevant one when we approach the corner*r* → 0 as the solution goes like *r* λ.

#### *5.2.2 Complex Solutions for λ (α <* **73***.***15◦***)*

Let us write λ = *p* + 1 + *iq*. The azimuthal velocity component is

$$v\_{\theta}(r,\theta) = -\frac{\partial \psi}{\partial r} = \mathfrak{R}\left(-\lambda r^{\lambda - 1} \, f(\theta)\right) \,. \tag{5.34}$$

On the symmetry axis of the corner, θ = 0, and therefore

$$v\_{\theta}(r,0) = \Re\left(-\lambda r^{\lambda - 1} \, f(0)\right) = \Re\left(r^{\lambda - 1}C\right) \, , \tag{5.35}$$

**Fig. 5.3** Real solutions for α = 80◦

where *C* = |*C*|*ei*β ≡ −λ*f* (0). Equation (5.35) yields

$$\upsilon\_{\theta}(r,0) = \Re\left(r^{p}|C|e^{iq\ln r}e^{i\beta}\right) = r^{p}|C|\cos(q\ln r + \beta)\,. \tag{5.36}$$

When *r* → 0, ln *r* → −∞ and the velocity vθ (*r*, 0) changes sign infinitely often. This behavior means that a string of counter-rotating vortices is present in the corner. The center of the *n*th corner eddy denoted by *rn* is the distance of this center to the origin. It is given by the relation vθ (*r*, 0) = 0 leading to

$$q \ln r\_n + \beta = -(2n+1)\frac{\pi}{2}, \ n = 0, 1, 2, \dots, \quad \text{or} \quad r\_n = e^{-(2n+1)\frac{\pi}{2q}}e^{-\frac{\beta}{q}}.\tag{5.37}$$

A simple calculation yields

$$\frac{r\_n}{r\_{n+1}} = \frac{r\_n - r\_{n+1}}{r\_{n+1} - r\_{n+2}} = e^{\frac{\pi}{q}}\,,\tag{5.38}$$

which shows that the sizes of the vortices fall off in geometrical progression with a common ratio *e*π/*q* , depending on the aperture angle of the corner. If we now inspect the velocity maxima, we find vθ ,*max* = *r p*|*C*| at points *rn*+ 1 2 = *e* −*n* π 2*q e* − β *q* . The maximum velocity will be called the intensity of the eddy. The corner vortices have their intensities falling off in geometrical progression with the common ratio

$$\left|\frac{v\_n}{v\_{n+1}}\right| = \left|\frac{r\_{n+\frac{1}{2}}}{r\_{n+\frac{3}{2}}}\right|^p = e^{\pi p/q} \,, \tag{5.39}$$

which also depends on α.

The numerical solution of Eq. (5.33) for the angle α is given in Table 5.1.

We observe that λ decreases when α increases. Furthermore, the imaginary part λ goes to zero when α reaches the value 73.15◦, meaning that λ then becomes real.

**Table 5.1** Main eigenvalue λ with respect to the corner half angle

**Fig. 5.4** Corner flow in a wedge of aperture α = 28.5◦. (Courtesy of E. Rønquist [81]. Taneda picture [102] (c) (1979) The Physics Society of Japan, is reprinted with permission)

Figure 5.4 shows a spectral element solution computed by Einar Rønquist [81] for the Stokes flow in a wedge. As the top lid moves at unit velocity, a series of Moffat corner eddies is generated. These eddies are stacked from top to bottom in an infinite cascade to the tip. The wedge shown has an aperture angle of 28.5◦. The asymptotic ratio of successive eddy intensities is 405. With the discretization comprising 30 elements of polynomial degree 8 shown in the figure, one obtains four eddies. The ratio of the strength of two successive eddies is from top to bottom 386, 406, and 411. We observe also that the computed results are very close to the experimental data provided by Taneda [102] and Van Dyke [114].

#### **5.3 Stokes Eigenmodes**

Most state-of-the-art numerical methods dealing with the Navier–Stokes equation rely on an implicit treatment of the Stokes operator and an explicit scheme for the nonlinearity (rejecting it as a source term). Hence it is essential to understand the structure of the Stokes operator. Furthermore if we can write the eigenmodes in closed form, or if we can compute them accurately, then we are able to use those modes as the basis for the approximation of the Stokes equation, see for example Batcho-Karniadakis [12].

The Stokes eigenproblem is defined by setting ∂*v* ∂*t* = λ*v* in Eq. (2.54) and assuming that *b* = **0**. The eigenvalue λ provides the growth or decay rate of the velocity field. Now the eigensystem becomes

$$
\lambda \boldsymbol{\nu} - \nu \Delta \boldsymbol{\nu} + \nabla p = \mathbf{0} \quad \text{in} \quad \Omega \text{ }, \tag{5.40}
$$

$$\text{div } \mathfrak{v} = 0 \qquad \text{in} \quad \Omega \text{ ,} \tag{5.41}$$

*v* = **0** on ∂
 , (5.42)

where *p* is normalized by ρ, the problem domain and ∂ its boundary.

#### *5.3.1 Periodic Stokes Eigenmodes*

Let us consider the fully periodic solutions of the transient Stokes problem Eqs. (2.54) in the open square domain =] − 1, +1[ 2 with the Fourier approximation

$$\begin{aligned} \boldsymbol{\nu}(\mathbf{x},t) &= \sum\_{\|\mathbf{k}\|<\infty} \widehat{\boldsymbol{\nu}}(\mathbf{k}) e^{i(\mathbf{k}\cdot\mathbf{x}+\boldsymbol{\alpha}t)}, & \boldsymbol{p} &= \sum\_{\|\mathbf{k}\|<\infty} \widehat{\boldsymbol{p}}(\mathbf{k}) e^{i(\mathbf{k}\cdot\mathbf{x}+\boldsymbol{\alpha}t)}, \\ \boldsymbol{b}(\mathbf{x},t) &= \sum\_{\|\mathbf{k}\|<\infty} \widehat{\boldsymbol{b}}(\mathbf{k}) e^{i(\mathbf{k}\cdot\mathbf{x}+\boldsymbol{\alpha}t)}, & \boldsymbol{b} &= \boldsymbol{b} \end{aligned} \tag{5.43}$$

where *v*, *p*, *b* are the complex Fourier coefficients, *k* the wavevector and ω a complex frequency. The notation *k* < ∞ is defined to mean −∞ < *ki* < +∞ for *i* = 1, 2.

Let us denote *ek* the unit vector in the direction of *k*. Then *k* = *k ek*. The solution is easily obtained

$$
\widehat{p} = -i\frac{\mathbf{e}\_k \cdot \mathbf{b}}{k}, \quad (i\omega + \nu k^2)\widehat{\mathbf{v}} = \left(\widehat{\mathbf{b}} - \mathbf{e}\_k(\mathbf{e}\_k \cdot \widehat{\mathbf{b}})\right). \tag{5.44}
$$

The resulting periodic eigenmode corresponding to λ = *i* ω is

$$
\widehat{p} = 0, \quad (\lambda + \nu k^2)\widehat{\mathfrak{v}} = \mathbf{0}, \quad k \cdot \widehat{\mathfrak{v}} = 0 \; . \tag{5.45}
$$

The periodic Stokes modes are constant pressure modes driven only by diffusion as λ = −ν*k*2. The incompressibility constraint div *v* = 0 does not influence the space configuration, except that the wavevector *k* must be orthogonal to the velocity. Geometrically speaking the velocity is contained in a plane perpendicular to *k*, while the pressure is aligned with the wavevector.

#### *5.3.2 Channel Flow Stokes Eigenmodes*

The problem is based on the plane channel flow between horizontal plates as treated in [69]. The flow is assumed periodic in the *x*2 direction while it is confined by rigid walls in *x*1 = ±1. We seek a solution of the two-dimensional Stokes equation in the form

$$\boldsymbol{\sigma}(\mathbf{x},t) = \left(\mu(\mathbf{x}\_{\mathrm{l}})e^{ik\boldsymbol{x}\_{2}+\lambda t}, \boldsymbol{\upsilon}(\mathbf{x}\_{\mathrm{l}})e^{ik\boldsymbol{x}\_{2}+\lambda t}\right), \quad \boldsymbol{p} = p(\mathbf{x}\_{\mathrm{l}})e^{ik\boldsymbol{x}\_{2}+\lambda t},\tag{5.46}$$

where *k* is a chosen wavenumber and *u*, v are complex functions. The Stokes equations satisfied by *u*, v, *p* are

#### 5.3 Stokes Eigenmodes 123

$$\begin{aligned} \lambda u &= -\frac{dp}{dx\_1} + \nu (\frac{d^2 u}{dx\_1^2} - k^2 u) \ , \\ \lambda v &= -ikp + \nu (\frac{d^2 v}{dx\_1^2} - k^2 v) \ , \\ \frac{du}{dx\_1} + ikv &= 0 \ , \end{aligned} \tag{5.47}$$

for −1 ≤ *x*1 ≤ +1. The boundary conditions are

$$\mathfrak{v}(\pm 1, x\_2, t) = \mathbf{0}\,,\tag{5.48}$$

for no-slip walls. The elimination of v and *p* in Eqs. (5.47) yields

$$
\lambda \left( D^2 - k^2 \right) u = \nu \left( D^2 - k^2 \right)^2 u \,, \tag{5.49}
$$

with *u*(±1) = *Du*(±1) = 0, where *D* = *d*/*dx*1. The solutions of Eq. (5.49) are either symmetric in *x*1

$$
\mu(\mathbf{x}\_1) = \cos \mu \cosh k \mathbf{x}\_1 - \cosh k \cos \mu \mathbf{x}\_1 \,, \tag{5.50}
$$

or antisymmetric

$$
\mu(\mathbf{x}\_1) = \sin \mu \sinh k \mathbf{x}\_1 - \sinh k \sin \mu \mathbf{x}\_1 \,\, . \tag{5.51}
$$

The eigenvalues are

$$
\lambda = -\nu(\mu^2 + k^2)\,, \tag{5.52}
$$

satisfying the relations

$$k \tanh k = -\mu \tan \mu \ , \tag{5.53}$$

for (5.50) and

$$k\coth k = \mu \cot \mu \text{ , }\tag{5.54}$$

for (5.51). For *k* = 1 and *k* = 10, the first symmetric eigenmode decays with λ/ν = −9.3137 and λ/ν = −103.0394, respectively, while, for *k* = 1 and *k* = 10, the antisymmetric eigenmode decays with λ/ν = −20.5706 and λ/ν = −112.0836. All the eigenvalues λ are real and negative, indicating strong damping by the viscous forces.

This problem is an excellent benchmark to test the accuracy of numerical methods treating the velocity-pressure formulation of the Navier–Stokes equations.

#### **5.4 Parallel Flow Around a Sphere**

A sphere of radius *R* is in a viscous steady state flow for which the velocity at infinity upstream is *U*. We assume a creeping flow such that we can have a solution of the Stokes equation (2.54). We place the Cartesian coordinate system such that the axis *x*3 is oriented in the direction of the flow incident on the sphere (Fig. 5.5). The boundary conditions expressed in spherical coordinates (see Fig. B.1) are

$$
\mathfrak{v} = \mathbf{0} \quad \text{at} \quad r = R,\tag{5.55}
$$

$$
\mathfrak{v} = U \mathfrak{e}\_3 \quad \text{at} \quad r = \infty \,. \tag{5.56}
$$

The problem thus posed is symmetric about the axis *Ox*3, and with respect to the longitude of the sphere. Consequently, ∂(•)/∂ϕ ≡ 0. Also then, vϕ = 0. Thus the mass conservation equation (B.20) reduces to

$$\frac{1}{r^2}\frac{\partial}{\partial r}(r^2v\_r) + \frac{1}{r\sin\theta}\frac{\partial}{\partial \theta}(v\_\theta \sin\theta) = 0 \,\,. \tag{5.57}$$

We deduce that a stream function ψ exists such that

$$v\_r = \frac{1}{r^2 \sin \theta} \frac{\partial \psi}{\partial \theta}, \quad v\_\theta = -\frac{1}{r \sin \theta} \frac{\partial \psi}{\partial r} \,. \tag{5.58}$$

Given the two-dimensional character of the flow, the vorticity will have a single component in the direction of the vector *e*ϕ that we denote ω. We write (cf. (B.5))

$$\rho(r,\theta) = -\frac{1}{r} \left[ \frac{1}{\sin\theta} \frac{\partial^2 \psi}{\partial r^2} + \frac{1}{r^2} \frac{\partial}{\partial \theta} \left( \frac{1}{\sin\theta} \frac{\partial \psi}{\partial \theta} \right) \right]. \tag{5.59}$$

In the case of the Stokes equation, vorticity is a harmonic function. We have (recall that the Laplacian of a vector is not equal to the Laplacian of its components, cf. (B.7))

$$\begin{split} \Delta\omega - \frac{\omega}{r^2 \sin^2 \theta} &= \frac{1}{r^2} \left( \frac{\partial}{\partial r} \left( r^2 \frac{\partial \omega}{\partial r} \right) + \frac{1}{\sin \theta} \frac{\partial}{\partial \theta} \left( \sin \theta \frac{\partial \omega}{\partial \theta} \right) - \frac{\omega}{\sin^2 \theta} \right) \\ &= \frac{1}{r} \frac{\partial^2}{\partial r^2} (r \omega) + \frac{1}{r^2} \frac{\partial}{\partial \theta} \left( \frac{1}{\sin \theta} \frac{\partial}{\partial \theta} \left( \omega \sin \theta \right) \right) \\ &= 0 \end{split} \tag{5.60}$$

The combination of relations (5.59)–(5.61) gives the following biharmonic equation:

$$
\left(\frac{\partial^2}{\partial r^2} + \frac{\sin\theta}{r^2} \frac{\partial}{\partial \theta} \left(\frac{1}{\sin\theta} \frac{\partial}{\partial \theta}\right)\right)^2 \psi = 0 \,. \tag{5.62}
$$

The boundary conditions (5.55) and (5.56), expressed in terms of the stream function, become

$$\begin{aligned} \frac{\partial \Psi}{\partial \theta} &= \frac{\partial \Psi}{\partial r} = 0, \quad \text{at} \quad r = R, \\\ v\_r &= U \cos \theta, \quad \frac{\partial \Psi}{\partial \theta} = Ur^2 \sin \theta \cos \theta, \quad \text{at} \quad r = \infty \end{aligned}$$

$$\begin{aligned} v\_\theta &= -U \sin \theta, \quad \frac{\partial \Psi}{\partial r} = Ur \sin^2 \theta. \end{aligned} \tag{5.63}$$

The condition at infinity can be easily integrated. It follows that

$$
\psi\_{\infty} = \frac{1}{2} U r^2 \sin^2 \theta \ . \tag{5.64}
$$

The form of this expression for ψ suggests that the stream function can be written in the general form

$$
\psi = \sin^2 \theta \, f(r) \,. \tag{5.65}
$$

Introducing (5.65) into (5.62), we find

$$\frac{d^4f}{dr^4} - \frac{4}{r^2}\frac{d^2f}{dr^2} + \frac{8}{r^3}\frac{df}{dr} - \frac{8}{r^4}f = 0\ . \tag{5.66}$$

Seeking a solution as a power series in *r n*, we obtain the characteristic polynomial

$$(n-2)(n-1)\left(n^2 - 3n - 4\right) = 0\,,$$

whose roots are *n* = −1, 1, 2, and 4. The function *f* (*r*) is thus

$$f = \frac{C\_{-1}}{r} + C\_1 r + C\_2 r^2 + C\_4 r^4 \,. \tag{5.67}$$

The imposition of the boundary condition at infinity, (5.64), requires *C*4 = 0 and *C*2 = 1 2*U*, while at the sphere, with v*r* = vθ = 0, we can determine *C*−1 = (1/4)*U R*3,*C*1 = −(3/4)*U R*. The stream function is then

$$\psi = \frac{U\,R^2}{2}\sin^2\theta \left(\frac{R}{2r} - \frac{3}{2}\frac{r}{R} + (\frac{r}{R})^2\right) \,. \tag{5.68}$$

We can easily deduce the velocities from (5.58)

$$v\_r = U\cos\theta \left[\frac{1}{2}(\frac{R}{r})^3 - \frac{3}{2}\frac{R}{r} + 1\right] \tag{5.69}$$

$$v\_{\theta} = U \sin \theta \left[ \frac{1}{4} (\frac{R}{r})^3 + \frac{3}{4} \frac{R}{r} - 1 \right]. \tag{5.70}$$

The Cartesian velocity components are linked to the spherical components by the transformation

$$
\begin{pmatrix} v\_1 \\ v\_2 \\ v\_3 \end{pmatrix} = \begin{pmatrix} \sin\theta\cos\varphi\cos\theta\cos\varphi - \sin\varphi \\\\ \sin\theta\sin\varphi\cos\theta\sin\varphi & \cos\varphi \\\\ \cos\theta & -\sin\theta & 0 \end{pmatrix} \begin{pmatrix} v\_r \\ v\_\theta \\ v\_\varphi \end{pmatrix} . \tag{5.71}
$$

Using (5.71), Eqs. (5.69) and (5.70) yield

$$\begin{aligned} v\_1 &= -\frac{3}{4} \frac{U R x\_1 x\_3}{r^3} \left( 1 - \frac{R^2}{r^2} \right), \\ v\_2 &= -\frac{3}{4} \frac{U R x\_2 x\_3}{r^3} \left( 1 - \frac{R^2}{r^2} \right), \\ v\_3 &= U \left[ 1 - \frac{3}{4} \frac{R x\_3^2}{r^3} \left( 1 - \frac{R^2}{r^2} \right) - \frac{1}{4} \frac{R}{r} \left( 3 + \frac{R^2}{r^2} \right) \right], \end{aligned} \tag{5.72}$$

as we have *x*1 = *r* sin θ cos ϕ *x*2 = *r* sin θ sin ϕ *x*3 = *r* cos θ .

The vorticity field is written (cf. (5.59))

$$
\omega = -\frac{3}{2} U R \left( \frac{\sin \theta}{r^2} \right) . \tag{5.73}
$$

The calculation of the pressure field can easily be accomplished by taking into account the vector identity (4.20) which leads to the Stokes equation

#### 5.4 Parallel Flow Around a Sphere 127

$$\nabla p = -\mu \,\,\mathbf{curl}\,\boldsymbol{\omega}\,. \tag{5.74}$$

With (B.5), we generate the system of equations

$$\frac{\partial p}{\partial r} = -\frac{\mu}{r \sin \theta} \frac{\partial}{\partial \theta} \left( \omega \, \sin \theta \right) = 3 \mu \, U \, R \frac{\cos \theta}{r^3} \,, \tag{5.75}$$

$$\frac{1}{r}\frac{\partial p}{\partial \theta} = \frac{\mu}{r}\frac{\partial (r\omega)}{\partial r} = \frac{3\mu}{2}U\,R\frac{\sin\theta}{r^3} \,. \tag{5.76}$$

Integration of (5.75) yields

$$p = -\frac{3\mu}{2} \, U \, R \frac{\cos \theta}{r^2} + q(\theta) \, .$$

Inserting this result in (5.76), we have

$$\frac{3\mu}{2}\,U\,R\frac{\sin\theta}{r^3} + \frac{q'(\theta)}{r} = \frac{3\mu}{2}\,U\,R\frac{\sin\theta}{r^3}\,.$$

The pressure field is then finally given by

$$p = -\frac{3\mu}{2} \, U \, R \frac{\cos \theta}{r^2} + p\_0 \,, \tag{5.77}$$

with *p*0 a constant reference pressure.

The uniform velocity flow around a sphere will generate pressure and shear forces. To calculate the pressure force in direction *Ox*3, we integrate the elementary forces over the surface of the sphere

$$dF\_{3,p} = - (\frac{3\mu}{2} \, U \, \frac{\cos \theta}{R} + p\_0) \cos \theta \, (2\pi R^2 \sin \theta) \, d\theta \, . \tag{5.78}$$

The factor 2π comes from the symmetry of the problem which allows us to take into account the longitudinal part of the integral. Integrating from θ = 0 to θ = π, we obtain

$$F\_{3,p} = -2\pi \,\mu \, U \, R \, . \tag{5.79}$$

Friction drag is obtained by integration over the sphere of the shear stress that acts on it, that is, σ*r*θ (cf. (B.19)) which is −3μ*U* sin θ/(2*R*) for *r* = *R*. This leads to

$$F\_{3,\sigma} = -\int\_{\theta=0}^{\theta=\pi} \left(\sigma\_{r\theta} \mid\_{r=R} \sin \theta\right) \left(2\pi R^2 \sin \theta\right) d\theta = -4\pi \,\mu U \, R \,. \tag{5.80}$$

Total drag, *F*3 = *F*3,*p* + *F*3,σ , known as *Stokes drag*, is the sum of the pressure force and the friction force

$$F\_3 = -6\pi \,\mu \, U \, R \, . \tag{5.81}$$

If we define the drag coefficient by

$$C\_D = \frac{|F\_3|}{\frac{1}{2}\rho U^2 \,\pi R^2},\tag{5.82}$$

we obtain

$$C\_D = \frac{24}{Re},\tag{5.83}$$

where *Re* = 2*U R*/ν. Note that the pressure drag represents one third of the total drag. Relation (5.83) is verified by experiments when *Re* < 1 which is valid in the neighborhood of the sphere. When we move away, the importance of the inertial terms grows and the Stokes solution diverges from the exact solution. Note that the solution that we have obtained is not applicable to the case of a set of spherical particles, as the presence of a spherical obstacle in the flow has impact relatively far away since the velocity profiles decrease as 1/*r*.

The solution for uniform flow around a fixed sphere can be transposed to the case of translation at uniform velocity *U* of a sphere of radius *R* in a fluid at rest at infinity. In this case, the coordinate system is still attached to the sphere and thus in translation at uniform velocity. This modifies the sign of *U* to become −*U* for the pressure and the vorticity. As for the velocity in the fluid, this is relative to the coordinate system, which leads to the following modifications: for the velocity and the stream function, *U* becomes −*U* and the uniform velocity field must also be subtracted from the corresponding relations.

#### *5.4.1 Oseen's Improvement*

The Stokes solution was improved by Oseen [70] who proposed the solution of the Navier–Stokes equations (1.74) as a sum of uniform velocity field and a perturbation such that

$$
\mathfrak{v} = U \mathfrak{e}\_3 + \mathfrak{v}'.\tag{5.84}
$$

In the case of the flow around a fixed sphere, the velocity *v* then takes into account the perturbation caused by the sphere in a flow uniform at infinity. With (5.84), the stationary inertial term takes the form

$$
\rho \frac{D\mathbf{v}}{Dt} = \rho \left( v\_j' \frac{\partial v\_i'}{\partial x\_j} + U \frac{\partial v\_i'}{\partial x\_3} \right) \,. \tag{5.85}
$$

Oseen'sassumptionamountsto neglectingthefirsttermwith respecttothe secondterm on the right-hand side of Eq. (5.85). We obtain a linearized Navier–Stokes equation

$$
\rho \, U \frac{\partial \mathbf{v}'}{\partial \mathbf{x}\_3} = -\nabla \, p + \mu \, \Delta \mathbf{v}' + \rho \mathbf{b} \, . \tag{5.86}
$$

The drag coefficient obtained with Oseen's solution is

$$C\_D = \frac{24}{Re} (1 + \frac{3}{16} Re) \,\text{.}\tag{5.87}$$

Experimental results show that (5.87) is approximately valid for *Re* < 5. Using matched asymptotic expansions [49], the corrected coefficient becomes

$$C\_D = \frac{24}{Re} (1 + \frac{3}{16} Re - \frac{19}{1280} Re^2 + O(Re^3))\,. \tag{5.88}$$

Another approximation for the drag coefficient is given as

$$C\_D = \frac{24}{Re} [1 + \frac{3}{16} Re + \frac{9}{160} Re^2 \ln Re + O(Re^2)]\,. \tag{5.89}$$

For a quantitative outline of the analysis leading to (5.89), we refer the reader to Proudman and Pearson [75].

A more accurate relation for the drag is proposed by Ockendon and Ockendon [64]:

$$C\_D = \frac{24}{Re} [1 + \frac{3}{16} Re + \frac{9}{160} Re^2 \ln Re + \frac{1}{160} (9\chi + 15 \ln 2 - \frac{323}{40}) Re^2 + \dotsb],$$

where γ = 0.5772156649 is Euler's constant.

#### **5.5 Parallel Flow Around a Cylinder**

After the detailed study of the flow around a sphere, we would like to investigate the case of a uniform steady state parallel flow impinging a fixed circular cylinder. We will unfold the same steps as in the previous section in order to compare both theoretical developments and draw conclusions about the impact of the geometrical configurations on the analytical results.

A cylinder of radius *R* is placed in a uniform parallel flow for which the velocity at infinite upstream is *U*. We assume again a creeping flow described by the Stokes equation (2.54). We place the Cartesian coordinate system such that the axis *x*3 is oriented in the direction of the flow incident on the cylinder (Fig. 5.6). The boundary conditions expressed in cylindrical or more precisely in polar coordinates are

$$
\mathbf{v} = \mathbf{0} \quad \text{at} \quad r = R,\tag{5.90}
$$

$$
\mathfrak{v} = U \mathfrak{e}\_3 \quad \text{at} \quad r = \infty \,. \tag{5.91}
$$

The flow is two-dimensional and therefore v*z* = 0. The velocity field is such that *v* = (v*r*(*r*, θ ), vθ (*r*, θ ), 0). The mass conservation equation (A.2) reduces to

**Fig. 5.6** Flow around a circular cylinder

$$\frac{1}{r}\frac{\partial}{\partial r}\left(rv\_r\right) + \frac{1}{r}\frac{\partial v\_\theta}{\partial \theta} = 0\ . \tag{5.92}$$

We deduce that a stream function ψ exists such that

$$v\_r = \frac{1}{r} \frac{\partial \Psi}{\partial \theta}, \quad v\_\theta = -\frac{\partial \Psi}{\partial r} \,. \tag{5.93}$$

Given the two-dimensional character of the flow, the vorticity will have a single component in the direction orthogonal to the plane that we denote ω. We can then write (see Eq. (A.6))

$$\omega(r,\theta) = -\frac{1}{r} \left( \frac{\partial}{\partial r} \left( r v\_{\theta} \right) - \frac{\partial v\_{r}}{\partial \theta} \right) \tag{5.94}$$

$$=-\left(\frac{1}{r}\frac{\partial\psi}{\partial r}+\frac{\partial^2\psi}{\partial r^2}+\frac{1}{r^2}\frac{\partial^2\psi}{\partial\theta^2}\right).\tag{5.95}$$

As vorticity is a harmonic function, its Laplacian is given by

$$
\Delta\omega = \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial\omega}{\partial r}\right) + \frac{1}{r^2}\frac{\partial^2\omega}{\partial\theta^2} = \frac{1}{r}\frac{\partial\omega}{\partial r} + \frac{\partial^2\omega}{\partial r^2} + \frac{1}{r^2}\frac{\partial^2\omega}{\partial\theta^2}.\tag{5.96}
$$

The combination of relations (5.95)–(5.96) gives the following biharmonic equation:

$$\left(\frac{\partial^2}{\partial r^2} + \frac{1}{r}\frac{\partial}{\partial r} + \frac{1}{r^2}\frac{\partial^2}{\partial \theta^2}\right)^2 \psi = 0 \,. \tag{5.97}$$

The boundary conditions (5.90) and (5.91), expressed in terms of the stream function, become

#### 5.5 Parallel Flow Around a Cylinder 131

$$\frac{\partial \psi}{\partial \theta} = \frac{\partial \psi}{\partial r} = 0, \quad \text{at} \quad r = R,\tag{5.98}$$

$$v\_r = U\cos\theta, \quad \frac{\partial\psi}{\partial\theta} = Ur\cos\theta, \quad \text{at} \quad r = \infty$$

$$v\_\theta = -U\sin\theta, \quad \frac{\partial\psi}{\partial r} = U\sin\theta \,. \tag{5.99}$$

The condition at infinity suggests that the stream function can be written in the general form

$$
\psi = f(r) \sin \theta \,\,. \tag{5.100}
$$

Introducing (5.100) into (5.97), we find

$$\frac{d^4f}{dr^4} + \frac{2}{r}\frac{d^3f}{dr^3} - \frac{3}{r^2}\frac{d^2f}{dr^2} + \frac{3}{r^3}\frac{df}{dr} - \frac{3}{r^4}f = 0 \,\,. \tag{5.101}$$