4.1. Initial helical tube

We search for an initial flow which is helically symmetric along the $z$-direction, of spatial period $2{\rm \pi} L$ and characterized by a compact vorticity field, such that $\boldsymbol {\omega }$ tends at least exponentially to zero away from a given helical line.

Let us first define this helical line. As illustrated in figure 1, the line intersects the $z=0$ plane (thereafter called $\varPi _0$) at a point ${\mathrm {A}\!^{\star }}$ defined by its cylindrical coordinates $r_{\mathrm {A}\!^{\star }}$ and $\theta _{\mathrm {A}\!^{\star }}$. The angle is set to $\theta _{\mathrm {A}\!^{\star }}=0$ without loss of generality. This helical line is thus located at

(4.1)\begin{equation} \boldsymbol{x}(\theta_s) =r_{\mathrm{A}\!^{{\star}}} \boldsymbol{e}_r(\theta_s) + {L} \theta_s \boldsymbol{e}_z,\quad \theta_s \in\mathbb{R}. \end{equation}

The plane perpendicular to this helical line at point ${\mathrm {A}\!^{\star }}$ is called the plane $\varPi ^\star _\perp$. The unit vector $\boldsymbol {e}_{{\scriptscriptstyle B}_{\mathrm {A}\!^{\star }}}$ defines the upward normal vector to plane $\varPi ^\star _\perp$

(4.2)\begin{equation} \boldsymbol{e}_{{\scriptscriptstyle B}_{\mathrm{A}\!^{{\star}}}}=\alpha(r_{\mathrm{A}\!^{{\star}}})\left[\boldsymbol{e}_{z}+\frac{r_{\mathrm{A}\!^{{\star}}}}{L} \boldsymbol{e}_{\theta}(0)\right]. \end{equation}

In this plane, a Cartesian basis $(\boldsymbol {e}_{r_{\mathrm {A}\!^{\star }}}, \boldsymbol {e}_{\varphi _{\mathrm {A}\!^{\star }}})$ is defined as well

(4.3a,b)\begin{equation} \boldsymbol{e}_{r_{\mathrm{A}\!^{{\star}}}} \equiv \boldsymbol{e}_{r}(0) ,\quad \boldsymbol{e}_{\varphi_{\mathrm{A}\!^{{\star}}}} \equiv \alpha(r_{\mathrm{A}\!^{{\star}}}) \left[\boldsymbol{e}_{\theta}(0)-\frac{r_{\mathrm{A}\!^{{\star}}}}{L} \boldsymbol{e}_{z}\right], \end{equation}

associated with local polar coordinates $(\rho, \psi )$ centred at point $A^\star$ with local basis $(\boldsymbol {e}_{\rho }, \boldsymbol {e}_{\psi })$

(4.4a,b)\begin{equation} \boldsymbol{e}_{\rho} =\cos \psi\; \boldsymbol{e}_{r_{\mathrm{A}\!^{{\star}}}}+\sin \psi\; \boldsymbol{e}_{\varphi_{\mathrm{A}\!^{{\star}}}} ,\quad \boldsymbol{e}_{\psi} =\cos \psi\; \boldsymbol{e}_{\varphi_{\mathrm{A}\!^{{\star}}}} -\sin \psi\; \boldsymbol{e}_{r_{\mathrm{A}\!^{{\star}}}} . \end{equation}

Figure 1.

Geometry for building the initial condition.

Practically, the initial condition is a thin-core vortex along the helical line: one chooses a compact profile for ${\omega _{{\scriptscriptstyle B}}}(\rho,\psi )$ and for ${u_{{\scriptscriptstyle H}}}(\rho,\psi )$ around point $A^\star$. This ensures that $\omega _\varphi$ and $\omega _r$ are also compact, a crucial issue. For instance, one assumes axisymmetric profiles similar to the Batchelor vortex

(4.5a,b)\begin{equation} \omega_{\scriptscriptstyle B}(\rho) = {\omega^\star_{\scriptscriptstyle B}} \exp\left( -\frac{\rho^2}{a^{\star2}}\right),\quad {u_{{\scriptscriptstyle H}}} (\rho )= {u^\star_{\scriptscriptstyle H}} \exp\left( -\frac{\rho^2}{a^{\star2}}\right), \end{equation}

where quantities ${\omega ^\star _{\scriptscriptstyle B}}$, ${u^\star _{\scriptscriptstyle H}}$ and core size $a^\star$ are three dimensional constants. From these profiles in the plane $\varPi ^\star _{\perp }$ as well as the values of the helix radius $r_{\mathrm {A}\!^{\star }}$ and reduced pitch $L$, it is possible to compute the initial fields $\omega _{\scriptscriptstyle B} (r,\varphi )$ and ${u_{{\scriptscriptstyle H}}}(r,\varphi )$ on the $\varPi _0$ plane used in simulations. To do so, we establish the connections between $(r,\varphi )$ given by a point in $\varPi _0$ and the polar radius $\rho$ of the corresponding point in $\varPi _\perp ^\star$, using invariances along the lines $\varphi = \theta - z/L=\mathrm {cst}$ (details of such a procedure can be found in Selçuk et al. Reference Selçuk, Delbende and Rossi2017, Appendix A). Let us consider an initial condition called case A characterized by circulation $\varGamma =1$, initial helix radius $r_{\mathrm {A}\!^{\star }}=1$, reduced pitch $L=0.3$, initial core size $a^\star =0.1$ and initial axial flow intensity ${u^\star _{\scriptscriptstyle H}}=1$. This yields $\omega _{\scriptscriptstyle B}^\star = 32.68$. On the $\varPi ^\star _{\perp }$ plane, the fields are hence axisymmetric and Gaussian (figure 2a) whereas in the $\varPi _0$ plane, they take a bean-like shape (figure 2b).

