Physics:Jeffery–Hamel flow

In fluid dynamics Jeffery–Hamel flow is a flow created by a converging or diverging channel with a source or sink of fluid volume at the point of intersection of the two plane walls. It is named after George Barker Jeffery(1915) and Georg Hamel(1917), but it has subsequently been studied by many major scientists such as von Kármán and Levi-Civita, Walter Tollmien, F. Noether, W.R. Dean, Rosenhead, Landau, G.K. Batchelor etc. A complete set of solutions was described by Edward Fraenkel in 1962.

Flow description

Consider two stationary plane walls with a constant volume flow rate [math]\displaystyle{ Q }[/math] is injected/sucked at the point of intersection of plane walls and let the angle subtended by two walls be [math]\displaystyle{ 2\alpha }[/math]. Take the cylindrical coordinate [math]\displaystyle{ (r,\theta,z) }[/math] system with [math]\displaystyle{ r=0 }[/math] representing point of intersection and [math]\displaystyle{ \theta=0 }[/math] the centerline and [math]\displaystyle{ (u,v,w) }[/math] are the corresponding velocity components. The resulting flow is two-dimensional if the plates are infinitely long in the axial [math]\displaystyle{ z }[/math] direction, or the plates are longer but finite, if one were neglect edge effects and for the same reason the flow can be assumed to be entirely radial i.e., [math]\displaystyle{ u=u(r,\theta),v=0,w=0 }[/math].

Then the continuity equation and the incompressible Navier–Stokes equations reduce to

[math]\displaystyle{ \begin{align} \frac{\partial (ru)}{\partial r} & =0, \\[6pt] u\frac{\partial u}{\partial r} & = - \frac{1}{\rho}\frac{\partial p}{\partial r} + \nu \left[\frac{1}{r} \frac{\partial}{\partial r} \left(r\frac{\partial u}{\partial r}\right) + \frac{1}{r^2}\frac{\partial^2 u}{\partial \theta^2}- \frac{u}{r^2}\right] \\[6pt] 0 & = - \frac{1}{\rho r} \frac{\partial p}{\partial \theta} + \frac{2 \nu}{r^2} \frac{\partial u}{\partial \theta} \end{align} }[/math]

The boundary conditions are no-slip condition at both walls and the third condition is derived from the fact that the volume flux injected/sucked at the point of intersection is constant across a surface at any radius.

[math]\displaystyle{ u(\pm \alpha) = 0, \quad Q= \int_{-\alpha}^\alpha u r \, d\theta }[/math]

Formulation

The first equation tells that [math]\displaystyle{ ru }[/math] is just function of [math]\displaystyle{ \theta }[/math], the function is defined as

[math]\displaystyle{ F(\theta) = \frac{r u}{ \nu}. }[/math]

Different authors defines the function differently, for example, Landau defines the function with a factor [math]\displaystyle{ 6 }[/math]. But following Whitham, Rosenhead the [math]\displaystyle{ \theta }[/math] momentum equation becomes

[math]\displaystyle{ \frac{1}{\rho}\frac{\partial p}{\partial\theta} = \frac{2 \nu^2}{r^2} \frac{dF}{d\theta} }[/math]

Now letting

[math]\displaystyle{ \frac{p-p_\infty}{\rho} = \frac{\nu^2}{r^2} P(\theta), }[/math]

the [math]\displaystyle{ r }[/math] and [math]\displaystyle{ \theta }[/math] momentum equations reduce to

[math]\displaystyle{ P = -\frac{1}{2} (F^2+F'') }[/math]

[math]\displaystyle{ P'= 2F', \quad \Rightarrow \quad P = 2F + C }[/math]

and substituting this into the previous equation(to eliminate pressure) results in

[math]\displaystyle{ F'' + F^2 + 4F + 2C =0 }[/math]

