Vortex Dynamics of Oscillating Flows
We employ the method of multiple scales (twotiming) to analyse the vortex dynamics of inviscid, incompressible flows that oscillate in time. Consideration of distinguished limits for Euler's equation of hydrodynamics shows the existence of two main asymptotic models for the averaged flows: strong vortex dynamics (SVD) and weak vortex dynamics (WVD). In SVD the averaged vorticity is 'frozen' into the averaged velocity field. By contrast, in WVD the averaged vorticity is 'frozen' into the 'averaged velocity + drift'. The derivation of the WVD recovers the CraikLeibovich equation in a systematic and quite general manner. We show that the averaged equations and boundary conditions lead to an energytype integral, with implications for stability.
Keywords
1 Introduction
Oscillating flows represent an important aspect of classical fluid dynamics and appear in various applications in medicine, biophysics, geophysics, engineering, astrophysics and acoustics. The term "oscillating flow" usually means that the fluid motion under consideration possesses a dominant frequency $\sigma$, which can be maintained either by an oscillating boundary condition, by an oscillating external force, or by selfoscillations of a flow. All other motions in the oscillating flow are considered as "slow", with the related timescale $T_{{\rm slow}}\gg 1/\sigma$. The scale $T_{{\rm slow}}$ also can be related to a boundary condition, to an external force, or it may characterise a natural intrinsic motion of the fluid. A related powerful mathematical approach is the twotiming method (see, for example, [Nayfeh1973]; [Kevorkian and Cole1996]). In this paper we use this method, together with the idea of distinguished limits, to provide an elementary, systematic, and justifiable procedure following the ideas proposed in [Vladimirov2005], [Yudovich2006], [Vladimirov2008]; [Vladimirov2010]. The results and contents of this paper may be summarised as:

1.
The development of a new analytical approach for the description of oscillatingintime flows. The fluid is assumed to be inviscid and incompressible, with the oscillations introduced via the boundary conditions.

2.
The analysis of distinguished limits for the Euler equation shows the existence of two asymptotic models for the averaged flows: "strong" or "standard" vortex dynamics (SVD), and weak vortex dynamics (WVD), the latter described by the CraikLeibovich equation (CLE). The CLE was originally derived for the description of the Langmuir circulations generated by surface waves ([Craik1985]; [Leibovich1983]; [Thorpe2004]). The derivation of the CLE demonstrates the remarkable fact that the Reynolds stresses can be expressed solely in terms of the drift velocity.

3.
In SVD the averaged vorticity is frozen into the averaged velocity, whereas in WVD the averaged vorticity is frozen into the averaged velocity $+$ drift velocity. It is important that in WVD the drift velocity has the same order of magnitude as the averaged velocity. Our derivation of the WVD and CLE is technically simpler than previous derivations. The formulation of the problem in its natural generality shows that the area of applicability of the CLE is broader than previously recognised. In particular, we consider flow domains that are threedimensional and of arbitrary shape; the oscillations are timeperiodic, but their spatial structure is arbitrary. We have also derived the averaged boundary conditions that are valid at the average positions of the boundaries.

4.
The slow timescale is uniquely linked to the magnitude of the prescribed velocity field at the boundary. Naturally, the higher the amplitude of velocity, the shorter the slow timescale.

5.
The WVD and CLE contain the drift velocity. The drift usually appears as the average velocity of Lagrangian particles (see [Stokes1847]; [Lamb1932]; [Batchelor1967]). In our consideration, drift velocity appears naturally as the result of an Eulerian averaging of the related PDEs without directly addressing the motions of particles.