Figure 2.

Contours of ${\omega _{{\scriptscriptstyle B}}}/{\omega _{{\scriptscriptstyle B}}}_{{max}}$ for case A (see text) (a) in plane $\varPi ^\star _\perp$ and (b) in plane $\varPi _0$. The dashed circle in (a) represents the initial vortex core size $\rho =a^\star$. Contour levels $\omega _p$ are given by $\sqrt {-\log (\omega _p)}= \frac {1}{10}p\sqrt {\log (10^{3})}$, with $p=1,\ldots,10$.

4.2. Reaching a viscous helical quasi-equilibrium state

The Navier–Stokes equations with helical symmetry are integrated over time starting from the previous initial condition (case A). The numerical code uses a polar grid in a circular sub-domain of plane $\varPi _0$; it implements the generalized $\varPsi$–$\omega _{\scriptscriptstyle B}$–${u_{{\scriptscriptstyle H}}}$ formulation described in § 2, with second-order finite differences along $r$, Fourier expansions along $\varphi$, second-order time discretization and fully implicit viscous terms (Delbende, Rossi & Daube Reference Delbende, Rossi and Daube2012). Kinematic viscosity is set to $\nu =1.59\times 10^{-5}$, so that the Reynolds number $\mbox {Re}\equiv \varGamma /(2{\rm \pi} \nu )$ is $10^{4}$. It can be observed (see figure 3 for $\omega _{\scriptscriptstyle B}$, and 4 for ${u_{{\scriptscriptstyle H}}}$) that filaments are emitted during the early relaxation stage (see $t=4$), and then quickly dissipate ($t=10$) while the vortex becomes quasi-steady ($t=20$).

Figure 3.

Case A at $t=4$, $10$ and $20$: contours of ${\omega _{{\scriptscriptstyle B}}}/{\omega _{{\scriptscriptstyle B}}}_{{max}}(t)$ in plane $\varPi _\perp (t)$. Contours levels are as defined in figure 2.

Figure 4.

Case A at $t=0$, $4$, $10$ and $20$: contours of ${u_{{\scriptscriptstyle H}}}/{u_{{\scriptscriptstyle H}}}_{{max}}$ in plane $\varPi _\perp (t)$. Contours levels are defined in the same way as in figure 2 for ${\omega _{{\scriptscriptstyle B}}}/{\omega _{{\scriptscriptstyle B}}}_{{max}}$.

We check how the DNS solution approaches an Euler equilibrium using conditions (3.15a,b). First, we determine the translating reference frame: velocity $U_0(t)$ is computed by the best correlation of the vorticity field between successive time steps. Second, the streamfunction $\varPsi ^{{(T)}}$ in the translating frame is computed using (3.6). From now on, we only consider streamfunction $\varPsi ^{{(T)}}$ which will be simply denoted as $\varPsi$. Figure 5(a,d,g) displays scatterplots of $u_{\scriptscriptstyle H}$ vs $\varPsi$ at several times for case A. As time increases, it is observed that the points converge to a single continuous curve, supporting the existence of the function $g$ such that $u_{\scriptscriptstyle H}=g(\varPsi )$. The second condition (3.15a,b) requires the derivative $\mathrm {d} u_{\scriptscriptstyle H} / \mathrm {d} \varPsi$, a quantity which is defined only if a steady state is reached. A procedure is proposed in Appendix A to provide an estimate value for this derivative (see formula (A4)) during the whole time evolution, that converges to the exact value when the quasi-steady state is reached. In figure 5(b,e,h), scatterplots of quantity $\varpi$ vs $\varPsi$ are found to support the existence of a function $f$ such that $\varpi =f(\varPsi )$.

Figure 5.

Case A: scatterplots at times $t=0$ (a–c), $t=2$ (d–f) and $t=20$ (g–i): (a,d,g) ${u_{{\scriptscriptstyle H}}}$ vs $\varPsi$ (the red curve is an estimate of $g(\varPsi )$ obtained via formula (A3)); (b,e,h) $\varpi$ vs $\varPsi$; (c, f,i) $\alpha \omega _{\scriptscriptstyle B}$ (black dots) and $\alpha ^2 \mathrm {d} (u_{\scriptscriptstyle H} + C)^2 / \mathrm {d} \varPsi$ (blue dots) vs $\varPsi$. Note that only grid points such that $\omega _{\scriptscriptstyle B}/\omega _{{\scriptscriptstyle B}\,max}>10^{-3}$ have been represented.

For case A, however, quantity $\alpha \omega _{\scriptscriptstyle B}$ dominates the value $\varpi$ (see figure 5c, f,i): when the jet component is weak, quantity $\varpi$ almost equals $\alpha {\omega _{{\scriptscriptstyle B}}}$, which visually becomes a function of $\varPsi$. We present a second simulation (case B) where the jet component is initially more intense: circulation $\varGamma =1$, initial helix radius $r_{\mathrm {A}\!^{\star }}=1$, reduced pitch $L=1.5$, initial core size $a^\star =0.2$ and initial axial flow intensity ${u^\star _{\scriptscriptstyle H}}=2$. This yields $\omega _{\scriptscriptstyle B}^\star = 9.52$. The Reynolds number is set to $\mbox {Re}=10^3$. As time evolves, velocity ${u_{{\scriptscriptstyle H}}}$ converges to a function of $\varPsi$ (see figures 6a,d,g and 6b,e,h). For case B, partial quantities $\alpha {\omega _{{\scriptscriptstyle B}}}$ and $\alpha ^2 \mathrm {d} (u_{\scriptscriptstyle H} + C)^2 / \mathrm {d} \varPsi$ have same orders of magnitude (see figure 6c, f,i) so that the corresponding points remain scattered even at late times. Quantity $\varpi =\alpha {\omega _{{\scriptscriptstyle B}}} - \alpha ^2 \mathrm {d} (u_{\scriptscriptstyle H} + C)^2 / \mathrm {d} \varPsi$, however, clearly converges towards a function of $\varPsi$ (see figure 6b,e,h). This fully confirms the new finding expressed by relations (3.15a,b).