Multiplying by [math]\displaystyle{ F' }[/math] and integrating once,

[math]\displaystyle{ \frac{1}{2} F'^2+ \frac{1}{3} F^3 +2F^2 +2CF = D, }[/math]

[math]\displaystyle{ \frac{1}{2} F'^2+ \frac{1}{3} (F^3 +6F^2 +6CF-3D) = 0 }[/math]

where [math]\displaystyle{ C,D }[/math] are constants to be determined from the boundary conditions. The above equation can be re-written conveniently with three other constants [math]\displaystyle{ a,b,c }[/math] as roots of a cubic polynomial, with only two constants being arbitrary, the third constant is always obtained from other two because sum of the roots is [math]\displaystyle{ a+b+c=-6 }[/math].

[math]\displaystyle{ \frac{1}{2}F'^2+\frac{1}{3}(F-a)(F-b)(F-c)=0, }[/math]

[math]\displaystyle{ \frac{1}{2}F'^2-\frac{1}{3}(a-F)(F-b)(F-c)=0. }[/math]

The boundary conditions reduce to

[math]\displaystyle{ F(\pm \alpha) = 0, \quad \frac{Q}{\nu}= \int_{-\alpha}^\alpha F \, d\theta }[/math]

where [math]\displaystyle{ Re=Q/\nu }[/math] is the corresponding Reynolds number. The solution can be expressed in terms of elliptic functions. For convergent flow [math]\displaystyle{ Q\lt 0 }[/math], the solution exists for all [math]\displaystyle{ Re }[/math], but for the divergent flow [math]\displaystyle{ Q\gt 0 }[/math], the solution exists only for a particular range of [math]\displaystyle{ Re }[/math].

Dynamical interpretation

The equation takes the same form as an undamped nonlinear oscillator(with cubic potential) one can pretend that [math]\displaystyle{ \theta }[/math] is time, [math]\displaystyle{ F }[/math] is displacement and [math]\displaystyle{ F' }[/math] is velocity of a particle with unit mass, then the equation represents the energy equation([math]\displaystyle{ K.E. + P.E.=0 }[/math], where [math]\displaystyle{ K.E. = \frac{1}{2} F'^2 }[/math] and [math]\displaystyle{ P.E. = V(F) }[/math] ) with zero total energy, then it is easy to see that the potential energy is

[math]\displaystyle{ V(F)=-\frac{1}{3}(a-F)(F-b)(F-c) }[/math]

where [math]\displaystyle{ V\leq 0 }[/math] in motion. Since the particle starts at [math]\displaystyle{ F=0 }[/math] for [math]\displaystyle{ \theta=-\alpha }[/math] and ends at [math]\displaystyle{ F=0 }[/math] for [math]\displaystyle{ \theta=\alpha }[/math], there are two cases to be considered.

• First case [math]\displaystyle{ b,c }[/math] are complex conjugates and [math]\displaystyle{ a\gt 0 }[/math]. The particle starts at [math]\displaystyle{ F=0 }[/math] with finite positive velocity and attains [math]\displaystyle{ F=a }[/math] where its velocity is [math]\displaystyle{ F'=0 }[/math] and acceleration is [math]\displaystyle{ F''=-dV/dF\lt 0 }[/math] and returns to [math]\displaystyle{ F=0 }[/math] at final time. The particle motion [math]\displaystyle{ 0\lt F\lt a }[/math] represents pure outflow motion because [math]\displaystyle{ F\gt 0 }[/math] and also it is symmetric about [math]\displaystyle{ \theta=0 }[/math].
• Second case [math]\displaystyle{ c\lt b\lt 0\lt a }[/math], all constants are real. The motion from [math]\displaystyle{ F=0 }[/math] to [math]\displaystyle{ F=a }[/math] to [math]\displaystyle{ F=0 }[/math] represents a pure symmetric outflow as in the previous case. And the motion [math]\displaystyle{ F=0 }[/math] to [math]\displaystyle{ F=b }[/math] to [math]\displaystyle{ F=0 }[/math] with [math]\displaystyle{ F\lt 0 }[/math] for all time([math]\displaystyle{ -\alpha\leq\theta\leq\alpha }[/math]) represents a pure symmetric inflow. But also, the particle may oscillate between [math]\displaystyle{ b\leq F\leq a }[/math], representing both inflow and outflow regions and the flow is no longer need to symmetric about [math]\displaystyle{ \theta=0 }[/math].

The rich structure of this dynamical interpretation can be found in Rosenhead(1940).

Pure outflow

For pure outflow, since [math]\displaystyle{ F=a }[/math] at [math]\displaystyle{ \theta=0 }[/math], integration of governing equation gives

[math]\displaystyle{ \theta = \sqrt{\frac{3}{2}} \int_F^a \frac{dF}{\sqrt{(a-F)(F-b)(F-c))}} }[/math]

and the boundary conditions becomes

[math]\displaystyle{ \alpha = \sqrt{\frac{3}{2}} \int_0^a \frac{dF}{\sqrt{(a-F)(F-b)(F-c))}}, \quad Re = 2\sqrt{\frac{3}{2}} \int_0^\alpha \frac{FdF}{\sqrt{(a-F)(F-b)(F-c))}}. }[/math]

The equations can be simplified by standard transformations given for example in Jeffreys.

• First case [math]\displaystyle{ b,c }[/math] are complex conjugates and [math]\displaystyle{ a\gt 0 }[/math] leads to

[math]\displaystyle{ F(\theta)=a-\frac{3M^2}{2}\frac{1-\operatorname{cn}(M\theta,\kappa)}{1+\operatorname{cn}(M\theta,\kappa)} }[/math]

[math]\displaystyle{ M^2 = \frac{2}{3} \sqrt{(a-b)(a-c)}, \quad \kappa^2=\frac{1}{2}+\frac{a+2}{2M^2} }[/math]

where [math]\displaystyle{ \operatorname{sn}, \operatorname{cn} }[/math] are Jacobi elliptic functions.

• Second case [math]\displaystyle{ c\lt b\lt 0\lt a }[/math] leads to

[math]\displaystyle{ F(\theta)=a-6k^2 m^2\operatorname{sn}^2(m\theta,k) }[/math]

[math]\displaystyle{ m^2 = \frac{1}{6} (a-c), \quad k^2=\frac{a-b}{a-c}. }[/math]

Limiting form

The limiting condition is obtained by noting that pure outflow is impossible when [math]\displaystyle{ F'(\pm\alpha)=0 }[/math], which implies [math]\displaystyle{ b=0 }[/math] from the governing equation. Thus beyond this critical conditions, no solution exists. The critical angle [math]\displaystyle{ \alpha_c }[/math] is given by

[math]\displaystyle{ \begin{align} \alpha_c &= \sqrt{\frac{3}{2}} \int_0^a \frac{dF}{\sqrt{F(a-F)(F+a+6))}},\\ &= \sqrt{\frac{3}{2a}} \int_0^1 \frac{dt}{\sqrt{t(1-t)\{1+(1+6/a)t\}}},\\ &= \frac{K(k^2)}{m^2} \end{align} }[/math]

where

[math]\displaystyle{ m^2 = \frac{3+a}{3}, \quad k^2 = \frac{1}{2}\left(\frac{a}{3+a}\right) }[/math]

where [math]\displaystyle{ K(k^2) }[/math] is the complete elliptic integral of the first kind. For large values of [math]\displaystyle{ a }[/math], the critical angle becomes [math]\displaystyle{ \alpha_c =\sqrt{\frac{3}{a}}K\left(\frac{1}{2}\right)=\frac{3.211}{\sqrt{a}} }[/math].

The corresponding critical Reynolds number or volume flux is given by

[math]\displaystyle{ \begin{align} Re_c = \frac{Q_c}{\nu} &= 2 \int_0^{\alpha_c} (a-6k^2 m^2\operatorname{sn}^2 m\theta) \, d\theta,\\ &= \frac{12k^2}{\sqrt{1-2k^2}} \int_0^K \operatorname{cn}^2 t dt,\\ &= \frac{12}{\sqrt{1-2k^2}}[ E(k^2) -(1-k^2)K(k^2)] \end{align} }[/math]

where [math]\displaystyle{ E(k^2) }[/math] is the complete elliptic integral of the second kind. For large values of [math]\displaystyle{ a, \left(\ k^2\sim \frac{1}{2}-\frac{3}{2a}\right) }[/math], the critical Reynolds number or volume flux becomes [math]\displaystyle{ Re_c=\frac{Q_c}{\nu} = 12 \sqrt{\frac{a}{3}} \left[E\left(\frac{1}{2}\right)-\frac{1}{2}K\left(\frac{1}{2}\right) \right]=2.934 \sqrt{a} }[/math].

Pure inflow

For pure inflow, the implicit solution is given by

[math]\displaystyle{ \theta = \sqrt{\frac{3}{2}} \int_b^F \frac{dF}{\sqrt{(a-F)(F-b)(F-c))}} }[/math]

and the boundary conditions becomes

[math]\displaystyle{ \alpha = \sqrt{\frac{3}{2}} \int_b^0 \frac{dF}{\sqrt{(a-F)(F-b)(F-c))}}, \quad Re = 2\sqrt{\frac{3}{2}} \int_\alpha^0 \frac{FdF}{\sqrt{(a-F)(F-b)(F-c))}}. }[/math]

