MIXING LAYER AND TURBULENT JET FLOW
IN A HELE–SHAW CELL

We perform the following scaling in Eqs. (1) and (3)

Then we neglect all terms of the order of $\varepsilon^{2}$. As a result, we obtain the equations of the horizontal-shear flow of a thin fluid layer in a closed channel

The last equation in (5) is a consequence of the first three ones. We supplement system (5) with boundary conditions (2).

Let us assume that the velocity profile $U(t,x,y)$ satisfies the condition

with the parameters $\beta\in(1,2]$ and $\gamma\in(0,1)$. To describe a three-layer flow, we introduce the averaged velocities in the layers and the root-mean-square deviation (shear velocity $q$) in the intermediate layer:

The process of fluid involvement in the vortex interlayer is assumed to be symmetric and is taken into account by kinematic conditions at the internal boundaries

This means that the rate of fluid entrainment into the vortex interlayer is proportional to the velocity of “large eddies” generated in the vicinity of the interface between the layers due to shear instability [4, 9]. According to [18], the proportionality coefficient $\sigma$ is approximately equal to $0.15$.

In averaging Eqs. (5) over the channel width the necessity of approximate calculation of the integrals of the functions $U^{n}$ ($n=2,3$) with respect to the variable $y$ arises. The following estimates are valid for weakly-shear layers:

As terms of the order of $\varepsilon^{2}$ are ignored in deriving model (5), the velocity $U$ is replaced in the course of averaging the equations of motion in the weakly shear layers by the corresponding averaged values of $u$ and $w$ determined by formulas (7). Using the obvious identity $U=v+(U-v)$, we calculate the integrals in the intermediate mixing layer:

Condition (6) implies the estimate $C_{3}=O(\varepsilon^{3\gamma})\ll\eta vq^{2}=O(\varepsilon^{2\gamma})$, $0<\gamma<1$. Therefore, the value of $C_{3}$ is small compared to other terms and can be omitted [26].

The fulfilment of conditions (6) imposed on the velocity profile $U(t,x,y)$ is sufficient to require at the initial moment of time $t=0$. Indeed, by virtue of (5), the variable $U_{y}$ satisfies the equation

and, consequently, $|U_{y}|$ decreases with time.

With allowance for the adopted assumptions, as a result of averaging Eqs. (5) over the channel width and using boundary conditions (2) and (8), we obtain system (4) for $k=1$ (Model II). The first three equations of system (4) express the mass balance in the layers. Since $\xi+\eta+\zeta=Y={\rm const}$ then the flow rate does not depend on the cross section, i.e. $Q=Q(t)$. The fourth and fifth equations are the local momentum conservation laws in the outer layers, the sixth and seventh equations are the total momentum and energy conservation laws. The last term in the energy equation with the factor $\kappa$ does not directly follow from the averaging over the channel width of the last equation in (5). It is added to take into account the energy dissipation. According to [18], the numerical value of the empirical parameter $\kappa$ lies in the interval from 2 to 8.

Governing equations (4) imply the fulfilment of the conditions

They are the consequence of the compatibility of the equations of balance of mass, momentum and energy in weakly-shear layers. A similar approach for determining the velocity at the interface was used in [4, 9]. Indeed, averaging Eqs. (5) over the width of one of the outer layers yields the system

Averaging over another weakly-shear layer gives the same equations up to the replacement of $u$, $\xi$ and $U|_{y=\xi}$ with $w$, $\zeta$ and $U|_{y=\xi+\eta}$. It is easy to show that the compatibility of system (10) for $\sigma q\neq 0$ is ensured by the condition $(U-u)^{2}|_{y=\xi}=0$.

At large Reynolds numbers, vortex motions of various scales arise as a result of the development of shear instability. Let a certain scale be chosen and a procedure for averaging the unknown variables be constructed. Then the velocity field and pressure can be represented as

where $\overline{U}\geq 0$, $\overline{V}$ and $\overline{p}$ are the average velocity components and pressure; $U^{\prime}$, $V^{\prime}$ and $p^{\prime}$ are the corresponding fluctuations relative to the average variables such that $\overline{U^{\prime}}=\overline{V^{\prime}}=\overline{p^{\prime}}=0$. We substitute representation (11) into system (1) and energy equation (3). After averaging over the selected scale, we obtain the Reynolds equations:

Here $\tau=-\overline{U^{\prime}V^{\prime}}$ is the Reynolds shear stress; $P=\overline{p}+\overline{{U^{\prime}}^{2}}$ is the total pressure; $\hat{q}^{2}=\overline{{U^{\prime}}^{2}}+\overline{{V^{\prime}}^{2}}$ is the specific kinetic energy of the fluctuating motion. On the right-hand sides of Eqs. (1), the true velocities are replaced by their average values without fluctuating terms. In deriving Eqs. (12) an assumption is made about the isotropy of the fluctuating motion, i.e. $\overline{{U^{\prime}}^{2}}=\overline{{V^{\prime}}^{2}}$. Moreover, the terms $\overline{U^{\prime}p^{\prime}}$, $\overline{V^{\prime}p^{\prime}}$ and $\overline{{U^{\prime}}^{n}{V^{\prime}}^{3-n}}$ ($n=0,...,3$) having a higher order of smallness compared to $\hat{q}^{3}$ are neglected. The last term in the fourth equation (12) is added to take into account the outflow of energy into small-scale motion. The characteristic frequency $\omega$ is responsible for the energy dissipation rate of turbulent motion. In view of condition (2) on the side walls of the channel $y=0$ and $y=Y$ the velocity component $\overline{V}$ vanishes.

In a developed turbulent flow the following dependence for the Reynolds shear stress was experimentally confirmed in [27]

Eqs. (12), (13) are the closed system for determining the functions $\overline{U}$, $\overline{V}$, $P$ and $\hat{q}$. Taking into account that the fluid moves mainly along the $Ox$ axis, we proceed to the approximation of a long channel ($\varepsilon^{2}\to 0$). We assume here that the fluctuating components $\overline{{U^{\prime}}^{2}}$, $\overline{U^{\prime}V^{\prime}}$ and the variable $\hat{q}^{2}$ are of the same order. Then the following scaling

in Eqs. (12) and neglecting the leading terms in powers of the small parameter $\varepsilon$ yields

where $\tau$ is defined by formula (13). Following [18], we assume $\omega=\sigma\kappa\hat{q}/\eta$. It should be noted that in scaling (14) the value $\varepsilon$ has the same order of smallness as the parameter $\sigma$, therefore, the terms $\sigma/\varepsilon$ are kept in Eqs. (15). Further, in Eqs. (15) we take $\varepsilon=1$.

Let us average Eqs. (15) over the channel width taking into account a three-layer representation of the flow (Fig. 1). We assume that in the outer weakly-shear layers of width $\xi(t,x)$ and $\zeta(t,x)$ there is no kinetic energy of fluctuating motion ($\hat{q}^{2}=0$), and the flow velocity $\overline{U}$ is equal to $u(t,x)$ and $w(t,x)$, respectively. In the intermediate mixing layer of width $\eta(t,x)$ we replace velocity $\overline{U}(t,x,y)$ by its average value $v(t,x)$ and denote $q^{2}(t,x)$ the average value of the function $\hat{q}^{2}(t,x,y)$. We recall that upon obtaining Eqs. (12) (and (15)), a certain averaging scale is fixed. If, as such a scale, we choose the width of the mixing layer $\eta$, then $\overline{U}=v$. In view of the second equation in (15) the total pressure is $P=p(t,x)$. At the interfaces $y=\xi$ and $y=\xi+\eta$, it is assumed that kinematic conditions (8) are satisfied (with velocities $\overline{U}$ and $\overline{V}$ instead of $U$ and $V$). As a result of averaging Eqs. (15) over the channel width, we obtain system (4) for $k=0$ (Model I).

For further analysis of the equations of motion and the construction of stationary solutions it makes sense to derive some consequences. By virtue of system (4) the variables $u-w$, $v$ and $q$ satisfy the equations

Using the sixth equation in system (4) we eliminate $p_{x}$, replace the energy equation with the differential consequence for the variable $q$, and also express $\eta$ and $v$ in terms of other variables. Then system (4) is transformed to the following equations to determine five unknowns $\xi$, $u$, $\zeta$, $w$ and $q$:

Here function $f_{q}$ is defined in the last formula (16), the variables $\eta$, $v$ and $G$ are expressed as follows:

The flow rate $Q(t)$ is considered to be given. The pressure $p$ on the rigid-lid is found from the fourth or fifth equation of system (4). As it was mentioned above the coefficient $k=0$ for Model I and $k=1$ in the case of Model II. Eqs. (17) is used below for an averaged description of the mixing layer.

Let us represent Eqs. (17) in the form