Figure 6.

Same as figure 5 but for case B (see text) plotted at times $t=0$ (a–c), $t=50$ (d–f) and $t=100$ (g–i).

Figure 7 represents the relation reached in the quasi-equilibrium state between $\varpi$ and ${u_{{\scriptscriptstyle H}}}$. When axial vorticity dominates over axial velocity (figure 7a) there is a linear relationship between $\varpi \approx \alpha {\omega _{{\scriptscriptstyle B}}}$ and ${u_{{\scriptscriptstyle H}}}$. When both effects have comparable orders of magnitude, this relation becomes highly nonlinear (figure 7b). Finding an analytical expression for such relations remains an open issue.

Figure 7.

Quantity $\varpi$ as a function of ${u_{{\scriptscriptstyle H}}}$ at quasi-equilibrium for (a) case A (see figure 5) and (b) case B (see figure 6).

4.3. Determination of the centre of the helical vortex

The above results suggest that, in order to characterize the structure of the vortex, the most convenient way is to select the point of maximum $\varPsi$ as a centre for the vortex. When this quantity is available, one looks for the point A where $\varPsi$ reaches a maximum in the plane $\varPi _0$ (for most single helical vortices, this point is unique), and one defines the associated plane $\varPi _\perp$ in a similar way as $\varPi _\perp ^\star$ is defined from ${\mathrm {A}\!^{\star }}$ (see above in § 4.1 and figure 1). Since the field $\varPsi$ depends on $(r,\varphi )$, there is a one-to-one relation between its values in the plane $\varPi _\perp$ and in the plane $\varPi _0$. Consequently, it also reaches a maximum in the $\varPi _\perp$ plane at point A. Other properties can be also useful to locate this central point. Since $r u_r = \partial _\varphi \varPsi$ and $u_\varphi = -\alpha (r)\partial _r\varPsi$, velocity at point A is such that $u_r=u_{\varphi }=0$ and is thus tangent to a helical line, i.e. $\boldsymbol {u}(\textrm {A})={u_{{\scriptscriptstyle B}}} {\boldsymbol {e}_{\scriptscriptstyle B}}$. Furthermore, if the solution is a quasi-steady equilibrium, then relation $u_{\scriptscriptstyle H}+C=g(\varPsi )$ holds and, necessarily, velocity $u_{\scriptscriptstyle H}$ has an extremum at point A because, at point A

(4.6a,b)\begin{equation} \partial_\varphi {u_{{\scriptscriptstyle H}}} = \frac{\mathrm{d} {u_{{\scriptscriptstyle H}}}}{\mathrm{d} \varPsi} \partial_\varphi \varPsi=0, \quad \partial_r {u_{{\scriptscriptstyle H}}} = \frac{\mathrm{d} {u_{{\scriptscriptstyle H}}}}{\mathrm{d} \varPsi} \partial_r\varPsi =0. \end{equation}

Since $\partial _\varphi {u_{{\scriptscriptstyle H}}}=r\omega _r$ and $\partial _r {u_{{\scriptscriptstyle H}}}=-\omega _{\varphi }/\alpha$, vorticity at point A is such that $\omega _r=\omega _{\varphi }=0$ and is thus tangent to a helical line, i.e. $\boldsymbol {\omega }(\textrm {A})={\omega _{{\scriptscriptstyle B}}}{\boldsymbol {e}_{\scriptscriptstyle B}}$. The points where $\alpha \omega _{\scriptscriptstyle B}$ and ${u_{{\scriptscriptstyle H}}}$ reach an extremum generally do not coincide. Only in the particular case when $u_{\scriptscriptstyle H}=0$, quantity $\alpha \omega _{\scriptscriptstyle B}$ is also a function of $\varPsi$ and thus reaches an extremum at the same point. This is illustrated in table 1: for case A with relatively small axial flow, these radial locations differ by less than 0.1 %, while they differ by 3 % for case B because of the larger axial flow.

Table 1.

Radial location of the maximum for different quantities, as measured in the DNS for cases A and B.

4.4. Quasi-equilibrium helical vortex with prescribed parameter values

The viscous quasi-equilibrium reached after a given simulation time depends on five dimensional parameters used to initiate the simulation as explained in § 4.1, namely circulation $\varGamma$, helix radius $r_{\mathrm {A}\!^{\star }}$, helical pitch $2{\rm \pi} L$, vortex core size $a^\star$, axial flow intensity ${u^\star _{\scriptscriptstyle H}}$. The vorticity amplitude ${\omega ^\star _{\scriptscriptstyle B}}$ is a function of circulation $\varGamma$ and velocity ${u^\star _{\scriptscriptstyle H}}$. The goal is to generate a quasi-equilibrium to be used as a base flow in an instability study. Hence, we would like to prescribe the final state obtained through the DNS rather than the initial state (4.5a,b). In order to characterize the final state, we need to define five such parameters. The centre A of the vortex can be determined as explained in the previous section, yielding the helix radius $r_{A}$. The helical pitch is a fixed value $2{\rm \pi} L$ and circulation $\varGamma$ is given by