Pure inflow is possible only when all constants are real [math]\displaystyle{ c\lt b\lt 0\lt a }[/math] and the solution is given by

[math]\displaystyle{ F(\theta)=a-6k^2 m^2\operatorname{sn}^2(K-m\theta,k)=b+6k^2 m^2 \operatorname{cn}^2(K-m\theta,k) }[/math]

[math]\displaystyle{ m^2 = \frac{1}{6} (a-c), \quad k^2=\frac{a-b}{a-c} }[/math]

where [math]\displaystyle{ K(k^2) }[/math] is the complete elliptic integral of the first kind.

Limiting form

As Reynolds number increases ([math]\displaystyle{ -b }[/math] becomes larger), the flow tends to become uniform(thus approaching potential flow solution), except for boundary layers near the walls. Since [math]\displaystyle{ m }[/math] is large and [math]\displaystyle{ \alpha }[/math] is given, it is clear from the solution that [math]\displaystyle{ K }[/math] must be large, therefore [math]\displaystyle{ k\sim 1 }[/math]. But when [math]\displaystyle{ k\approx 1 }[/math], [math]\displaystyle{ \operatorname{sn} t\approx \tanh t, \ c\approx b, \ a\approx -2b }[/math], the solution becomes

[math]\displaystyle{ F(\theta) = b\left\{3\tanh^2 \left[\sqrt{-\frac{b}{2}}(\alpha-\theta)+\tanh^{-1}\sqrt{\frac{2}{3}}\right]-2\right\}. }[/math]

It is clear that [math]\displaystyle{ F\approx b }[/math] everywhere except in the boundary layer of thickness [math]\displaystyle{ O\left(\sqrt{-\frac{b}{2}}\right) }[/math]. The volume flux is [math]\displaystyle{ Q/\nu\approx 2\alpha b }[/math] so that [math]\displaystyle{ |Re|=O(|b|) }[/math] and the boundary layers have classical thickness [math]\displaystyle{ O\left(|Re|^{1/2}\right) }[/math].