where $\mathbf{U}=(\xi,u,\zeta,w,q)^{\rm T}$ is the vector of unknown variables, $\mathbf{A}$ and $\mathbf{F}$ is the corresponding matrix of size $5\times 5$ and the right-hand side vector. The eigenvalues $\lambda=\lambda_{i}$ of matrix $\mathbf{A}$ determine the propagation velocity of the characteristics. By virtue the last equation in system (17), Model I ($k=0$) has the contact characteristic $dx/dt=v$. Model II ($k=1$) has the same characteristic since equation $s_{t}+vs_{x}=f_{s}$ holds on for the variable $s=q/\eta$ (here the right-hand side term $f_{s}$ does not contain derivatives of the unknown functions). To determine the others characteristics of system (17), an algebraic equation of the fourth degree arises

The case $\xi=\zeta\to Y/2$ corresponds to the initial point of confluence of two flows of the same width. The polynomial $\chi(\lambda)$ is simplified and takes the form

For both models ($k=0$ and $k=1$), equation $\chi(\lambda)=0$ has two complex and two real roots. It is necessary to note that in the framework of Model I, all real characteristics are contact and propagate with the average velocity $v$ in the intermediate mixing layer. Thus, the system under consideration is not hyperbolic in the vicinity of the confluence of flows. It is associated with the development of the Kelvin–Helmholtz instability at the interface.

The equations describing a stationary mixing layer are derived above. To construct the solution, it is necessary to impose the conditions at the point of mixing layer formation, i.e., to find the asymptotic of the stationary solution (18) as $\eta\to 0$. Without loss of generality, we assume that $\eta=0$ at $x=0$, and the values of the functions at this point are marked by the zero subscript. Let at $x=0$ the parameters of weakly shear layers $u_{0}$, $w_{0}$, $\xi_{0}$ and $\zeta_{0}=Y-\xi_{0}$ be specified (all values are positive). We assume that there are finite limits of the functions $v\to v_{0}>0$, $q\to q_{0}>0$ as $x\to 0$ and find these values.

By virtue of the last two equations in (16) and (in the case $k=1$) the second equation in (4) we have the relations

For model I ($k=0$) it follows from the first equation in (19) that $v_{0}=(u_{0}+w_{0})/2$ and from the second equation in (19) the value of $q_{0}^{2}$ is determined. For model II ($k=1$) relations (19) are reduced to the quadratic equation for determining $r=v_{0}/u_{0}$

Here we denote $r_{0}=w_{0}/u_{0}$. This equation has two real roots $r=r_{1,2}$ for all $r_{0}>0$ if $\kappa>\kappa_{*}$, where $\kappa_{*}\approx 7.657$. From the two roots $r_{1}$ and $r_{2}$, we choose the closest to the value of $(1+r_{0})/2$. If $r_{0}$ tends to unity (the velocities of the merging flows $u_{0}$ and $w_{0}$ are close), then the region of existence of the quadratic equation solution with respect to the parameter $\kappa$ expands. The indicated restriction on the choice of the dissipation parameter $\kappa$ should be taken into account in calculations according to Model II.

We obtain a numerical solution of Eqs. (18) for the parameters of the incoming flow taken from [2] (case A): $u_{0}=0.118$ m/s, $w_{0}=0.238$ m/s. The calculations are performed for a straight channel having 16 m long and 3 m wide. Its depth $h$ is equal to 0.052 m. We assume $\xi_{0}=\zeta_{0}=1.5$ m, $c_{f}=0.003$ and $\kappa=4$. For this and all subsequent tests, we choose $\sigma=0.15$. The solution of the stationary mixing layer problem is not difficult, since it is described by a system of ordinary differential equations. The calculations are carried out in the MATLAB environment using the ode45 function that implements the fourth-order Runge–Kutta method.