(4.7)\begin{equation} \varGamma \equiv \iint_{\varPi_0} \omega_{z} \, r\, \mathrm{d} r \,\mathrm{d} \theta. \end{equation}

The core size can be obtained through an integral in the plane $\varPi _\perp$ associated with point A

(4.8)\begin{equation} a^{2} \equiv \frac{1}{\varGamma} \iint_{\varPi_{{\perp}}} \boldsymbol{\omega} \boldsymbol{\cdot} \boldsymbol{e}_{{\scriptscriptstyle B}_{A}} \left(\boldsymbol{x}-\boldsymbol{x}_{A} \right)^{2} d S_{{\perp}}, \end{equation}

where $\boldsymbol {e}_{{\scriptscriptstyle B}_{A}} \equiv \boldsymbol {e}_{{\scriptscriptstyle B}}(\textrm {A})$ and $\boldsymbol {x}_{A}\equiv \boldsymbol {x}(\textrm {A})$. The last parameter is the axial velocity parameter

(4.9)\begin{equation} W_B \equiv \frac{ \alpha_{A}^2}{{\rm \pi} a^2} \iint_{\varPi_{0}} {u_{{\scriptscriptstyle H}}} \,\mathrm{d} S \quad \text{with } \alpha_{A} \equiv \alpha(r_{A}). \end{equation}

An inverse swirl number can be defined as the ratio of the axial velocity parameter $W_B$ and the typical azimuthal velocity in the core $\varGamma /(2{\rm \pi} a)$

(4.10)\begin{equation} \bar W_B \equiv \frac{ 2{\rm \pi} a\, W_B}{ \varGamma }. \end{equation}

If the profiles are close to a thin curved Batchelor vortex (4.5a,b), the core size $a$ is almost equal to $a^\star$ and the quantity $W_B$ is close to the amplitude of helical velocity $\alpha _{A}{u^\star _{\scriptscriptstyle H}}$. Indeed, for a thin core, $\mathrm {d} S_{\perp } \approx \alpha \,\mathrm {d} S \approx \alpha _{A}\, \mathrm {d} S$ and

(4.11)\begin{align} W_B \approx \frac{\alpha_{A}}{{\rm \pi} a^2} \iint_{\varPi_{{\perp}}} {u_{{\scriptscriptstyle H}}} \mathrm{d} S_\perp=\frac{\alpha_{A}}{{\rm \pi} a^2} {{u^\star_{\scriptscriptstyle H}} } 2 {\rm \pi}\int_{0}^{\infty} \mathrm{e}^{-(\rho/{a^\star})^{2}} \rho \,\mathrm{d} \rho =\frac{\alpha_{A}}{{\rm \pi} a^2}{u^\star_{\scriptscriptstyle H}} {\rm \pi}a^{\star2}\approx \alpha_{A} {u^\star_{\scriptscriptstyle H}}. \end{align}

Another possibility for obtaining a value for $a$ and $W_B$ is to use a Gaussian fit of the axisymmetric component of fields ${\omega _{{\scriptscriptstyle B}}}$ and/or ${u_{{\scriptscriptstyle H}}}$.

From now on in § 4, all variables are put in dimensionless form using as characteristic dimensional length scale $r_{{A}}$ and time scale $r^2_{{A}}/\varGamma$ of the final state. Since the dimensionless radius at point A and dimensionless circulation are equal to one, the quasi-equilibrium state depends only on two dimensionless lengths $L$ and $a$, and on the inverse swirl number $\bar W_B$. In order to simulate the dynamics of the vortex towards the quasi-equilibrium, we also need to select a value for the Reynolds number $\mbox {Re}\equiv \varGamma /(2{\rm \pi} \nu )$ and a simulation time denoted as $T_{sim}$. In the two-dimensional framework, it is known that the characteristic time necessary for the vortex to reach a quasi-equilibrium is of order $a^2\mbox {Re}^{{1}/{3}}$ with a pre-factor of order 40 (Bernoff & Lingevitch Reference Bernoff and Lingevitch1994). We here select the final simulation time $T_{sim}\sim 60a^2 \mbox {Re}^{{1}/{3}}$ and, in practice, check a posteriori that it is large enough to reach a quasi-equilibrium. For $\mbox {Re}=10^4$, for instance, we selected $T_{sim}=20, 40, 100$ for $a=0.11, 0.174$ and $0.3$, respectively.

The reduced pitch $L$ and circulation $\varGamma$ of the initial condition are identical to those of the quasi-equilibrium we are looking for (hereafter called the final state). By contrast, initial parameters $r_{\mathrm {A}\!^{\star }}$, $a^\star$ and ${u^\star _{\scriptscriptstyle H}}$ differ from the final quantities $r_{A}=1$, $a$ and $W_B/\alpha _{{A}}$. We use an iterative procedure to determine these three unknown parameters. First, we select guess values for $a^\star$ and $r_{\mathrm {A}\!^{\star }}$, so that we link the two initial parameters ${u^\star _{\scriptscriptstyle H}}$ and ${\omega ^\star _{\scriptscriptstyle B}}$ to parameters ${a}$, $L$ and $\bar W_B$. This connection is obtained via two ‘conservation’ laws (B6) and (B12) obtained in Appendix B. They read in dimensionless form

(4.12a,b)\begin{equation} \mathcal{H}_1(t)-\frac{2}{L}\mathcal{H}_2(t)=1,\quad \mathcal{H}_4(t) = \mathcal{H}_4(0) -\frac{2t}{L\,\mbox{Re}}, \end{equation}