6.
The CLE leads to an energytype integral for the averaged flows, which allows us to consider "Arnoldtype" results, such as the generalized "isovorticity conditions", the energy variational principle, the first and second variation of energy, and several (nonlinear and/or linear) stability criteria for averaged flows.
2 TwoTiming Problem and Distinguished Limits
We study the motion of a homogeneous inviscid incompressible fluid in a timedependent threedimensional domain $Q(t)$ with oscillating boundary $\partial Q(t)$ (see Fig. 1), which is prescribed as
$F^{\dagger}({\boldsymbol{x}}^{\dagger},t^{\dagger})=0.$  (1) 
The velocity field ${\boldsymbol{u}}^{\dagger}={\boldsymbol{u}}^{\dagger}({\boldsymbol{x}}^{% \dagger},t^{\dagger})$ and vorticity ${\boldsymbol{\omega}}^{\dagger}\equiv\nabla^{\dagger}\times{\boldsymbol{u}}^{\dagger}$ satisfy the equations
$\displaystyle\frac{\partial{\boldsymbol{\omega}}^{\dagger}}{\partial{t}^{% \dagger}}+[{\boldsymbol{\omega}}^{\dagger},{\boldsymbol{u}}^{\dagger}]^{% \dagger}=0,\quad\nabla^{\dagger}\cdot{\boldsymbol{u}}^{\dagger}=0,$  
$\displaystyle\nabla^{\dagger}\equiv\left(\partial/\partial x_{1}^{\dagger},% \partial/\partial x_{2}^{\dagger},\partial/\partial x_{3}^{\dagger}\right),$  (2) 
where ${\boldsymbol{x}}^{\dagger}=(x_{1}^{\dagger},x_{2}^{\dagger},x_{3}^{\dagger})$ are Cartesian coordinates, ${t}^{\dagger}$ is time, daggers denote dimensional variables, and square brackets stand for the commutator of two vector fields $[{\boldsymbol{a}},{\boldsymbol{b}}]\equiv({\boldsymbol{b}}\cdot\nabla){% \boldsymbol{a}}({\boldsymbol{a}}\cdot\nabla){\boldsymbol{b}}$. The kinematic boundary condition at $\partial Q$ is
${\mathrm{d}}F^{\dagger}/{\mathrm{d}}t^{\dagger}=0\quad\text{at}\quad F^{% \dagger}({\boldsymbol{x}}^{\dagger},t^{\dagger})=0.$  (3) 
We consider oscillating flows that possess characteristic scales of velocity and length, together with two additional timescales:
$U^{\dagger},\quad L^{\dagger},\quad T_{{\rm fast}}^{\dagger},\quad T_{{\rm slow% }}^{\dagger}.$  (4) 
There are therefore two independent dimensionless parameters,
$T_{{\rm fast}}\equiv T_{{\rm fast}}^{\dagger}/T^{\dagger},\quad T_{{\rm slow}}% \equiv T_{{\rm slow}}^{\dagger}/T^{\dagger},\quad\text{where}\quad T^{\dagger}% \equiv L^{\dagger}/U^{\dagger},$  (5) 
which represent the dimensionless timescales. The scale $T_{{\rm fast}}$ characterises the given period of oscillations; hence the dimensional and dimensionless frequencies of oscillation are
$\sigma^{\dagger}\equiv 1/T_{{\rm fast}}^{\dagger},\quad\sigma\equiv T^{\dagger% }/T_{\text{fast}}^{\dagger}.$  (6) 
We choose the dimensionless independent variables as
$\displaystyle{\boldsymbol{x}}\equiv{\boldsymbol{x}}^{\dagger}/L^{\dagger},% \quad t\equiv t^{\dagger}/T^{\dagger}.$  (7) 
The dimensionless 'fast time' $\tau$ and 'slow time' $s$ are defined as:
$\tau\equiv t/T_{{\rm fast}}=\sigma t,\qquad s\equiv t/T_{{\rm slow}}\equiv St,% \quad\text{with}\quad S\equiv T^{\dagger}/T^{\dagger}_{{\rm slow}}.$  (8) 
We seek oscillatory solutions of Eq. (1) in the form
${\boldsymbol{u}}^{\dagger}=U^{\dagger}{\boldsymbol{u}}({\boldsymbol{x}},s,\tau);$  (9) 
the $\tau$dependence is $2\pi$periodic whereas, in general, the $s$dependence is not periodic. Transforming Eq. ((1)) to dimensionless variables and using the chain rule gives
$\left(\frac{\partial}{\partial\tau}+\frac{S}{\sigma}\frac{\partial}{\partial s% }\right){\boldsymbol{\omega}}+\frac{1}{\sigma}[{\boldsymbol{\omega}},{% \boldsymbol{u}}]=0.$  (10) 
The natural small parameter in our consideration is $1/\sigma$. The essence of the twotiming method is based on the assumption that the ratio $T_{\text{slow}}/T_{\text{fast}}=S/\sigma$ also represents a small parameter. As a result, Eq. (10) contains two independent small parameters, $\varepsilon_{1}$ and $\varepsilon_{2}$:
${\boldsymbol{\omega}}_{\tau}+\varepsilon_{1}{\boldsymbol{\omega}}_{s}+% \varepsilon_{2}[{\boldsymbol{\omega}},{\boldsymbol{u}}]=0;\quad\varepsilon_{1}% \equiv\frac{S}{\sigma}\ll 1,\quad\varepsilon_{2}\equiv\frac{1}{\sigma}\ll 1,$  (11) 
where the subscripts $s$ and $\tau$ denote partial derivatives.
Then, in the twotiming method, we make the standard auxiliary (but technically essential) assumption that the variables $s$ and $\tau$ are (temporarily) considered to be mutually independent. Its justification can be given a posteriori after solving (11), rewriting the solution in terms of the original variable $t$, and estimating the errors/residuals in the original equation (1), also expressed in terms of $t$ ([Yudovich2006]).
Let us temporarily forget about the definitions of $\varepsilon_{1}$ and $\varepsilon_{2}$ in (11) and treat them as abstract small parameters. In order to construct a rigorous asymptotic procedure with $(\varepsilon_{1},\varepsilon_{2})\to(0,0)$ we have to consider the various paths approaching the origin in the $(\varepsilon_{1},\varepsilon_{2})$plane. One may expect that there are infinitely many different asymptotic solutions to (11) corresponding to different paths (the usual sequence of the limits $\varepsilon_{1}\to 0$ and then $\varepsilon_{2}\to 0$, or with the order reversed, correspond to the 'broken' paths). However, for (11) (as well as for many other equations) one can find a few exceptional paths, which we shall call the distinguished limits. The notion of a distinguished limit is imprecisely defined (see, for example, [Nayfeh1973]; [Kevorkian and Cole1996]), varying between different books and papers. We suppose that a distinguished limit is given by a path that allows us to build a selfconsistent asymptotic procedure, leading to a finite/valid solution in any approximation. No systematic procedure of finding all possible distinguished paths is known, and so this may be regarded as still a problem of experimental mathematics.
We have considered in detail a number of different paths parametrized by
$\varepsilon_{1}=\delta^{k}\quad\text{and}\quad\varepsilon_{2}=\delta^{l}\quad% \text{with}\quad\delta\to 0,$  (12) 
where $k$ and $l$ are integers. From our search, we have found only two distinguished paths; these allow us to build the solutions
$\displaystyle\varepsilon_{1}=\delta,\ \varepsilon_{2}=\delta:$  $\displaystyle{\boldsymbol{\omega}}_{\tau}+\delta{\boldsymbol{\omega}}_{s}+% \delta[{\boldsymbol{\omega}},{\boldsymbol{u}}]=0,$  (13)  
$\displaystyle\varepsilon_{1}=\delta^{2},\ \varepsilon_{2}=\delta:$  $\displaystyle{\boldsymbol{\omega}}_{\tau}+\delta^{2}{\boldsymbol{\omega}}_{s}+% \delta[{\boldsymbol{\omega}},{\boldsymbol{u}}]=0.$  (14) 
The solutions may be expressed as the regular series
$({\boldsymbol{\omega}},{\boldsymbol{u}})=\sum_{k=0}^{\infty}\delta^{k}({% \boldsymbol{\omega}}_{k},{\boldsymbol{u}}_{k}),\quad k=0,1,2,\ldots.$  (15) 
The difference between the cases (13) and (14) appears in the main (zeroth order) approximation. For (13), the averaged 'standard' vortex dynamics takes place in the leading order approximation,
$\overline{{\boldsymbol{\omega}}}_{0}\neq 0,\qquad\overline{{\boldsymbol{u}}}_{% 0}\neq 0,$  (16) 
subsequent approximations producing various 'oscillatory' and 'mean' corrections. This is the case of Strong Vortex Dynamics (SVD). In contrast, for Eq. (14), the fluid motion is purely oscillatory in the main approximation,
$\overline{{\boldsymbol{\omega}}}_{0}\equiv 0,\quad\overline{{\boldsymbol{u}}}_% {0}\equiv 0.$  (17) 
Hence for the case (14) we consider only a relatively weak vorticity developing on the background wave motion. This leads to the CraikLeibovich equation and to Weak Vortex Dynamics (WVD). All other cases (12) that we have considered can be transformed either to one of these two main cases, or else they produce inconsistent/unsolvable systems of successive approximations, or else they lead to secular growth in $s$.
Equations (13) or (14) must be complemented by the boundary condition (3), with the same ordering of small parameters. This leads respectively to:
$\displaystyle F_{\tau}+\delta F_{s}+\delta{\boldsymbol{u}}\cdot\nabla F=0\quad% \text{at}\quad F({\boldsymbol{x}},s,t)=0,$  (18)  
$\displaystyle F_{\tau}+\delta^{2}F_{s}+\delta{\boldsymbol{u}}\cdot\nabla F=0% \quad\text{at}\quad F({\boldsymbol{x}},s,t)=0,$  (19) 
where the prescribed deformed oscillating boundary (1) , in its dimensionless form, is given by the exact expression
$F=\overline{F}_{0}({\boldsymbol{x}},s)+\delta\widetilde{F}_{1}({\boldsymbol{x}% },s,\tau)=0,$  (20) 
with given functions $\overline{F}_{0}$ and $\widetilde{F}_{1}$.
In order to make analytic progress we require a number of specific assumptions. Thus we assume that any dimensionless function $f({\boldsymbol{x}},s,\tau)$:

