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)^[1] and Georg Hamel(1917),^[2] but it has subsequently been studied by many major scientists such as von Kármán and Levi-Civita,^[3] Walter Tollmien,^[4] F. Noether,^[5] W.R. Dean,^[6] Rosenhead,^[7] Landau,^[8] G.K. Batchelor^[9] etc. A complete set of solutions was described by Edward Fraenkel in 1962.^[10]

Flow description

Consider two stationary plane walls with a constant volume flow rate ${\displaystyle Q}$ is injected/sucked at the point of intersection of plane walls and let the angle subtended by two walls be ${\displaystyle 2\alpha }$. Take the cylindrical coordinate ${\displaystyle (r,\theta ,z)}$ system with ${\displaystyle r=0}$ representing point of intersection and ${\displaystyle \theta =0}$ the centerline and ${\displaystyle (u,v,w)}$ are the corresponding velocity components. The resulting flow is two-dimensional if the plates are infinitely long in the axial ${\displaystyle z}$ 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., ${\displaystyle u=u(r,\theta ),v=0,w=0}$.

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

${\displaystyle {\begin{aligned}{\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{aligned}}}$

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.

${\displaystyle u(\pm \alpha )=0,\quad Q=\int _{-\alpha }^{\alpha }ur\,d\theta }$

Formulation

The first equation tells that ${\displaystyle ru}$ is just function of ${\displaystyle \theta }$, the function is defined as

${\displaystyle F(\theta )={\frac {ru}{\nu }}.}$

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

${\displaystyle {\frac {1}{\rho }}{\frac {\partial p}{\partial \theta }}={\frac {2\nu ^{2}}{r^{2}}}{\frac {dF}{d\theta }}}$

Now letting

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

the ${\displaystyle r}$ and ${\displaystyle \theta }$ momentum equations reduce to

${\displaystyle P=-{\frac {1}{2}}(F^{2}+F'')}$

${\displaystyle P'=2F',\quad \Rightarrow \quad P=2F+C}$

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

${\displaystyle F''+F^{2}+4F+2C=0}$

Multiplying by ${\displaystyle F'}$ and integrating once,

${\displaystyle {\frac {1}{2}}F'^{2}+{\frac {1}{3}}F^{3}+2F^{2}+2CF=D,}$

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

where ${\displaystyle C,D}$ are constants to be determined from the boundary conditions. The above equation can be re-written conveniently with three other constants ${\displaystyle a,b,c}$ 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 ${\displaystyle a+b+c=-6}$.

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

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

The boundary conditions reduce to

${\displaystyle F(\pm \alpha )=0,\quad {\frac {Q}{\nu }}=\int _{-\alpha }^{\alpha }F\,d\theta }$

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

Dynamical interpretation^[13]

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

${\displaystyle V(F)=-{\frac {1}{3}}(a-F)(F-b)(F-c)}$

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

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

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

Pure outflow

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

${\displaystyle \theta ={\sqrt {\frac {3}{2}}}\int _{F}^{a}{\frac {dF}{\sqrt {(a-F)(F-b)(F-c))}}}}$

and the boundary conditions becomes

${\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))}}}.}$

The equations can be simplified by standard transformations given for example in Jeffreys.^[14]

• First case ${\displaystyle b,c}$ are complex conjugates and ${\displaystyle a>0}$ leads to

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

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

where ${\displaystyle \operatorname {sn} ,\operatorname {cn} }$ are Jacobi elliptic functions.

• Second case ${\displaystyle c<b<0<a}$ leads to

${\displaystyle F(\theta )=a-6k^{2}m^{2}\operatorname {sn} ^{2}(m\theta ,k)}$

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

Limiting form

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

${\displaystyle {\begin{aligned}\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{aligned}}}$

where

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

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

The corresponding critical Reynolds number or volume flux is given by

${\displaystyle {\begin{aligned}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}tdt,\\&={\frac {12}{\sqrt {1-2k^{2}}}}[E(k^{2})-(1-k^{2})K(k^{2})]\end{aligned}}}$

where ${\displaystyle E(k^{2})}$ is the complete elliptic integral of the second kind. For large values of ${\displaystyle a,\left(\ k^{2}\sim {\frac {1}{2}}-{\frac {3}{2a}}\right)}$, the critical Reynolds number or volume flux becomes ${\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}}}$.

Pure inflow

For pure inflow, the implicit solution is given by

${\displaystyle \theta ={\sqrt {\frac {3}{2}}}\int _{b}^{F}{\frac {dF}{\sqrt {(a-F)(F-b)(F-c))}}}}$

and the boundary conditions becomes

${\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))}}}.}$

Pure inflow is possible only when all constants are real ${\displaystyle c<b<0<a}$ and the solution is given by

${\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)}$

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

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

Limiting form

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

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

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

References