where

(4.13a–c)\begin{equation} \mathcal{H}_1 \equiv \iint_{\varPi_0} \alpha {\omega_{{\scriptscriptstyle B}}} \, \mathrm{d} S,\quad \mathcal{H}_2\equiv \iint_{\varPi_0} \alpha^4 {u_{{\scriptscriptstyle H}}} \, \mathrm{d} S~~~\mbox{and}~~~ \mathcal{H}_4 \equiv \iint_{\varPi_0} {u_{{\scriptscriptstyle H}}} \, \mathrm{d} S. \end{equation}

Let us introduce the quantities

(4.14)\begin{equation} {I}_p^\star (a^\star,r_{\mathrm{A}\!^{{\star}}} ,L) \equiv \iint_{\varPi_0} \alpha^p \exp(-(\rho/a^\star)^2) \,r\,\mathrm{d} r\,\mathrm{d}\theta, \end{equation}

with $p$ positive integers. Such quantities depend on $L$ via the quantity $\rho$, itself depending on $r$, $\theta$, $r_{\mathrm {A}\!^{\star }}$ and $L$ (see figure 1). Inserting the initial condition (4.5a,b) into (4.12a,b) allows one to recast variables $\varGamma (t=0)$ and $\mathcal {H}_4(t=T_{sim})$ as

(4.15a,b)\begin{equation} {\omega^\star_{\scriptscriptstyle B}} I_1^\star{-} \frac{2}{L}{u^\star_{\scriptscriptstyle H}} I_4^\star=1,\quad \mathcal{H}_4(T_{sim}) = {u^\star_{\scriptscriptstyle H}} I_0^\star{-} \displaystyle\frac{2T_{sim} }{ L\, \mbox{Re}}. \end{equation}

Introducing relation (4.10) leads to determination of ${u^\star _{\scriptscriptstyle H}}$ and ${\omega ^\star _{\scriptscriptstyle B}}$

(4.16a,b)\begin{equation} {u^\star_{\scriptscriptstyle H}} = \frac{1}{I_0^\star}\left[\frac{ \bar W_B a}{2\alpha_{A}^2 } + \displaystyle\frac{2 T_{sim}}{L\, \mbox{Re}}\right],\quad {\omega^\star_{\scriptscriptstyle B}} = \frac{1}{I_1^\star}\left(1+\frac{2}{L}{u^\star_{\scriptscriptstyle H}} I_4^\star\right). \end{equation}

Of course, since $a^\star$ and $r_{A}^\star$ have only been guessed at this stage, the DNS initiated by the above procedure does not lead at once to the prescribed values of $a$ and $r_{A}=1$ at $t=T_{sim}$. We have to use a standard iterative procedure to search for the proper pair $(a^\star,r_{\mathrm {A}\!^{\star }})$ leading, after a simulation time of duration $T_{sim}$, to the prescribed values $(a, 1)$.

4.5. Effect of Reynolds number

The ability to prescribe the final state parameters presented in the previous section makes it possible to use several distinct Reynolds numbers to try to achieve the same quasi-equilibrium state. Figure 8 shows the influence of $\mbox {Re}$ on the two curves ${u_{{\scriptscriptstyle H}}}=g(\varPsi )$ and $\varpi =f(\varPsi )$ obtained numerically. The cases $\mbox {Re}=5000$ and $10^4$ cannot be visually distinguished, yet at lower Reynolds numbers (see $\mbox {Re}=1000$), the equilibrium curves shift more significantly as the value of $\varPsi _{max}$ also decreases. In the following, we always will adopt a sufficiently large Reynolds number (typically $\mbox {Re}=10^4$) to make sure that the quasi-equilibrium state obtained is $\mbox {Re}$-independent.

Figure 8.

Relations (a) ${u_{{\scriptscriptstyle H}}}=g(\varPsi )$ and (b) $\varpi =f(\varPsi )$ for quasi-equilibria obtained at different Reynolds numbers: $\mbox {Re}=500$ (red dots, $T_{sim}=4$), $\mbox {Re}=10^3$ (blue dots, $T_{sim}=4$), $\mbox {Re}=5\times 10^3$ (magenta dots, $T_{sim}=10$) and $\mbox {Re}=10^4$ (black solid line, $T_{sim}=20$). The dimensionless parameters of these states are $L=0.3$, $a=0.11$ and $\bar W_B = 0.2$.

4.6. Helical vortex cases

Using the procedure presented, we generate several base flow configurations with various values of parameters $L$, $a$ and $\bar W_B$. For a helical vortex of reduced pitch $L$, curvature is given by $\kappa =1/(1+L^2)$ and torsion by $\tau =L/(1+L^2)$. As plotted in figure 9, curvature ${\kappa }$ decreases from $1$ to $0$ as $L$ increases from $0$, while torsion $\tau$ first increases from $0$ to $0.5$ (reached for $L=1$), and then slowly decreases to $0$. The asymptotic theory presented in § 5.1 uses parameter $\varepsilon \equiv a\kappa$ as a small parameter to introduce curvature effects on a straight axisymmetric vortex. A given value of $\varepsilon$ can be achieved from different pair values $(L, a)$, thus corresponding to different curvature and torsion levels. We choose five helical vortex states, the parameters of which are presented in table 2.

Figure 9.

Non-dimensional curvature $\kappa$ (solid line) and torsion $\tau$ (black dashed line) as functions of the reduced pitch $L$. Pitches $L$ and $1/L$ lead to the same value of the torsion but to different curvature levels, as illustrated for $L=0.7$ and $1/0.7\approx 1.43$.

Table 2.