•
is of order one, $f={O}(1)$; and that all required ${\boldsymbol{x}}$, $s$, and $\tau$derivatives of $f$ are also ${O}(1)$;

•
is $2\pi$periodic in $\tau$, i.e. $f({\boldsymbol{x}},s,\tau)=f({\boldsymbol{x}},s,\tau+2\pi)$;

•
has an average given by
$\displaystyle\langle{f}\,\rangle\equiv\frac{1}{2\pi}\int_{\tau_{0}}^{\tau_{0}+% 2\pi}f({\boldsymbol{x}},s,\tau)\,{\mathrm{d}}\tau\equiv\overline{f}({% \boldsymbol{x}},s)\qquad\forall\ \tau_{0};$ 
•
can be split into averaged and purely oscillating parts, $f({\boldsymbol{x}},s,\tau)=\overline{f}({\boldsymbol{x}},s)+\widetilde{f}({% \boldsymbol{x}},s,\tau)$; the tildefunctions (or purely oscillating functions) are such that $\langle\widetilde{f}\,\rangle=0$ and the barfunctions are $\tau$independent. Furthermore, we introduce the tildeintegration which keeps the result in the tildeclass:
$\widetilde{f}^{\tau}\equiv\int_{0}^{\tau}\widetilde{f}({\boldsymbol{x}},s,% \sigma)\,{\mathrm{d}}\sigma\frac{1}{2\pi}\int_{0}^{2\pi}\left(\int_{0}^{\mu}% \widetilde{f}({\boldsymbol{x}},s,\sigma)\,{\mathrm{d}}\sigma\right)\,{\mathrm{% d}}\mu.$ (21)
3 Weak Vortex Dynamics (WVD)
In WVD we seek the solution of Eq. (14) ,
${\boldsymbol{\omega}}_{\tau}+\delta[{\boldsymbol{\omega}},{\boldsymbol{u}}]+% \delta^{2}{\boldsymbol{\omega}}_{s}=0,\quad\delta\to 0,$  (22) 
in the form of the regular series (15) . We restrict the class of possible solutions by imposing (17) .
The equations for successive approximations show that the zeroth order approximation of (22) is $\widetilde{{\boldsymbol{\omega}}}_{0\tau}=0$; its unique solution (within the tildeclass) is $\widetilde{{\boldsymbol{\omega}}}_{0}\equiv 0$. Together with (17) it shows that the full vorticity vanishes,
${\boldsymbol{\omega}}_{0}\equiv 0,$  (23) 
which means that the velocity field at leading order is purely oscillatory and potential. Then, similarly, the equation of the first order approximation of (22) yields $\widetilde{{\boldsymbol{\omega}}}_{1\tau}=0$. Its unique solution (within the tildeclass) is $\widetilde{{\boldsymbol{\omega}}}_{1}\equiv 0$, while the mean value $\overline{{\boldsymbol{\omega}}}_{1}$ remains undetermined. We write this symbolically as
$\widetilde{{\boldsymbol{\omega}}}_{1}\equiv 0,\quad\overline{{\boldsymbol{% \omega}}}_{1}=\ {\rm?}.$  (24) 
The second order approximation that takes into account both (23) and (24) is
$\widetilde{{\boldsymbol{\omega}}}_{2\tau}=[\overline{{\boldsymbol{\omega}}}_{% 1},\widetilde{{\boldsymbol{u}}}_{0}],$  (25) 
which after the use of tildeintegration (21) yields
$\widetilde{{\boldsymbol{\omega}}}_{2}=[\widetilde{{\boldsymbol{u}}}_{0}^{\tau}% ,\overline{{\boldsymbol{\omega}}}_{1}],\quad\overline{{\boldsymbol{\omega}}}_{% 2}=\ {\rm?}.$  (26) 
The third order approximation that takes into account both (23) and (24) is
$\widetilde{{\boldsymbol{\omega}}}_{3\tau}+\overline{{\boldsymbol{\omega}}}_{1s% }+[{\boldsymbol{\omega}}_{2},\widetilde{{\boldsymbol{u}}}_{0}]+[\overline{{% \boldsymbol{\omega}}}_{1},{\boldsymbol{u}}_{1}]=0.$  (27) 
Its average (barred) part is
$\overline{{\boldsymbol{\omega}}}_{1s}+[\overline{{\boldsymbol{\omega}}}_{1},% \overline{{\boldsymbol{u}}}_{1}]+\langle[\widetilde{{\boldsymbol{\omega}}}_{2}% ,\widetilde{{\boldsymbol{u}}}_{0}]\rangle=0,$  (28) 
which can be transformed with the use of (26) and the Jacobi identity to
$\displaystyle\overline{{\boldsymbol{\omega}}}_{1s}+[\overline{{\boldsymbol{% \omega}}}_{1},\overline{{\boldsymbol{u}}}_{1}+\overline{{\boldsymbol{V}}}_{0}]% =0,$  (29)  
$\displaystyle\overline{{\boldsymbol{V}}}_{0}\equiv\frac{1}{2}\langle[% \widetilde{{\boldsymbol{u}}}_{0},\widetilde{{\boldsymbol{u}}}_{0}^{\tau}]\rangle.$  (30) 
It can be seen that if $\widetilde{{\boldsymbol{u}}}_{0}$ is solenoidal then the drift velocity $\overline{{\boldsymbol{V}}}_{0}$ is also solenoidal, i.e. $\nabla\cdot\overline{{\boldsymbol{V}}}_{0}=0$.
After dropping subscripts and bars in $\overline{{\boldsymbol{u}}}_{1}$ and $\overline{{\boldsymbol{\omega}}}_{1}$, Eq. (29) can be used as the WVD model for the evolution of the averaged vorticity:
${\boldsymbol{\omega}}_{s}+[{\boldsymbol{\omega}},{\boldsymbol{w}}]=0,\quad% \text{where }\quad{\boldsymbol{w}}\equiv{\boldsymbol{u}}+\overline{{% \boldsymbol{V}}}_{0},$  (31) 
which shows that the averaged vorticity is frozen into the 'velocity $+$ drift'. This result is known as the CraikLeibovich equation (CLE) (see, for example, [Craik1985]). The derivation of the CLE here is much simpler technically than previous derivations, and minimises the number of assumptions needed (e.g. those on the flow geometry). We should emphasize that the drift velocity here is not considered to be small; it is of the same order of magnitude as the Eulerian averaged velocity. Equation ((31) may be integrated (in space) as
${\boldsymbol{u}}_{s}+({\boldsymbol{u}}\cdot\nabla){\boldsymbol{u}}+{% \boldsymbol{\omega}}\times\overline{{\boldsymbol{V}}}_{0}=\nabla p,\quad% \nabla\cdot{\boldsymbol{u}}=0,$  (32) 
where $p$ is a function of integration (a modified pressure) and the second equation follows from the continuity equation in ((1)${\Lsh}$).
Next, we should derive the averaged kinematic boundary condition for LABEL:boundsps6a)((20) . Following similar steps to those used when deriving (22)(31) leads to the averaged equation
$\overline{F}_{0s}+{\boldsymbol{w}}\cdot\nabla\overline{F}_{0}=0,\quad{% \boldsymbol{w}}\equiv{\boldsymbol{u}}+\overline{{\boldsymbol{V}}}_{0}.$  (33) 
When the averaged boundary does not depend on $s$, it is given by the equation
$F=\overline{F}_{0}({\boldsymbol{x}})+\delta\widetilde{F}_{1}({\boldsymbol{x}},% \tau)=0,$  (34) 
and (33) gives the is given for both Eqs. (35) and (40), the citation of equations has been changed in text. Kindly check and confirm the citation of Eqs. (35) and (40) throughout the article. averaged 'noleak' condition:
${\boldsymbol{w}}\cdot\overline{{\boldsymbol{n}}}_{0}=0\quad\text{at}\quad% \overline{F}_{0}({\boldsymbol{x}})=0,$  (35) 
where $\overline{{\boldsymbol{n}}}_{0}$ is the main approximation to the unit normal to $\partial Q$,
${\boldsymbol{n}}({\boldsymbol{x}},s,\tau)=\nabla F/\nabla F=\overline{{% \boldsymbol{n}}}_{0}({\boldsymbol{x}},s)+\varepsilon\widetilde{{\boldsymbol{n}% }}_{1}({\boldsymbol{x}},s,\tau)+\cdots.$  (36) 
Of course, for the averaged velocity the 'noleak' condition (35) corresponds to the presence of a leak:
${\boldsymbol{u}}\cdot\overline{{\boldsymbol{n}}}_{0}=\overline{{\boldsymbol{V% }}}_{0}\cdot\overline{{\boldsymbol{n}}}_{0}\quad\text{at}\quad\partial Q_{0}.$  (37) 
The effective boundary $\partial Q_{0}$ for this averaged flow is given by the equation $\overline{F}_{0}({\boldsymbol{x}},s)=0$; it means that the boundary conditions are prescribed not at the real boundary, but at its averaged position. Equations (31) , (32) and (33) or (37) form the closed model describing the averaged WVD flow.
The drift velocity $\overline{{\boldsymbol{V}}}_{0}$ is to be calculated from (30)$\Lsh$, where $\widetilde{{\boldsymbol{u}}}_{0}$ represents the solution of the previous approximation, $\widetilde{{\boldsymbol{u}}}_{0\tau}=\nabla\widetilde{p}_{0}$ and $\nabla\cdot\widetilde{{\boldsymbol{u}}}_{0}=0$, together with the boundary condition $\widetilde{F}_{1\tau}+\widetilde{{\boldsymbol{u}}}_{0}\cdot\nabla\overline{F}_% {0}$ at $\overline{F}_{0}=0$.
4 Alternative Scaling for the Slow Time
It is possible to find alternative scalings for the slow time scale $s$, while respecting the constraints given by the distinguished limits (13) , (14) . The slow time is defined in such a way that intervals of order one in $s$ correspond to changes of order one in the physical fields. In the SVD the mean velocity (16) is $O(1)$, and so we must have $s=t$. Physically, this means that in order to transport an admixture a dimensionless distance of order one, we need a dimensional time of order one. In the WVD the mean velocity (17) is $O(\delta)$. Then advection with $\overline{{\boldsymbol{u}}}_{0}=O(\delta)$ requires the slow timescale $s=\delta\,t$ (for $s=1$ the interval of the original 'physical' time is $1/\delta$).
The new formal small parameter $\delta$, introduced to describe the distinguished limit, can be related to the fast time scale $1/\sigma$ in a number of different ways. It is instructive to rewrite Eq. (22) as
$\sigma{\boldsymbol{\omega}}_{\tau}+\frac{1}{\sigma^{\alpha}}{\boldsymbol{% \omega}}_{s}+\sigma^{\beta}[{\boldsymbol{\omega}},{\boldsymbol{u}}]=0,\quad% \sigma\to\infty,$  (38) 
with constants $\alpha$ and $\beta$. In order to make (22) coincide with (38) we require
$\beta=(1\alpha)/2,\quad\delta=1/\sigma^{(\alpha+1)/2}.$  (39) 
Equations (38) , (39) can be interpreted in the following way: the related slow timescale is $s=t/\sigma^{\alpha}$ (where $\alpha>1$, which means that $s$ is a 'slow' variable in comparison with $\tau$) and the velocity is $\sigma^{\beta}{\boldsymbol{u}}$, not ${\boldsymbol{u}}$.
This transformation allows one to vary the slow timescale. Consider, for example, the case when [instead of (20) ] the boundary is prescribed as
$F=\overline{F}_{0}({\boldsymbol{x}},t)+\delta\widetilde{F}_{1}({\boldsymbol{x}% },t,\tau)=0,$  (40) 
which can appear in many practical applications. In this case the slow timescale $s=t$ is prescribed by the boundary condition; in fact we have $\alpha=0$ in (39) (and so for WVD we have the same time scales as for SVD) namely
$\tau=\sigma t,\qquad s=t.$  (41) 
However, if this is to be the correct scaling in the WVD then we must have $\beta=1/2$, so that the velocity at the boundary is $O(\sqrt{\sigma})$, not $O(1)$; also, from (39) , the small parameter of the decomposition should be chosen as $\delta=1/\sqrt{\sigma}$. Another interesting possibility corresponds to $\alpha=\beta=1/3$. In this case, $s=t/\sqrt[3]{\sigma}$, $\delta=\sigma^{2/3}$ and the velocity is $\sqrt[3]{\sigma}{\boldsymbol{u}}$. Although such an asymptotic scaling may look exotic, it would be required if the particular slow timescale $s=t/\sqrt[3]{\sigma}$ were prescribed by the boundary conditions. The original case (11) , which corresponds to an $O(1)$ velocity, corresponds to $\alpha=1$ and $\beta=0$. The general tendency is physically natural: to shorten the slow timescale (decreasing $\alpha$), one needs to increase the amplitude of the boundary oscillations (increasing $\beta$) (see Eq. (39) ).
5 Stokes Drift and Langmuir Vortices
In order to connect the above model equations (31) with classical areas of fluid dynamics, let us show that for a plane surface wave $\overline{{\boldsymbol{V}}}_{0}$, from (30)$\Lsh$, gives the classical Stokes drift and that this then leads to an understanding of the nature of Langmuir vortices.
The dimensional solution for a plane potential travelling gravity wave is
$\displaystyle{\boldsymbol{u}}^{\dagger}_{0}=U^{\dagger}\widetilde{{\boldsymbol% {u}}}_{0},\quad\widetilde{{\boldsymbol{u}}}_{0}=\exp(k^{\dagger}z^{\dagger})% \left(\begin{array}[]{c}\cos(k^{\dagger}x^{\dagger}\tau)\\ \sin(k^{\dagger}x^{\dagger}\tau)\end{array}\right),$ 
where $U^{\dagger}={k^{\dagger}g^{\dagger}a^{\dagger}}/{\sigma^{\dagger}}$ with $\sigma^{\dagger}$, $a^{\dagger}$, and $g^{\dagger}$ the dimensional frequency, spatial wave amplitude, and gravity (see [Stokes1847]; [Lamb1932]; [Debnath1994]). Then Eq. (30)$\Lsh$ yields
$\displaystyle\widetilde{{\boldsymbol{u}}}_{0}=e^{z}\left(\begin{array}[]{c}% \cos(x\tau)\\ \sin(x\tau)\end{array}\right),\quad\overline{{{\boldsymbol{V}}}}_{0}=e^{2z}% \left(\begin{array}[]{c}1\\ 0\end{array}\right).$ 
The dimensional version is
$\overline{{{\boldsymbol{V}}}}_{0}^{\dagger}=\frac{U^{2}k^{\dagger}}{\sigma^{% \dagger}}e^{2k^{\dagger}z^{\dagger}}\left(\begin{array}[]{c}1\\ 0\end{array}\right),$  (42) 
which agrees with the classical expression for the drift velocity ([Lamb1932]; [Debnath1994]; [Batchelor1967]). To obtain (42) one should take into account that the transformation to the physical formula for drift includes a move from the slow time $s=t/\sigma$ to the physical time $t$.
The structure of the CLE and WVD can be seen as a relatively passive alteration of the original Euler equations, since we still have frozenin vorticity dynamics. However, the additional terms (which contain the drift velocity) make a qualitative change to the properties of the solutions. One example of such a new property is related to the Langmuir circulations (see Fig. 2). In order to illustrate such a qualitative change, let us consider the class of translationally invariant averaged flows. Let the zeroth approximation (23) take the form of a plane potential travelling gravity wave with the drift velocity (30). Let Cartesian coordinates $(x,y,z)$ be such that $\overline{{\boldsymbol{V}}}_{0}=(U,0,0)$, $U=e^{2z}$, $\overline{{\boldsymbol{u}}}_{1}=(u,v,w)$, where all components are $x$independent (translationallyinvariant).
Then the component form of ((32) is
$\displaystyle u_{s}+vu_{y}+wu_{z}=0,$  
$\displaystyle v_{s}+uv_{y}+wv_{z}Uu_{y}=\overline{p}_{y},$  
$\displaystyle w_{s}+vw_{y}+ww_{z}Uu_{z}=\overline{p}_{z},$  
$\displaystyle v_{y}+w_{z}=0,$ 
which can be rewritten (see [Vladimirov1985a]; [Vladimirov1985b]) as
$\displaystyle v_{s}+vv_{y}+wv_{z}=\overline{P}_{y}\rho\Phi_{y},$  (43)  
$\displaystyle w_{s}+vw_{y}+ww_{z}=\overline{P}_{z}\rho\Phi_{z},$  
$\displaystyle v_{y}+w_{z}=0,$  
$\displaystyle\rho_{s}+u\rho_{x}+v\rho_{y}=0,$ 
where $\rho\equiv u$, $\Phi\equiv U=e^{2z}$ and $\overline{P}$ is a new modified pressure. One can see that Eq. (43) are mathematically equivalent to the system of equations for an incompressible stratified fluid in the Boussinesq approximation. The effective 'gravity field' ${\boldsymbol{g}}=\nabla\Phi=(0,0,2e^{2z})$ is nonhomogeneous, which makes the analogy with the 'standard' stratified fluid incomplete. Nevertheless, longitudinal vortices should naturally appear in (43) as a RayleighTaylor type instability of an inversely stratified equilibrium corresponding to $(u,v,w)=(u(z),0,0)$ with any increasing function $u(z)\equiv\rho(z)$ (see Fig. 2).
6 Energy Variational Principle and Arnold Stability
The energy integral for the averaged motion can be written as:
$E=E(s)=\frac{1}{2}\int_{Q}({\boldsymbol{u}}+\overline{{\boldsymbol{V}}}_{0})^{% 2}{\mathrm{d}}{\boldsymbol{x}}={\mathrm{const}}.$  (44) 
One can show, by virtue of (32) , that its $s$derivative can be written as
$\frac{{\mathrm{d}}E}{{\mathrm{d}}s}=\int_{Q}\left(p+\frac{{\boldsymbol{u}}^{2% }}{2}\right)({\boldsymbol{u}}+\overline{{\boldsymbol{V}}}_{0})\cdot\overline{{% \boldsymbol{n}}}_{0}\,{\mathrm{d}}{\boldsymbol{x}}=0,$  (45) 
which is zero owing to (35) .
According to (32) , vorticity is frozen into ${\boldsymbol{u}}+\overline{{\boldsymbol{V}}}_{0}$. We may then use a slightly modified Arnold isovorticity condition ([Arnold and Khesin1999]) in its differential form,
$\displaystyle{\boldsymbol{u}}_{\theta}={\boldsymbol{f}}\times{\boldsymbol{% \omega}}+\nabla\alpha,\quad{\mathrm{div}\,}{\boldsymbol{u}}=0,\quad{\mathrm{% div}\,}{\boldsymbol{f}}=0\quad\text{in}\quad Q_{0};$  (46)  
$\displaystyle({\boldsymbol{u}}+\overline{{\boldsymbol{V}}}_{0})\cdot\overline{% {\boldsymbol{n}}}_{0}=0,\quad{\boldsymbol{f}}\cdot\overline{{\boldsymbol{n}}}_% {0}=0\quad\text{at}\quad\partial Q_{0};$ 
where ${\boldsymbol{u}}({\boldsymbol{x}},\theta)$ is the unknown function, ${\boldsymbol{f}}={\boldsymbol{f}}({\boldsymbol{x}},\theta)$ is an arbitrary given solenoidal function, $\theta$ is a scalar parameter along an isovortical trajectory, and subscript $\theta$ denotes a partial derivative. The function $\alpha({\boldsymbol{x}},\theta)$ is to be determined from the condition $\nabla\cdot{\boldsymbol{u}}=0$. The initial data at $\theta=0$ for ${\boldsymbol{u}}({\boldsymbol{x}},\theta)$ in (46) corresponds to a steady flow
${\boldsymbol{u}}({\boldsymbol{x}},0)={\boldsymbol{U}}({\boldsymbol{x}}),\quad{% \boldsymbol{\omega}}({\boldsymbol{x}},0)={\boldsymbol{\Omega}}({\boldsymbol{x}% }),$  (47) 
where ${\boldsymbol{U}}({\boldsymbol{x}})$ and ${\boldsymbol{\Omega}}({\boldsymbol{x}})$ represent the steady solutions ($\partial/\partial s=0$) of (31) and (32) with the noleak boundary conditions (35) .
Differentiation of $E$ with respect to $\theta$ produces the first variation
$E_{\theta}\Big_{\theta=0}=\int_{Q_{0}}{\boldsymbol{f}}\left({\boldsymbol{% \Omega}}\times{\boldsymbol{W}}\right){\mathrm{d}}{\boldsymbol{x}}=0,\quad{% \boldsymbol{W}}\equiv{\boldsymbol{U}}+\overline{{\boldsymbol{V}}}_{0},$  (48) 
which vanishes for any solenoidal function ${\boldsymbol{f}}$ by virtue of the equations of motion and the boundary conditions for steady flow. This equality gives us the variational principle: any steady flow represents a stationary point on the isovortical sheet. The only difference from Arnold's classical result is the boundary conditions in the definition of the isovorticity sheet (46) .
The calculation of the second variation yields:
$E_{\theta\theta}\Big_{\theta=0}=\int_{Q_{0}}\left({\boldsymbol{u}}_{\theta}% ^{2}+({\boldsymbol{W}}\times{\boldsymbol{f}})\cdot{\boldsymbol{\omega}}_{% \theta}\right)_{\theta=0}{\mathrm{d}}{\boldsymbol{x}},$  (49) 
where ${\boldsymbol{W}}\equiv{\boldsymbol{U}}+\overline{{\boldsymbol{V}}}_{0}$. This is analogous to Arnold's result; expression (49) shows that the stationary point of the energy functional in the threedimensional case always represents a saddle point.
However, the stability conditions can be obtained for steady plane flows in the case when a stream function $\Psi(x_{1},x_{2})$ for the combined velocity ${\boldsymbol{W}}(x_{1},x_{2})$ can be introduced as $W_{1}=\partial\Psi/\partial x_{2}$, $W_{2}=\partial\Psi/\partial x_{1}$. For the plane flow the second variation (49) , combined with the standard Casimir integral containing an arbitrary function of vorticity ([Arnold and Khesin1999]) takes the form
$E_{\theta\theta}\Big_{\theta=0}=\int_{Q_{0}}\left({\boldsymbol{u}}_{\theta}% ^{2}\frac{d\Psi}{d\Omega}\omega_{\theta}^{2}\right)_{\theta=0}\,{\mathrm{d}}x% _{1}{\mathrm{d}}x_{2},$  (50) 
where $\omega({\boldsymbol{x}},\theta)$ and $\Omega({\boldsymbol{x}})$ are the $x_{3}$components of the full vorticity and the steady vorticity at $\theta=0$, and the functional dependence $\Psi=\Psi(\Omega)$ characterises the plane steady flow under consideration. Then, similarly to the Arnold cases, the inequalities with two positive (or two negative) constants $C^{}$, $C^{+}$ satisfying
$C^{}<\frac{d\Psi}{d\Omega}<C^{+}$  (51) 
give both sufficient linear and nonlinear stability conditions for the positively (and negatively) defined energycasimir functional. One should take into account that these stability conditions determine stability with respect to arbitrary perturbations, not only isovortical ones. However, this part of the analysis is similar to Arnold's wellknown results of 1966 (see [Arnold and Khesin1999]) and is not presented here.
Finally, we note that for plane WVD flows we have derived sufficient stability conditions that differ from the classical ones only by replacing the streamfunction for the velocity field ${\boldsymbol{U}}$ by the streamfunction for the combined velocity ${\boldsymbol{W}}$. Therefore an important conclusion for the stability of the plane WVD flows is that virtually any plane steady flow can be made stable by the choice of the 'proper' field of drift velocity.
7 Discussion
Our main achievement in this paper is a significant simplification of the derivation of the CraikLeibovich equation (CLE). Its most known derivation ( [Craik1985] ) is performed with the use of the Generalized Lagrangian Mean theory (GLM) (see [Andrews and McIntyre1978] ; [Craik1985] ; [Buhler2009] ) and further theoretical studies are often performed in GLM terms ( [Holm1996] ). In contrast, here we introduce the CLE in its natural simplicity and generality and use only the standard Eulerian description, making our derivation much more accessible.
Examples of oscillating flows are not restricted by a deformed domain as in Fig. 1 ; they can be extended to many oscillating flows that appear in practical applications, such asoscillating or rotating rigid bodies, moving pistons and acoustics ( [Vladimirov and Ilin2013] ); some of these cases are illustrated in Fig. 3 .
The two small parameters that we have used represent two ratios of three timescales. It should be noted that in the derivation of CLE we have not used the amplitude as a small parameter. The main field of velocity oscillations $\widetilde{{\boldsymbol{u}}}_{0}$ is of dominant order and so is not small. Nevertheless, the oscillatory spatial amplitude of material particles and the related spatial amplitude of the deformation of the boundary $\partial Q(t)$ are both small.
It is instructive to derive the averaged equations of the second approximation (while the CLE (29) ) appears in the first approximation). This appears as the linearized equation at the first approximation and contains an additional 'force term' depending on the previous approximations. It means that some additional instabilities are possible, beyond classical instabilities of the linearized problem. This raises interesting questions about the meaning of linearization and its nonuniqueness.
It is remarkable that the same CLE (29) describes the WVD in the case of acoustics (see [Vladimirov and Ilin2013]), when $\widetilde{{\boldsymbol{u}}}_{0}$ represents a given acoustic wave that satisfies the wave equations and cannot be solenoidal. An important qualitative addition of the 'acoustical CLE' is that the drift velocity can be an arbitrary solenoidal function. It gives greater general significance to studies of the CLE, since an arbitrary solenoidal function $\overline{{\boldsymbol{V}}}_{0}$, Eq. (30) now has a practical meaning.
A viscous term can be straightforwardly added to the right hand side of the CLE (see [Craik1985] ). Our derivation (22)  (30) $\Lsh$ shows that in order to accommodate such an addition the dimensionless viscosity (or the inverse Reynolds number) should be of order $\delta^{3}$.
The same analysis as above is valid for stratified fluids in the Boussinesq approximation where the generalization of the CLE is straightforward (see [Craik1985] ). At the same time, the analogy with stratification (43) discussed earlier leads to Richardsontype stability criteria even in the case of the CLE for homogeneous fluid.
The CLE is Hamiltonian as is immediately clear from its form. This question was considered by [Holm1996] for the CLE in the GLM form introduced by [Andrews and McIntyre1978], [Craik1985], [Buhler2009], which is somewhat different from ours. However, the investigation of the Hamiltonian structure of our equations is beyond the scope of this paper. Similarly, it would be of interest to study the possibility of a finitetime vorticity singularity for the CLE.
A generalization of the CLE and WVD has been obtained for MHD by [Vladimirov2013]. In this case the drift velocity appears in both the equation for the advection of vorticity and the equation for the advection of magnetic field. Similar results restricted to kinematic MHD have also been obtained in [Vladimirov2010] and [Herreman and Lesaffre2011]. There remains the challenge of developing a selfconsistent theory of the full MHD equations ([Moffatt1978]), which may be viewed as a complementary approach to that of [Courvoisier et al.2010].
Acknowledgements
We thank Profs A.D.D. Craik, K.I. Ilin and H.K. Moffatt for helpful discussions.References
 [Andrews and McIntyre1978] Andrews, D., McIntyre, M.E.: An exact theory of nonlinear waves on a Lagrangianmean flow. J. Fluid Mech. 89, 609646 (1978)
 [Arnold and Khesin1999] Arnold, V.I. and Khesin, B.A.: Topological Methods in Hydrodynamics. Springer, Berlin (1999)
 [Batchelor1967] Batchelor, G.K.: An Introduction to Fluid Dynamics. CUP, Cambridge (1967)
 [Buhler2009] Buhler, O.: Waves and Mean Flows. CUP, Cambridge (2009)
 [Craik1985] Craik, A.D.D.: Wave Interactions and Fluid Flows. CUP, Cambridge (1985)
 [Courvoisier et al.2010] Courvoisier, A., Hughes, D.W., Proctor, M.R.E.: Selfconsistent meanfield magnetohydrodynamics. Proc. R. Soc. A 466, 583601 (2010)
 [Debnath1994] Debnath, L.: Nonlinear Water Waves. Academic Press, Boston (1994)
 [Herreman and Lesaffre2011] Herreman, W., Lesaffre, P.: Stokes drift dynamos. J. Fluid Mech. 679, 3257 (2011)
 [Holm1996] Holm, D.D.: The ideal CraikLeibovich equation. Phys. D 98, 415441 (1996)
 [Kevorkian and Cole1996] Kevorkian, J., Cole, J.D.: Multiple Scale and Singular Perturbation Method. Applied Mathematical Sciences, vol. 114. Springer, New York (1996)
 [Lamb1932] Lamb, H.: Hydrodynamics, 6th edn. CUP, Cambridge (1932)
 [Leibovich1983] Leibovich, S.: The forms and dynamics of Langmuir circulations. Annu. Rev. Fluid Mech. 15, 391427 (1983)
 [Moffatt1978] Moffatt, H.K.: Magnetic Field Generation in Electrically Conducting Fluids. CUP, Cambridge (1978)
 [Nayfeh1973] Nayfeh, A.H.: Perturbation Methods. Wiley, NY (1973)
 [Stokes1847] Stokes, G.G.: On the theory of oscillatory waves. Trans. Camb. Philos. Soc. 8, 441455 (1847). (Reprinted in Math. Phys. Papers, 1, 197219)
 [Thorpe2004] Thorpe, S.A.: Langmuir circulation. Annu. Rev. Fluid Mech. 36, 5579 (2004)
 [Vladimirov and Ilin2013] Vladimirov, V.A., Ilin, K.I.: An asymptotic model in acoustics: acoustic drift equations. J. Acoust. Soc. Am. 134(5), 34193424 (2013)
 [Vladimirov2013] Vladimirov, V.A.: MHD drift equation: from Langmuir circulations to MHD dynamo? J. Fluid Mech. 698, 5161 (2013)
 [Vladimirov2005] Vladimirov, V.A.: Vibrodynamics of pendulum and submerged solid. J. Math. Fluid Mech. 7, S397412 (2005)
 [Vladimirov2008] Vladimirov, V.A.: Viscous flows in a halfspace caused by tangential vibrations on its boundary. Studies Appl. Math. 121(4), 337367 (2008)
 [Vladimirov2010] Vladimirov, V.A.: Admixture and drift in oscillating fluid flows. (physics, fludyn) (2010)
 [Vladimirov1985a] Vladimirov, V.A.: An example of equivalence between density stratification and rotation. Sov. Phys. Dokl. 30(9), 748750 (1985). (translated from Russian)
 [Vladimirov1985b] Vladimirov, V.A.: Analogy between the effects of density stratification and rotation. J. Appl. Mech. Tech. Phys. 26(3), 353362 (1985). (translated from Russian)
 [Yudovich2006] Yudovich, V.I.: Vibrodynamics and vibrogeometry of mechanical systems with constrains. Uspehi Mekhaniki 4(3), 26158 (2006). (in Russian)