Let two streams moving with the above-given velocities $u_{0}$ and $w_{0}$ merge in the cross section $x=0$, where an intermediate mixing layer is formed. For the chosen flow parameters at $x=0$, formulas (19) uniquely determine the values $v_{0}\in(u_{0},w_{0})$ and $q_{0}$. For Model I we have $v_{0}\approx 0.178$, for Model II — $v_{0}\approx 0.171$; in both cases $q_{0}\approx 0.035$. The results of calculation using Eqs. (18) are shown in Fig. 2 by solid curves (Model I) and dashed lines (Model II). As can be seen from the graphs, in this example, Models I and II give almost the same result. A slight difference in the calculations for these models is observed only for the velocity $v$ at the initial stage of the mixing layer formation. This is due to the difference in definitions of $v_{0}$ and $q_{0}$ by formulas (19) for $k=0$ and $k=1$. Dashed-dotted lines in Fig. 2 correspond to the similar solution (in the framework of Model II) for the flow in an open channel. The solution is obtained in [4] and verified by experimental data [2]. Comparison of the calculating results of the mixing layer evolution for shear flows in open and closed channels is not entirely correct. Nevertheless, at the initial stage of the mixing layer formation, a good agreement between the solutions of the equations describing the flows under a rigid-lid and with a free surface is observed.

The calculation results using Eqs. (18) for Models I ($k=0$) and II ($k=1$) practically coincide if the ratio of the flow velocities in the input section satisfies the inequality $1/4<w_{0}/u_{0}<4$. As a rule, for problems of jet flows one of the outer fluid layers is at rest. Therefore, the relation $w_{0}/u_{0}\to 0$ (or $w_{0}/u_{0}\to\infty$). In this case, the determination of $v_{0}$ and $q_{0}$ from relations (19) for $k=1$ is possible only for sufficiently large values of the dissipation parameter $\kappa>\kappa_{*}$. In addition, in the framework of Model II, the value of $v_{0}$ differs markedly from the half-sum of the outer layers velocities $(u_{0}+w_{0})/2$. For this reason, in such problems to apply Model I is preferable. To avoid the singularity in Eqs. (18), the fluid velocity in the resting layer is assumed to be positive, but close to zero.

The mixing layer boundaries obtained by the proposed model are in a good agreement with the experimental data and the results of direct numerical simulation performed in [22, 23] for jet Hele–Shaw flows at the Reynolds numbers of the order of $10^{4}$. We consider jet flow in a closed channel of the size $L\times Y\times h$ (corresponds to the directions $x\times y\times z$), where $h=2.4$ mm, $L=267\,h$ and $Y=200\,h$. The channel is filled with fluid at rest. Through a slot of width $9.6\,h$ located in the centre of the cross-section $x=0$ a fluid with velocity $u_{0}=3.6$ m/s is supplied [23]. Below we use this velocity as a scale for the mixing layer reconstruction. Therefore, in dimensionless variables, the jet velocity is chosen equal to one. To visualize the fluid flow at the boundaries of the inlet section, the colouring fluid is injected using needles. As a result of the development of Kelvin–Helmholtz instability, turbulent mixing layers are formed at the interfaces between the jet and the fluid at rest.

Let us choose the following parameters $L=200\,h$, $Y=100\,h$, $\xi_{0}=4.8\,h$, $u_{0}=1$, $w_{0}=0.001$ and compute the mixing layer boundaries $y=\xi$ and $y=\xi+\eta$ using Eqs. (18) with $k=0$ (Model I). The channel thickness $h=1$ is chosen as the flow scale. The values of the empirical parameters are $\sigma=0.15$, $\kappa=6$. According to [8] the friction coefficient of the two-dimensional rectangular duct flow can be found as

Therefore, for flows with the Reynolds numbers of the order of $10^{4}$ considered in [23], the friction coefficient is $c_{f}=0.007$.

Boundaries of the mixing layer $y=\xi$ and $y=\xi+\eta$ obtained by Eqs. (18), as well as a symmetric about the axis $y=0$ solution corresponding to the conditions $\zeta_{0}=4.8\,h$, $w_{0}=1$ and $u_{0}=0.001$, are shown in Fig. 3. Solid curves correspond to the dissipation parameter $\kappa=6$, dashed lines — $\kappa=8$. The results of these calculations are overlaid on the grey-scaled picture of a turbulent jet presented in [23] (Fig. 3, d). It should be noted that the region of intense mixing mainly propagates into the layer of the fluid at rest and the distance between the internal boundaries of the mixing layers does not decrease until the stability of the jet is lost. As can be seen from Fig. 3, Eqs. (18) make it possible to determine the boundaries of intense mixing region in a Hele–Show jet flow without time-consuming calculations on the basis of 2D or 3D equations.

Let a solution of Eqs. (18) be constructed in the framework of Model I or II. Below we show that this solution can be applied to determine the velocity profile and Reynolds shear stress across the mixing layer. For this goal, we use stationary equations (15) taking into account the intermittency of the flow (irregularity of the layer boundaries in a turbulent flow), we set (see [15] and [18], Ch. 7)