Parameters $(\bar W_B, L, a)$ of the DNS cases investigated, with corresponding values of parameter $\varepsilon$, curvature $\kappa$ and torsion $\tau$.

The first three states correspond to the same value of $\varepsilon =0.1$ but to different geometrical parameter values. State BS1 is characterized by relatively small values of the helical pitch and core size, leading to high curvature and low torsion. State BS2 is characterized by larger values of the helical pitch and core size, leading to low curvature but near-maximum torsion. State BS3 is the same as BS1 except that the inverse swirl parameter $\bar W_B$ is twice as much that of BS1 (or BS2). State BS4 is a state with a larger value $\varepsilon =0.16$ obtained through a larger core size. State BS5 is similar to BS4, yet with reversed axial flow (hence negative inverse swirl).

5.1. Asymptotic expansion

The leading-order term $\bar {\boldsymbol {u}}^{(0)}$ is assumed to be a Batchelor vortex. Indeed, this vortex model is known to be a self-similar attracting viscous solution for two-dimensional vortices with axial flow. Such a profile is also found experimentally and numerically in the core of helical vortices (Ali & Abid Reference Ali and Abid2014; Selçuk et al. Reference Selçuk, Delbende and Rossi2017; Okulov et al. Reference Okulov, Kabardin, Mikkelsen, Naumov and Sørensen2019; Shtork et al. Reference Shtork, Gesheva, Kuibin, Okulov and Alekseenko2020). This means that the velocity distribution in the plane $\varPi _{\perp }$ is

(5.4a–c)\begin{equation} \bar u_{\rho}^{(0)} = 0,\quad \bar u_{\psi}^{(0)} = \frac{1- \exp(-\bar\rho^2)}{\bar\rho},\quad \bar u_{b}^{(0)} = \bar W_0 \exp(-\bar\rho^2) + \bar W_{00} . \end{equation}

The parameter $\bar W_0$ indicates the jet strength in the axial direction and the parameter $\bar W_{00}$ is a correction at leading order deduced from the background velocity outside the vortex. Based on these velocity fields, the vorticity field at zeroth order in $\varepsilon$ is

(5.5)\begin{equation} \bar{\boldsymbol{\omega}}^{(0)}= \bar{\omega_{{\scriptscriptstyle B}}}^{(0)} \boldsymbol{e}_b + \bar\omega_{\psi}^{(0)} \boldsymbol{e}_\psi, \end{equation}

with

(5.6a,b)\begin{equation} \bar{\omega_{{\scriptscriptstyle B}}}^{(0)} = 2 \exp{(- \bar\rho^2)},\quad \bar\omega_{\psi}^{(0)} = 2 \bar W_0 \bar\rho \exp{(- \bar\rho^2)}. \end{equation}

For any point M in plane $\varPi _\perp (\textrm {A})$ except A itself, the velocity component $u_{\scriptscriptstyle B}$ along direction $\boldsymbol {e}_{\scriptscriptstyle B}(\textrm {M})$ is different from the velocity component $u_b$ along the direction $\boldsymbol {e}_b \equiv \boldsymbol {e}_{\scriptscriptstyle B}(\textrm {A})$ orthogonal to the plane. Quantity $u_{\scriptscriptstyle B}/\alpha$ is given by

(5.7)\begin{equation} \frac{u_{\scriptscriptstyle B}}{\alpha}= \frac{u_\rho}{\alpha} \;\boldsymbol{e}_{\rho} \boldsymbol{\cdot} \boldsymbol{e}_{\scriptscriptstyle B} + \frac{u_\psi}{\alpha} \;\boldsymbol{e}_{\psi} \boldsymbol{\cdot} \boldsymbol{e}_{\scriptscriptstyle B} + \frac{u_b}{\alpha} \;\boldsymbol{e}_{b} \boldsymbol{\cdot} \boldsymbol{e}_{\scriptscriptstyle B}, \end{equation}

where $\alpha \equiv \alpha (\textrm {M})$. Using relations (deduced from relations (C3)–(C11) in Appendix C)

(5.8a–c)\begin{equation} \boldsymbol{e}_\rho \boldsymbol{\cdot} \boldsymbol{e}_{\scriptscriptstyle B} = 0,\quad \boldsymbol{e}_{\psi}\boldsymbol{\cdot} \boldsymbol{e}_{\scriptscriptstyle B} = \alpha \alpha_{A} \frac{\rho}{L},\quad \boldsymbol{e}_b \boldsymbol{\cdot} \boldsymbol{e}_{\scriptscriptstyle B} = \alpha \alpha_{A}\left(1 + \frac { r^2_{A}}{L^2} -\frac{r_{A}}{L}\frac{\rho}{L} \cos \psi \right), \end{equation}

one obtains, at zeroth order in $\varepsilon$

(5.9)\begin{align} \frac{\bar u_{\scriptscriptstyle B}^{(0)}}{\alpha} = \frac{\bar u^{(0)}_b(\bar\rho)}{\alpha_{A}} + {\bar u^{(0)}_\psi(\bar\rho)}{\alpha_{A}}\frac{\bar \rho}{\bar L} = \frac{1}{\alpha_{A}}[\bar W_0 \exp{(- \bar\rho^2)} +\bar W_{00}]+\frac{\alpha_{A}}{\bar L}[1- \exp{(- \bar\rho^2)}]. \end{align}

Quantity $\bar W_0$ is, by construction, given by $\bar W_0 = \bar u^{(0)}_b(0) -\bar u^{(0)}_b(\infty )$. Using relation (5.9) twice at $\bar \rho =0$ and at the limit $\bar \rho \to \infty$, one deduces, at zeroth order in $\varepsilon$

(5.10)\begin{equation} \bar W_0 = \alpha_{A}\left(\frac{\bar u^{(0)}_{\scriptscriptstyle B}(0)}{\alpha_{A}}-\frac{\bar u^{(0)}_{\scriptscriptstyle B}(\infty)}{\alpha}\right) + \frac{\alpha_{A}^2}{\bar L} = \alpha_{A}\bar u^{(0)}_{\scriptscriptstyle H}(0) + \frac{\alpha_{A}^2}{\bar L}. \end{equation}

Quantity $\bar W_{00}$ can be obtained by matching the value reached by $\bar u_{\scriptscriptstyle B}^{(0)} / \alpha$ outside the vortex to the constant introduced in (3.10)

(5.11)\begin{equation} \bar C = \frac{2{\rm \pi} a C}{\varGamma} = \frac{1}{\bar L}+\bar U_z^{\infty{(T)}}. \end{equation}

As $\bar \rho \to \infty$, expression (5.9) tends to $\bar W_{00}/\alpha _{A}+\alpha _{A}/\bar L$. The matching eventually yields

(5.12)\begin{equation} \bar W_{00} = \alpha_{A}\left(\bar C - \frac{\alpha_{A}}{\bar L}\right) = \alpha_{A}\bar U_z^{\infty{(T)}} + \frac{\alpha_{A}(1 -\alpha_{A})}{\bar L}. \end{equation}

In Blanco-Rodríguez et al. (Reference Blanco-Rodríguez, Le Dizès, Selçuk, Delbende and Rossi2015), the simplified expression $\bar W_{00} \approx \alpha _A U_z^{\infty {(T)}}$ was used. Expression (5.12) shows a better agreement with the numerical results (see § 5.2).

The two constants $\bar W_0$ and $\bar W_{00}$ can thus be related to the parameter values of the base helical vortex using (5.10) and (5.12) with $\bar U_z^{\infty {{(T)}}}=-\bar U_0 = \bar \varOmega _0\bar L$

(5.13a,b)\begin{equation} \bar W_0 \approx \bar W_B + \frac{\alpha_{A}^2}{\bar L},\quad \bar W_{00} =\alpha_{A} \bar\varOmega_0\bar L + \frac{\alpha_{A}(1-\alpha_{A})}{\bar L} . \end{equation}

The evaluation of $\bar W_{00}$ requires the knowledge of $\bar \varOmega _0$. Although this latter quantity could be computed theoretically (Okulov Reference Okulov2004) as done by Brynjell-Rahkola & Henningson (Reference Brynjell-Rahkola and Henningson2020), we use here the value deduced from the DNS results.

Going now to order $\varepsilon$, the perturbation is purely dipolar (Blanco-Rodríguez et al. Reference Blanco-Rodríguez, Le Dizès, Selçuk, Delbende and Rossi2015)

(5.14a,b)\begin{gather} \bar u_{\rho}^{(1)} ={-} \frac{\bar\varPsi^{(1)}}{\bar\rho} \sin \psi,\quad \bar u_{\psi}^{(1)}= \left[ \bar\rho\bar u^{(0)}_\psi -\frac{\mathrm{d} \bar\varPsi^{(1)}}{\mathrm{d} \bar\rho} \right] \cos \psi, \end{gather}

(5.15)\begin{gather} \bar u_{b}^{(1)} = \left[\bar\rho\bar u^{(0)}_b - \frac{1}{\bar u^{(0)}_\psi}\frac{\mathrm{d}\bar u^{(0)}_b}{\mathrm{d}\bar\rho}\bar\varPsi^{(1)} \right] \cos \psi, \end{gather}

where the first-order streamfunction perturbation $\bar \varPsi ^{(1)}(\bar \rho )$ satisfies

(5.16)\begin{equation} \bar\varPsi^{(1)}= \bar\varPsi^{(1)}_0 + \bar W^2_0 \bar\varPsi^{(1)}_1 + \bar W_{00} \bar W_0 \bar\varPsi^{(1)}_2. \end{equation}

The three contributions $\bar \varPsi ^{(1)}_0$, $\bar \varPsi ^{(1)}_1$ and $\bar \varPsi ^{(1)}_2$ are solutions of

(5.17a–c)\begin{equation} \mathcal{L} \bar\varPsi^{(1)}_0 = \bar K_0,\quad \mathcal{L}\bar\varPsi^{(1)}_1 = \bar K_1,\quad \mathcal{L} \bar\varPsi^{(1)}_2 =\bar K_2, \end{equation}

where the linear operator $\mathcal {L} (\bullet )$ is defined as

(5.18)\begin{equation} \mathcal{L} ({\bullet}) \equiv \frac{ \partial^2 ({\bullet}) }{ \partial \bar\rho^2}+ \frac{1 }{\bar\rho} \frac{ \partial ({\bullet}) }{ \partial \bar\rho} - \left(\frac{1 }{\bar\rho^2} + H_0(\bar\rho)\right) ({\bullet}) \end{equation}

and

(5.19)\begin{equation} \left. \begin{aligned} \bar H_0(\bar\rho) & \equiv \frac{1 }{\bar u_{\psi}^{(0)} } \frac{ \mathrm{d} \bar{\omega_{{\scriptscriptstyle B}}}^{(0)} }{ \mathrm{d} \bar\rho}= \frac{-4\bar\rho^2}{\exp{(\bar\rho^2)} -1},\\ \bar K_0(\bar\rho) & \equiv 2 \bar\rho\bar\omega^{(0)}_b +\bar u_{\psi}^{(0)} = \frac{1- \exp{(- \bar\rho^2)}}{\bar\rho} + 4 \bar\rho \exp{(- \bar\rho^2)} ,\\ \bar K_1(\bar\rho) & \equiv{-} \frac{4 \bar\rho^3 }{\exp{( \bar\rho^2)} -1 } \exp{(- \bar\rho^2)},\quad \bar K_2(\bar\rho) \equiv{-}\frac{4 \bar\rho^3 }{\exp{( \bar\rho^2)} -1 } . \end{aligned} \right\} \end{equation}