The functions $q(x)$, $\eta(x)$ and $P=p(x)$ are considered as known ones. They are determined by the solution of stationary Eqs. (4). Under the above-mentioned assumptions, Eqs. (15) for steady flows take the form

To analyse this system, it is convenient to pass to the Mises variables $(x,\psi)$, where $\psi$ is the stream function associated with the velocity field by relations $\overline{U}=\psi_{y}$, $\overline{V}=-\psi_{x}$.

We denote

and take into account that

(derivatives of the functions $\overline{V}$ and $\hat{q}$ are calculated in the similar way). Then Eqs. (21) in the variables $(x,\psi)$ are transformed to the semi-linear hyperbolic system with respect to the variable $x$:

The characteristics of system (22) are given by the equations $d\psi/dx=\pm\sigma q$. The problem is similar to that considered in [18, 5] for a plane-parallel shear flow in a channel of finite thickness (flow with a pressure gradient). This is formulated for Eqs. (22) in the form of the Goursat problem with data on the characteristics. The difference consists only in the modification of the right-hand side (22) related to the consideration of friction. Solution of Eqs. (22) is constructed numerically using a Godunov-type scheme. Then, the velocity profile $\overline{U}(x,y)$ is obtained as a result of integration of the equation $y_{\psi}=1/\tilde{u}(x,\psi)$ and the transition to the original variables $(x,y)$.

The results of experimental measurements of the velocity profile and Reynolds shear stress are given in [22]. This work considers a turbulent jet flow between two parallel plates (separated by a distance $h$) issuing from a rectangular nozzle (with a span-wise size $d$) in the stream-wise direction $x$. Fig. 3 in [22] presents the measurement results for jet flows with ${\rm Re}=10^{4}$, $2\cdot 10^{4}$ and $5\cdot 10^{4}$ obtained for $d/h=8.75$ at the cross-section $(x-x_{0})/h=75$. Here $x_{0}$ is the virtual origin of the jet ($x_{0}=0$ for ${\rm Re}=5\cdot 10^{4}$). The velocity profiles and Reynolds stresses presented in [22] are normalized to $U_{C}$ and $U_{C}^{2}$, respectively. Here $U_{C}$ is the centreline velocity in the far jet field taken in the form [10]

with $c_{f}=0.0075$ and $c_{t}=0.007$.

To carry out the calculations based on Eqs. (18) (Model I) and then using Eqs. (22), we choose $h=1$, $Y=100\,h$, $\xi_{0}=d/2=4.375$, $u_{0}=1$ and $w_{0}=0.025$. The empirical parameters are as follows: $\sigma=0.15$, $\kappa=8$ and $c_{f}=0.0075$. We also assume $x_{0}=0$. The calculation results (functions $\overline{U}(x,y)$ and $-\sigma q(x)\hat{q}(x,y)$ normalized to the jet velocity $u(x)$ and $u^{2}(x)$) are shown in Fig. 4 by solid curves at the cross-section $x/h=75$. The value $y_{0.5}$ is chosen so that $\overline{U}=0.5\,u$ for $y=y_{0.5}$. As can be seen from Fig. 4, the calculation results for the proposed models are in a good agreement with experimental data [22]. We note that outside the mixing layer, $\hat{q}=0$ and $\overline{U}$ coincides with the external velocity $u$ (for $y<\xi$) or $w$ (for $y>\xi+\eta$). Since we take $w_{0}=0.025>0$ to avoid singularity in the governing equations, on the right border of Fig. 4 (a) the calculated velocity profile is slightly higher than the experimental values. The condition $w>0$ is fulfilled everywhere in the flow. The Reynolds shear stress practically coincides with the data from [22] for ${\rm Re}=2\cdot 10^{4}$ in the region $y/y_{0.5}>0.75$. The differences for $y/y_{0.5}<0.5$ can be associated with meandering and penetration of vortices into the central part of the jet.

Fig. 6 presents the stream-wise velocity in the centre of the jet. It follows from Fig. 6 that the centreline velocity in the far-field $U_{C}(x)$ given by the empirical formula (23) has the similar behaviour as the velocity $u(x)$ calculated by Eqs. (18) for $x/h>50$. At the considered cross-section $x/h=75$ these values almost coincide. Thus, using the solution of Eqs. (18), we can obtain the transverse distribution of the velocity profile and Reynolds shear stress in the mixing layer on the basis of hyperbolic Eqs. (22).