The boundary conditions at $\bar \rho =0$ are

(5.20)\begin{equation} \bar\varPsi^{(1)}_{i}(0)=0,\quad \frac{ \partial \bar\varPsi^{(1)}_{i} }{ \partial \bar\rho}(0)=0,\quad i=0, 1, 2. \end{equation}

We hence compute the first-order streamfunction perturbation $\bar \varPsi ^{(1)}$ numerically. For large $\bar \rho$, $\bar \varPsi ^{(1)}(\rho )$ expands as (Fukumoto & Miyazaki Reference Fukumoto and Miyazaki1991)

(5.21)\begin{equation} \bar\varPsi^{(1)}(\bar\rho) \sim \frac{1}{2}\bar\rho\log \bar\rho + \bar\rho \bar A + O\left(\frac{1}{\bar\rho}\right), \end{equation}

where, for the Batchelor vortex profile, $\bar A$ reads as

(5.22)\begin{equation} \bar A = \frac{1}{4}[1-\bar W_0(\bar W_0+4 \bar W_{00})+ \gamma - \log 2],\quad \gamma\approx0.577. \end{equation}

This formula corrects the typo $2\bar W_{00}$ in (3.12) of Blanco-Rodríguez et al. (Reference Blanco-Rodríguez, Le Dizès, Selçuk, Delbende and Rossi2015).

5.2. Comparison with DNS results

Quasi-equilibrium states listed in table 2 are obtained as final states of simulations through the procedure explained in § 4.4. Figure 10(a) (respectively figure 10b) shows the axial vorticity in the plane $\varPi _\perp$ for BS1 (respectively BS2). In addition, the internal structure can be compared with the results obtained by the asymptotic theory presented in § 5.1. To do so, one needs not only to rewrite the DNS results in the reference frame and in the basis used in the theory, but also to use the same non-dimensionalization process. More precisely, this means:

1. (i) adopting the frame ${(T)}$ translating at velocity $U_0 \boldsymbol {e}_z$ with the vortex. This amounts to transforming two of the velocity components: $u_{\scriptscriptstyle B} \longrightarrow u_{\scriptscriptstyle B}-\alpha U_0$ and $u_\varphi \longrightarrow u_\varphi +\alpha U_0 r/L$;

2. (ii) projecting the velocity field onto the plane $\varPi _\perp (\textrm {A})$ in order to get in-plane components $(u_\rho, u_\psi )$ and deduce out-of-plane component $u_b$, as explained in Appendix C;

3. (iii) performing an azimuthal decomposition in $\varPi _\perp (\textrm {A})$ of such components using

   (5.23)\begin{equation} u_{\rho,\psi,b}^{(m)}(\rho_{i})=\frac{1}{N_{\psi}} \sum_{j=1}^{N_{\psi}} u_{\rho,\psi,b}(\rho_{i}, \psi_{j}) \mathrm{e}^{-\mathrm{i} m \psi_{j}}\quad (i=1,\ldots,N_\rho), \end{equation}
   where $(\rho _{i}, \psi _{j})$ are the locations of $N_\rho \times N_\psi$ nodes of a polar mesh in plane $\varPi _\perp$; and

4. (iv) making the various quantities dimensionless using the reference scales $a_0$ and $2{\rm \pi} a_0^2/\varGamma$ as for the theory.

Based on the above procedure, the parameters presented in table 2 correspond to the values of effective core size $a_0$, inverse swirl $\bar W_0$ and constant $\bar W_{00}$ reported in table 3. Each case has been analysed, showing the adequacy of the asymptotic theory up to relatively large values of $\varepsilon$. Comparisons between asymptotic theory and DNS of helical vortices are shown in figures 11 and 12 for BS1 ($\varepsilon =0.1$) and BS5 ($\varepsilon =0.16$), respectively. In each graph, the $m=0$ and $m=1$ contributions to $(u_\rho ^{(m)},u_\psi ^{(m)},u_b^{(m)})$ are plotted. Excellent agreement is found for the axisymmetric part. A fairly small mismatch subsists for the components of the $m=1$ contribution.

Figure 10.

(a) Contours of $\omega _{\scriptscriptstyle B}$ in plane $\varPi _\perp$ for the quasi-equilibrium state BS1. (b) Same for BS2, with the same spatial range represented. Contour levels are as in figure 2. The Reynolds number is $\mbox {Re} = 10^4$, and the simulation time $T_{sim}=20$ (respectively 100) for BS1 (respectively BS2).

Figure 11.

Case BS1: radial distributions of velocity components $(\bar u_\rho ^{(m)},\bar u_\psi ^{(m)},\bar u_b^{(m)})$ for $m=0$ (a–c) and $m=1$ (d–f). Blue solid line: DNS results; red dashed line: theory.

Figure 12.

Same as figure 11 for case BS5.

Table 3.

For the different helical vortex states under study: case parameters $a$ and $\bar W_B$; effective core size $a_0$ and inverse swirl $\bar W_{0}$ obtained by DNS; values of $\bar W_B+\alpha _A^2/\bar L$ and $\bar W_{00}$ obtained through (5.13a,b).