Streamwise-travelling viscous waves in channel flows

The unsteady viscous flow induced by streamwise-travelling waves of spanwise wall velocity in an incompressible laminar channel flow is investigated. Wall waves belonging to this category have found important practical applications, such as microfluidic flow manipulation via electro-osmosis and surface acoustic forcing and reduction of wall friction in turbulent wall-bounded flows. An analytical solution composed of the classical streamwise Poiseuille flow and a spanwise velocity profile described by the parabolic cylinder function is found. The solution depends on the bulk Reynolds number R, the scaled streamwise wavelength \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document}λ, and the scaled wave phase speed U. Numerical solutions are discussed for various combinations of these parameters. The flow is studied by the boundary-layer theory, thereby revealing the dominant physical balances and quantifying the thickness of the near-wall spanwise flow. The Wentzel–Kramers–Brillouin–Jeffreys (WKBJ) theory is also employed to obtain an analytical solution, which is valid across the whole channel. For positive wave speeds which are smaller than or equal to the maximum streamwise velocity, a turning-point behaviour emerges through the WKBJ analysis. Between the wall and the turning point, the wall-normal viscous effects are balanced solely by the convection driven by the wall forcing, while between the turning point and the centreline, the Poiseuille convection balances the wall-normal diffusion. At the turning point, the Poiseuille convection and the convection from the wall forcing cancel each other out, which leads to a constant viscous stress and to the break down of the WKBJ solution. This flow regime is analysed through a WKBJ composite expansion and the Langer method. The Langer solution is simpler and more accurate than the WKBJ composite solution, while the latter quantifies the thickness of the turning-point region. We also discuss how these waves can be generated via surface acoustic forcing and electro-osmosis and propose their use as microfluidic flow mixing devices. For the electro-osmosis case, the Helmholtz–Smoluchowski velocity at the edge of the Debye–Hückel layer, which drives the bulk electrically neutral flow, is obtained by matched asymptotic expansion.


Introduction
In this paper, we study the unsteady laminar flow generated in a Poiseuille flow channel by the following wall waves of sinusoidal spanwise velocity travelling along the streamwise direction: where w * w is the spanwise wall velocity, x * indicates the streamwise direction, t * denotes time, and A * denotes the oscillation amplitude. Two parameters define the wall motion: the streamwise wavelength λ * and the phase speed U * . The third parameter defining the physical system is the Reynolds number R based on the bulk velocity of the streamwise Poiseuille flow U * b and the half channel height h * . Our study is theoretical and numerical and aims at a complete characterization of the viscous flow in the parameter space. In the following, a variety of flow configurations where these waves may play an important role are discussed, including flow mixing in microfluidic systems and turbulent drag reduction.

Travelling waves in microfluidic systems
Small-scale oscillating flows often feature in microfluidic and micro-electromechanical systems. A benefit of the oscillations is the promotion of mixing in the flow, which, as a result of the small length scales and the small velocities involved, is essentially laminar. At such Reynolds numbers, R ≈ 0.1-10, oscillating flows such as that induced by the wall motion (1) can be enforced by surface acoustic waves or by electro-osmosis waves.
Surface acoustic waves (SAWs) have been utilized extensively for a wide range of microfluidic applications, such as micromixing, micropumping, drop transport, cell handling, and microejectors [1][2][3][4][5], although shear-horizontal waves have never been employed in microfluidic flow mixing. SAWs are created by piezoelectric transduction within a thin solid substrate below a fluid, so that electric power causes the mechanical deformation of the substrate, which, in turn, leads to the motion of the fluid.
Wall-normal Rayleigh acoustic waves have been used for mixing of microfluidic flows [6][7][8][9]. However, they generate compression waves in a liquid and suffer from energy dissipation (leaky waves) [6,10]. When instead in-plane motion occurs, thanks to the mismatch of the sound speeds and densities of the substrate material and the fluid, the acoustic propagation is confined within the substrate, while the fluid flow is incompressible. This is a relevant simplification for the analysis of SAWs because the incompressible Navier-Stokes equations with an imposed slip velocity describe the dynamics (i.e. an analytical solution for streamwise standing waves is found on p. 133 of the book by Bruus [11]). Tan et al. [4], in their Fig. 2c, show interfacial standing shear waves which are in-plane and sinusoidal as an interesting variant of SAWs. Although no bulk flow is present in the system studied by Neumann et al. [12], this example shows that it is possible to generate shear-horizontal acoustic waves in a thin solid substrate to affect an overlying liquid layer.
Shear-horizontal surface waves, also called Love waves when a layer of lower acoustic velocity is used for increased sensitivity, have also been studied extensively as efficient biosensors and chemical sensors for flowing solutions because of their low dissipation when compared with wall-normal Rayleigh waves [13,14]. However, these studies are mainly experimental and Lange et al. [13] indeed remark that improved design of these biosensors can be achieved by studying the fluid dynamics generated by the interaction of the spanwise waves and the overriding streamwise flow. Figure 1 depicts a schematic of a biosensor based on shear-horizontal travelling waves, an excellent technological application of the waves studied herein. Shear-horizontal SAWs waves have also been employed by Neumann et al. [12] to manipulate proteins attached to supported lipid bilayers. Love waves have also been used more recently as non-intrusive rheometers. In particular, the interaction between the SAWs and the liquid has been studied to extract the relationship between the wave attenuation and the viscosity [15,16].
SAWs are usually characterized by maximum surface velocities of 1 m/s and frequencies in the broad range of 100 kHz-100 MHz. The characteristic wavelength can be up to 100 µm, i.e. comparable with the microfluidic Fig. 1 Schematic of a shear-horizontal SAW biosensor, adapted from Lange et al. [13]. The large black arrows indicate the liquid sample flowing over the piezoelectric crystal substrate, the small black arrow denotes the electric signal exchanged between the two interdigital transducers, and the blue sinusoidal line shows the travelling surface acoustic wave channel height. Therefore, for Reynolds numbers of the order of unity, realistic ratios between the streamwise wavelength and the channel height can be in the range 0.1-1. The ratio between the wall wave speed and the bulk fluid velocity can vary greatly, namely from null to the order of 10 4 . The latter scenario is more common as it corresponds to a single acoustic wave, travelling at the sound speed of about 1500 m/s below a liquid flowing at about 100 µm/s. Nonetheless, standing acoustic waves can be generated by the interference of travelling waves [8].
In addition to mechanical wall motions, travelling wave on the walls of a microchannel of the form (1) can also be engendered via electro-osmosis [17][18][19]. Surface electrodes driven by AC currents below a fluid can generate a uniform plug flow within a very thin charged Debye-Hückel layer, which drags the overlying uncharged fluid by shear stresses. As remarked by Ajdari [18], since these layers are usually much thinner than the radii of curvature of the surface and the channel height, it can be assumed that the uncharged fluid is simply affected by an imposed effective slip velocity, which is linearly related to the electric field, and that the bulk flow is described by the incompressible Navier-Stokes equations. A concise explanation of the physics of wall-based electro-osmosis in microchannels is found in Chang and Yang [20]. These shear motions have also been proposed as an electro-osmotic pumping device to drive fluid along a channel [21,22]. Micromixing has also been successfully achieved through electro-osmotic wall forcing [23][24][25][26][27].
In line with these microfluidic mixing applications, we complement our theoretical/numerical results with ideas on the laboratory realization of the waves engendered by (1) through surface acoustic forcing by piezoelectric crystals and through electro-osmosis actuation for microfluidic flow mixing (refer to Sect. 8).

Travelling waves for turbulent drag reduction
The wall wave motion given by (1) has also been studied beneath wall-bounded turbulent flows, first via direct numerical simulations in a turbulent channel flow by Quadrio et al. [31] and Quadrio and Ricco [32], and experimentally in a pipe with rotating sections by Auteri et al. [28] and in a wind-tunnel flow over a deformable Kagome lattice surface by Bird et al. [29,30]. Drag reduction or drag increase occurs depending on the forcing parameters λ and U . For the pipe with rotating sections [28] and the wind-tunnel flow over the Kagome surface [29,30], the bulk Reynolds number is obviously much larger than unity and λ is either comparable or a few times larger than the pipe radius in the case of Auteri et al. [28] or the boundary-layer thickness in the case of Bird et al. [29,30].
It is obviously very different to investigate the flow engendered by the waves given by (1) in the laminar regime or in the fully developed turbulent regime. Nevertheless, there is ample evidence that the laminar profiles generated by spanwise wall motion are very useful to study various aspects of the corresponding fully turbulent flow. Choi et al. [33] and Quadrio and Ricco [32] have indeed verified that the unsteady space-averaged spanwise profile may closely match the corresponding laminar solution. The good agreement occurs when the wall forcing acts on a time scale which is much shorter than the life time of the near-wall turbulent structures. Under these conditions, the drag reduction scales with the thickness of the spanwise boundary layer, which is computed through the laminar solution. Furthermore, the near-wall laminar solutions have been instrumental for the accurate computation of the power spent for moving the wall against the viscous flow resistance, the optimal layer thickness which leads to maximum drag reduction, or the smallest period of wall forcing which guarantees drag reduction [32]. Choi et al. [33] and Ricco et al. [34] have also utilized the laminar Stokes layer solution to define a scaling parameter for drag reduction prediction and Choi [35] has taken advantage of the spanwise laminar flow behaviour to interpret the changes of the near-wall turbulent structures.
1.3 Objectives and structure of the paper Motivated by the possibility of microfluidic flow manipulation offered by shear-horizontal waves, by their extensive use as bio-and chemical sensors, and by the importance of the laminar solutions for the study of turbulent drag reduction by spanwise forcing, a complete study on the laminar spanwise flow engendered by the wall motion given by (1) is presented herein. The investigation is based on numerical calculations and on asymptotic analysis. The spanwise momentum equation is first simplified to a second-order ordinary differential equation and solved numerically by a second-order finite-difference scheme. Its solution is also expressed analytically through the parabolic cylinder function (hereinafter referred to as PCF).
The Reynolds number R, the wave speed U , and the wavelength λ are treated as asymptotic parameters, thus deriving asymptotic analytical solutions utilizing the boundary-layer and the Wentzel-Kramers-Brillouin-Jeffreys (WKBJ) theories. Employing the boundary-layer scaling, the thicknesses of the near-wall viscous layers are quantified, while the WKBJ solution gives the correct flow structure across the whole channel. When the streamwise diffusion is negligible, the convection term due to the wave motion may balance the streamwise Poiseuille flow convection term in a specified range of wave phase speeds. For these cases, alternative solutions found using WKBJ turning-point composite expansions and the method of Langer [36,37] are derived. All the asymptotic solutions show excellent agreement with the numerical solutions and offer the further advantage over the numerical approach that the dominant physical balances are revealed. The final aim of our work is to discuss how these waves can be generated in a laboratory via electro-osmosis and surface acoustic forcing, and to propose them as microfluidic flow mixers.
The scaling and the simplification of the spanwise momentum equation are discussed in Sect. 2, while Sect. 3 presents the analytical solution in terms of the PCF. The numerical results are shown in Sect. 4, the boundary-layer theory results are discussed in Sect. 5, while Sect. 6 presents the WKBJ results and Sect. 7 contains the Langer solution. A discussion on future applications of the travelling waves for flow mixing is contained in Sect. 8 and a summary is found in Sect. 9.

Governing equation
Laminar flow confined between two infinite parallel plates at a distance 2h * is considered. The superscript * hereinafter denotes a dimensional quantity. The streamwise, wall-normal, and spanwise directions are indicated by x * , y * , and z * , respectively, while t * denotes time. The walls move along the spanwise direction with velocity given by (1), where without loss of generality, the wavelength λ * > 0 because the flow is invariant to a change of λ * to −λ * . The flow is governed by the incompressible Navier-Stokes equations: where u * = {u * , v * , w * } is the velocity vector with components along x * , y * , and z * , p * is the pressure, ρ * is the density of the fluid, ν * is the kinematic viscosity of the fluid, and ∇ = {∂/∂ x * , ∂/∂y * , ∂/∂z * }. At the walls y * = 0 and y * = 2h * , the no-slip and no-penetration boundary conditions are u * = v * = 0, and w * = A * R e 2π i(x * −U * t * )/λ * , where R indicates the real part. The phase speed U * can be positive (forward-travelling wave), null (standing wave), or negative (backward-travelling wave). Figure 2 shows the flow domain for forward-travelling waves.
As the boundary conditions (3) depend only on x * and t * and a pressure gradient is present only along x * , all terms containing z * derivatives vanish. The continuity equation (2a) and the x-momentum equation (2b) are thus independent of w * . The streamwise flow is the classical Poiseuille flow, while w * satisfies the simplified zmomentum equation The x * coordinate is now scaled by the wall streamwise wavelength λ * , while the y * coordinate is scaled by h * to enable the definition of the dimensionless coordinates x = x * /λ * = O(1) and y = y h * 0 u * (y * ) dy * , while the spanwise velocity w * is scaled by the wave amplitude A * , i.e. w = w * /A * = O(1). The time is non-dimensionalized by the period of the wall motion, i.e. t = t * U * /λ * = O (1). In terms of these non-dimensional quantities Eq. (4) becomes Spanwise viscous flow Eq. (5) simplifies to subject to A prime indicates differentiation with respect to y and use has been made of the Poiseuille solution, u = 3y(1− y/2). Three parameters appear in Eq. (7): the Reynolds number R = U * b h * /ν * , the scaled phase speed U = U * /U * b , and the scaled wavelength λ = λ * / h * .

Analytical solution in terms of the parabolic cylinder function
By introducing the wall-normal coordinateŷ = φ(1 − y), where φ = √ 2[3π R/(λi)] 1/4 , and the spanwise velocity W ŷ = W(y), Eq. (7) simplifies to where a = √ iπ/(3λR)[iR(U − 3/2) − 2π/λ] − 1/2. The boundary conditions (8) become W (±φ) = 1. Equation (9) is in the form of Eq. (3.5.11) on p. 96 in Bender and Orszag [38] and therefore the solution of Eq. (5) can be written in terms of two linearly independent PCFs: While this expression provides an analytical solution for w, its practical value for determining the actual spanwise flow velocity is somewhat limited. This is because, to the best of our knowledge, there does not exist yet a robust numerical code which solves for the PCF in the entire complex plane. Temme [39] offers a complete list of journal articles on the numerical algorithms, outlining the restrictions on the complex argument. He states that "...constructing reliable software for all possible combinations of the complex parameters is a challenging problem".
Mathematica solves the complex-to-complex PCF, but it produces erroneous results for spanwise viscous layers extending to the whole channel width. As a full computational implementation of a numerical code for complexto-complex PCF lies outside of the scope of the present work, numerical and asymptotic methods will now be employed to study the flow.

Numerical results
Flow profiles are first obtained by solving Eq. (7) numerically through a second-order finite-difference scheme with uniform mesh size along y [40]. We have chosen an implicit scheme to avoid stiffness-related problems. As the flow is symmetric with respect to the centreline, only half of the domain along y needs to be considered, and the boundary conditions become W (0) = 1, W (1) = 0. In the following figures, the numerical profiles are shown for eight ξ values, equally spaced across one oscillation period. Increasing the phase speed U renders the spanwise viscous region thinner, as shown in Fig. 3 for R = 1, λ = 1, U = 10 (left), and U = 100 (right). This behaviour is the same as for λ/R 1, where the boundary-layer thickness becomes thinner as the frequency increases (refer to contour plot in Quadrio and Ricco [32,Fig. 8]). In the cases shown in Fig. 3, all four terms in Eq. (4) play a role, as the unsteadiness and the transport due to the Poiseuille flow are balanced by viscous diffusion along x and y. These profiles are very similar to the classical Stokes problem profiles. Note that, for λ = R = 1 and U = 10 (and smaller), a boundary layer does not exist as viscous effects diffuse across the whole wall-normal extent of the channel.
Flows for R = 1 with U = 0 and wavelengths λ = 0.01, 0.1, 1 are shown in Fig. 4(left). A boundary layer emerges as λ decreases, and the flow is self-similar as the wall-normal coordinate is rescaled by the streamwise wavelength,ȳ = y/λ = O(1), in the limit λ 1. The rescaled profiles show excellent agreement, especially near the wall when the wall velocity is large. By introducingȳ in (7) and taking the limit λ 1, the solution to Eq. (7) is w = R{exp[2π (iξ − y/λ)]}, as clearly shown in Fig. 4(left). This asymptotic scaling is confirmed by the analysis given in Sect. 5. No convective transport is at work: the physical balance is between the wall-normal and the streamwise viscous diffusion. As the wavelength increases to λ = 10, Fig. 4(right) shows that the velocity tends to a constant value along y.  (7), are negligible. Like the flows at order-one Reynolds number shown in Fig. 3, the spanwise flow occupies the entire channel, with the viscous region becoming thinner as λ/R decreases. For the case with λ = 1000 (Fig. 5, right), the spanwise velocity is finite at the centreline.
Profiles for λ R resemble the classical Stokes solution for high positive U and for negative U . Trends similar to the exponentially decaying profiles shown in Fig. 4 are found when λ 1. However, the profiles vary significantly from the classical Stokes solution when U ≈ 1 or smaller (and positive) and λ is sufficiently larger than 1 and smaller than R. Figure 6 shows four sets of profiles for λ = 50, R = 1000, and 0.5 ≤ U ≤ 1.25. The trends are similar to the profiles of the Stokes layer only in the upper portion of the viscous layer, while at lower locations the trends show oscillatory behaviour which is distinctly different from the Stokes layer. For example, profiles for U = 0.75 with w = 0.7 at y = 0 may decay and change their curvature as y increases without crossing the w = 0 line for the whole channel extent. The spanwise velocity may therefore be positive along the whole channel, which is never the case for the Stokes layer.
Solving the governing equation (7) represents a simple numerical exercise. However, it is clear from the numerical results presented that restricting the analysis to a computational endeavour severely limits the understanding of the physical problem. Therefore, asymptotic methods, i.e. the boundary-layer theory in Sect. 5, the WKBJ theory in Sect. 6, and the Langer theory in Sect. 7, are employed. These analyses are useful because approximate analytical solutions of (7) are obtained and because insight on the physics is gained, which cannot be revealed either through the full numerical solution or the PCF analytical expression (10). In particular, the asymptotic approach quantifies the thickness of the spanwise viscous layer, highlighting the physical balance very near the wall, and explains the occurrence of the wiggly behaviour for wave speeds comparable with the bulk velocity, shown in Fig. 6 (refer to Sect. 6). The theoretical analysis precisely identifies the parameter range for this turning-point regime, i.e. R −1 λ R and 0 ≤ U ≤ 3/2, and quantifies the thickness of the thin turning-point layer and of the other two order-one regions which confine this layer. The physical balances in these three layers are revealed, which explains the mathematical forms of their asymptotic solutions. The solid lines in Fig. 6, located at y = 1 − √ 1 − 2U/3, indicate the turning-point location. It will be further shown that the asymptotic analysis is also useful for the design of the proposed micromixer based on the travelling waves. The boundary-layer theory indeed identifies the cases where the spanwise flow is confined to a very thin wall-bounded layer, which are clearly not candidates as efficient mixers because the bulk flow, where the mixing is required, is largely unaffected by the wall motion.

Boundary-layer theory
To expedite the analysis, the asymptotic parameter ε ≡ λ/R is defined, which can be written as As also clear from multiplying each term of Eq. (5) by A * U * /λ * , ε represents the ratio between the wall-normal viscous effects and the convection effects due to the transport of the Poiseuille flow on the streamwise gradient of the spanwise flow. When U = O(1) these convection effects are also comparable to the unsteadiness due to the wave motion. In terms of ε, Eq. (7) is written as The cases for ε 1 are studied in Sect. 5.1 and the cases for ε = O(1) and ε 1 are studied in Sect. 5.2.

The small-ε regime: λ R
In the limit ε 1, the solution of Eq. (11) is W = 0, which satisfies the boundary condition W (1) = 0, but not the wall condition W (0) = 1. A boundary layer therefore exists in the proximity of the wall, where the solution varies rapidly. When y 1, the coordinate y is rescaled as Y = y/δ = O(1), where δ = ε β is the boundary-layer thickness and β is an unknown positive number. Within the boundary layer, Eq. (11) becomes where W (Y ) = W (y). There are two possible distinguished limits in (12): β = 1/2 and β = 1/3. This scenario is typically encountered in nested boundary-layer problems (refer to Bender and Orszag [38, example 6 on p. 453]), where an inner-inner boundary layer must exist on one side of the domain because the inner boundary layer does not satisfy the boundary condition. However, only one boundary layer exists in our case and the two β values correspond to distinguished flow regimes which match asymptotically in the parameter space (λ, R). Equation (12) is written as The case with β = 1/2 is studied first. The terms in (13a) that are proportional to Y and Y 2 may be neglected when the y diffusion balances either the convection term due to the wave transport, i.e. U = O(1), or the spanwise viscous diffusion term, i.e. λR = O(1), or both. Equation (13a) can then be simplified to The solution to Eq. (14) is where the branch cut in the square root is taken as the negative real axis. The convection effects brought about by the travelling wave are balanced by the viscous effects along the x and y directions, and the boundary-layer thickness which identifies the oscillating-wall regime. The viscous effects along x are negligible. Expression (16), combined with (6), leads to the solution of the Stokes problem for a flat plate oscillating sinusoidally in time beneath a stationary fluid. This is because in this high-frequency limit the streamwise-dependent term in the exponent in (6) is negligible with respect to the unsteady term in (6). The boundary-layer thickness is . This flow was also studied by McHale et al. [41] for a shear-horizontal shear wave travelling below a stationary fluid as an idealized model of a quartz crystal microbalance used, among many applications, for in situ monitoring of film deposition and analysis of polymer coatings.
When U (λR) −1 , expression (15) simplifies to which identifies the regime where the x and y viscous diffusion effects balance each other to give a boundary-layer thickness δ = O(λ) (refer to the numerical solution in Fig. 4(left)). This flow is shown to be steady by rescaling time by the average time taken by a fluid particle to cover a distance λ * along x * , i.e.t = t * U * b /λ * . The rescaling is necessary because the wave speed is now small. It follows that the exponent in (6) is 2π i x − Ut , which simplifies to 2π ix for small U .
In Eq. (13a), the term involving Y 2 is always negligible with respect to the term involving Y . However, if neither of the other two terms multiplying W (the convective term containing U and the x-diffusion term containing (λR) −1 ) is O(1) and at least one balances the term ε 1/2 Y , there is no term to balance the y viscous diffusion term. This may , the term proportional to Y alone balances the y viscous diffusion term, leading to the steady-wave case studied by Viotti et al. [42]. The scaling with β = 1/2 therefore breaks down because the term proportional to W is always smaller than the wall-normal viscous term and the scaling with β = 1/3 applies.
The β = 1/3 case is now investigated. The solution to (13b) is written in terms of the Airy function of the first kind as which was found by Quadrio and Ricco [32]. Expression (18) (18), may be neglected when λ R −1/2 , namely when the last term in (13b) is negligibly small. When U (λ/R) 1/3 and λ R −1/2 , the steady-wave regime with no x diffusion effects, studied by Viotti et al. [42], is retained as the last two terms in the argument of the Airy function in (18) are negligible. The steady-wave regime with streamwise diffusion is found when U (λ/R) 1/3 and λ = O(R −1/2 ). For large U , the asymptotic limit of (18) has been studied by Quadrio and Ricco [32]: the classical solution (16) to the classical Stokes problem is recovered. Analogously, the thin-layer steady-wave regime solution (17) The following asymptotic formula therefore applies: where c k are given in Abramowitz and Stegun [43, 10.4.59 on p. 448]. By substituting (19) into Eqs. (20) and (18), one finds Expression (21) reduces to (17) by expressing the algebraic term in (21) through a Taylor series expansion with respect to λ 1, i.e.
It is thus shown that the asymptotic solutions corresponding to the two distinguished limits of (12), for β = 1/2 and β = 1/3, match asymptotically in the parameter space. Figure 7 shows a schematic of the asymptotic regimes for λ R. When U 1 and λ = O(U −1/2 ), Eq. (7) can be simplified in similar fashion to the λ R regime. It reduces to whose solution is (15). The two simplified solutions, i.e. (16) and (17), are obtained in the limits λ U −1/2 and λ U −1/2 , respectively. Figure 8 shows a schematic of the asymptotic regimes for λ = O(R). Note that a boundary layer exists only when the convection due to the wall motion is negligible (areas enclosed by dashed line), whereas the viscous effects extend throughout the entire channel when the PCF solutions apply. The trapezoid in the top right corner of Fig. 8 indicates that the classical Stokes layer occurs when λ U −1/2 , and hence, for every value of U , there exists a very large λ above which the classical Stokes layer is always recovered. (λR) −1 , respectively. We close this section by pointing out that the boundary-layer analysis reveals no information on the physics behind the peculiar profiles shown in Fig. 6. This is studied through the WKBJ theory in Sect. 6.

Turning-point solution by matched asymptotic expansion
In Eq. (11), the streamwise diffusion is negligible when the last term in the parenthesis is smaller than all the other terms. When U = O(1), this occurs when λ R −1/2 . Equation (11) simplifies to εW (y) − iP(y)W(y) = 0, where P(y) = −3π y 2 + 6π y − 2πU . In this limit, the WKBJ asymptotic solution (25) is not valid for 0 ≤ U ≤ 3/2 because P(y) = 0 when y = y 0 = 1 − √ 1 − 2U/3, which causes the WKBJ solution (25) to become unbounded there. In this case, a turning point occurs at y = y 0 . The WKBJ theory reveals that the profiles shown in Fig. 6 belong to this category and is able to define precisely the parameters for which this behaviour occurs. For 0 < U 1, the turning point is located at y 0 ∼ U/3. The relationship between y 0 and U is shown in Fig. 12(left) of Appendix C.
Following Bender and Orszag [38], a composite WKBJ expansion can be constructed in the half channel 0 ≤ y ≤ 1, i.e. where and the branch cut is chosen so that The solution W m (y) in (27) occupies y − y 0 > ε 1/3 and is bounded above by the centreline, while W w (y) in (27) occupies y − y 0 < −ε 1/3 and is bounded below by the channel wall. The solution W 0 (y) in (27) is about the turning point at y = y 0 . W c 0w (y) is the common part matching the behaviour near the turning point to the WKBJ behaviour near the wall, and W c 0m (y) is the common part matching the behaviour near the turning point to the WKBJ behaviour near the centreline. These common parts and W 0 (y) are derived in Appendix A. As detailed in Appendices B and C, the constants a m , b m , a 0 , b 0 , a w , and b w are determined by matching the solutions asymptotically and applying the boundary conditions W(0) = 1 and W (1) = 0.
As shown in Appendix C, the asymptotic solution (27) shows excellent agreement with the numerical solution. The changing nature of the flow on either side of the turning point can be explained by the changes in the dominant balances in (11). In this parameter range, the dominant behaviour close to the wall is governed by a balance between the unsteady convection and the wall-normal viscous stresses. Moving away from the wall, the streamwise Poiseuille-driven convection increases. Between the wall and the turning point, the streamwise Poiseuille-driven convection remains smaller than the convection resulting from unsteady wall wave forcing. At the turning point, a constant-viscous-stress balance exists between the streamwise Poiseuille-driven convection and the convection due to the unsteady wave forcing. Above the turning point the convection due to the Poiseuille flow is more significant than the contribution due to the wall wave forcing and the y viscous diffusion is present again.

Langer theory
An asymptotic approximation which is alternative to the WKBJ turning-point solution can be constructed via the Langer transformation: η = η(y) and v = η (y)W(y) [44]. This is uniformly valid between the walls and the centreline. Equation (26) can be written as For ε 1 and Δ = O(1), the approximate solution to (30) can be calculated from the related equation Following Langer [36,37], we choose iP/ η 2 = η, and hence, upon integration, where κ(y) is given by (29). A uniform asymptotic approximation to W(y) is given by where and η 0 = ε −1/3 η(0) and η 1 = ε −1/3 η(1). The real and imaginary parts of the Langer solution (31), shown in Fig. 10 for two cases investigated in Sect. 6.1, show excellent agreement with the numerical profiles.

Application to flow mixing
It is our hope that our theoretical results will spur experimental research for which the fluid-structure interaction driven by shear-horizontal waves is important. In this section, we put forward ideas for their possible future applications.
-Active mixing of laminar flows by surface acoustic waves. Rayleigh waves, i.e. wall-normal displacement streamwise-travelling acoustic waves, have been utilized for microfluidic mixing. Sritharan et al. [1] and Tseng et al. [45] experimentally verified that Rayleigh waves can mix cofluent streams with very different passive scalar concentrations. However, in liquids these waves suffer from severe energy dissipation due to the compression waves engendered by the wall-normal displacement (leaky-wave phenomenon). Also, although the induced small-scale secondary recirculatory motion is beneficial for mixing, it may interfere with the smoothness of the streamwise flow and create additional pressure gradients and therefore additional losses. We instead propose to use shear-horizontal waves to mix the cofluent streams studied by Sritharan et al. [1] and Tseng et al. [45], which, to the best of our knowledge, have never been employed in microfluidic mixing. A schematic of the microfluidic SAW mixer is shown in Fig. 11(left).
The main advantages over the Rayleigh waves would be (i) less energy dissipation (higher efficiency) because the shear-horizontal waves do not suffer from acoustic streaming energy loss due to the absence of wall-normal motion, as thoroughly discussed by Lange et al. [13] and Rocha et al. [46], which implies (ii) mixing along longer streamwise distances [46]. Furthermore, (iii) as the mixing occurs through the spanwise velocity, the streamwise flow remains smooth and no additional induced pressure gradients must be accounted for. A fourth advantage is that (iv) the spanwise waves would be better mixers than two-dimensional Rayleigh waves given the concentration distribution at the inlet, shown in Fig. 11(left), which is uniform along the wall-normal direction, but strongly varying along the spanwise direction. Although the streamlines of the streamwise flow are unchanged when the SAWs are implemented, the mixing is required primarily along the spanwise direction where the concentration variations are most intense. The passive scalar equation to be solved is where θ is the passive scalar concentration (mass or temperature, for example), Pe = U * b h * /α * p is the Peclet number, and α * p is the diffusion coefficient for mass transfer and the thermal diffusivity for heat transfer. The spanwise waves w(x, y, z, t) would act, through the boxed term in (33), on the spanwise gradient of θ , which would be most intense at the upstream channel location, where the two cofluent streams start interacting. The Reynolds numbers in the experiments of Sritharan et al. [1] and Tseng et al. [45] are in the range 10 −2 −1 and the forcing wavelengths are comparable or smaller than the channel height, which leads to ε = λ/R in the range 1-100. This corresponds to the order-one-ε and large-ε regimes with U belonging to any columns of Fig. 8 because a wide range of wave speeds can be generated by wave interference, as proved by the standing-wave study by Ding et al. [8]. We finally note that the boundary-layer analysis of Sect. 5 is also here useful because it allows identifying the regimes where the spanwise flow is confined very near the wall, which are bound not to be efficient mixers as the mixing is required across the whole channel height. The turning-point regimes, for which the wave speed is comparable with the bulk velocity, could instead qualify as candidates for good mixing because the flow may extend along the channel height (refer to Fig. 6).
-Active mixing of laminar flows by electro-osmotic waves.
The travelling wave flow produced by (1) could also be generated by electro-osmotic waves and used for microfluidic flow mixing, as shown in Fig. 11(right). To the best of our knowledge, these waves have never been created through electro-osmosis, but proper design of an unsteady and spatially inhomogeneous electric field could achieve this purpose. They would be an unsteady and streamwise-modulated version of the waves suggested by Ajdari [18], a streamwise-modulated variant of the oscillatory flow studied numerically by Dutta and Beskok [47], or an optimized variant of the micromixer proposed by Oddy et al. [23], where the wall forcing would be spanwise and streamwise-modulated instead of simply oscillatory and spatially uniform. Extending the analyses of Dutta and Beskok [47] and Meisel and Ehrhard [25], the governing equation for the spanwise velocity with electro-osmotic effects included is subject to the no-slip boundary condition at the wall and to the zero-gradient condition at the centreline (we only consider half channel for simplicity). Here * d is the thickness of the Debye-Hückel layer, E * z = E * z,0 R e 2π i(x * −U * t * )/λ * is the spanwise electric field, ε * is the permittivity, and ζ * is the zeta potential (we have set the ionic energy parameter equal to unity, refer to Dutta and Beskok [47, p. 5098], and for simplicity assumed an exponentially decaying potential as in Meisel and Ehrhard [25] rather than more complex potential functions as in Afonso et al. [48] and Wang et al. [49]). The details of the formulation for the Stokes layer case, i.e. for u * = 0 and E * z = E * z,0 R e iω * t * , are found in Dutta and Beskok [47]. The above problem can be conveniently simplified under the assumption that Debye-Hückel layer * d is much thinner than the channel height and the viscous layers studied in Sect. 5. As also discussed by Qiao and Aluru [50], this hypothesis is amply verified as * d is very small, i.e. of the order of 100 nm [47,51], therefore about three orders of magnitude smaller than the channel height and two orders smaller than the viscous layers generated by the travelling waves. This means that the electric potential is confined in this very thin near-wall Debye-Hückel layer, while the bulk flow is electrically neutral and driven underneath by the electro-osmotic motion of the Debye-Hückel layer. The scaled form of (34) is where (1). Note that here w = w * /U * b as A * cannot be used for scaling like in Sect. 2 because it is not defined. It is found in the following through asymptotic matching. By defining the Debye-Hückel-layer coordinate Y d = y/δ d and velocity W d = w, the Debye-Hückel-layer equation is found at leading order whose solution is obtained by use of the boundary conditions W d (0) = 0 and W d (∞) = 0. The Helmholtz-Smoluchowski velocity is obtained as follows: In dimensional form, (38) becomes W * hs = − ε * ζ * E * z,0 /μ * R e 2π i(x * −U * t * )/λ * . This is the velocity that drives the bulk electrically neutral flow which we have studied in the previous sections. We can now quantify the amplitude of the wall travelling waves defined in (3), i.e. A * = −ε * ζ * E * z,0 /μ * . In summary, the bulk spanwise flow, which is relevant for mixing, is thus described by Eq. (4) and driven by the unsteady and streamwise-modulated Helmholtz-Smoluchowski velocity W hs in (38), while the spanwise velocity in the very thin Debye-Hückel layer is brought to zero at the wall through the no-slip condition. The composite solution is w c (x, y, t) = w(x, y, t) As for the SAW mixer, the passive scalar Eq. (33) is to be solved. Similar electro-osmotic microfluidic mixers have been studied by Sasaki et al. [52] and Huang et al. [53]. Sasaki et al. [52] employed meandering electrodes to mix two microstreams and stress the importance of obtaining analytical results for the fluid flow in order to optimize the mixing performance, while Huang et al. [53] generated in-plane microvortices to prove that up to 30-fold mixing enhancement can be achieved compared to mixing due to diffusion only. The spanwise spatial pattern displayed in Fig. 11(right) can be achieved by utilizing thin strips of different glass coatings and spatially modulated electric fields [18,20]. Typical frequencies of electro-osmosis actuators are of the order of 10 Hz [23,53] and the streamwise length of the microelectrode arrays, which define the forcing wavelength, can be in the range of 100-500 µm. The flow parameters therefore correspond to ratios ε = λ/R in the range of 0.1-1 and to U of order unity. The small-ε and order-one-ε regimes characterize these flows. The turning-point WKBJ analysis and the Langer analysis are thus relevant for these flow regimes, where the wave speeds are positive and comparable with the bulk velocity. The streamwise wall-shear stress would play a crucial role in the spanwise flow dynamics and thus for the mixing performance. The design proposed in Fig. 11(right) would generate a square streamwise-travelling wave of spanwise velocity, a micro-scale analogous version of the pipe-flow wave employed by Auteri et al. [28]. As explained in Sect. 9, thanks to the linearity of the problem, Fourier decomposition will allow the use of our sinusoidal-wave results to construct the base flow for this mixing problem.
As mentioned in the Introduction, the shear-horizontal waves have also been used for turbulent drag reduction. For this application, the phase speeds span a wide range of values from null to several times the bulk streamwise velocity. It therefore follows that the regime with ε = λ/R 1 properly describes these flows. It also occurs that λ R −1/2 and therefore the streamwise viscous diffusion due to the waves is negligible. We reiterate that, in fully developed turbulent channel flows, the laminar flow solutions may only be representative of the spanwise-averaged velocity profile and only under strict conditions of the wave parameters, as amply discussed in Quadrio and Ricco [32]. It is interesting to note that the drag-increase regime discovered by Quadrio et al. [31] for a specific range of positive phase speeds is included within the range for which the turning-point regime occurs. Further investigation is also required to generalize the stability analyses for the classical Stokes layer [54,55] to the travelling wave case.

Summary
Channel flows with spanwise wall forcing consisting of in-phase sinusoidal travelling waves of spanwise wall velocity have been investigated. A novel three-dimensional time-dependent solution of the incompressible Navier-Stokes equations is constructed, with the solution represented as a linear combination of complex-to-complex parabolic cylinder functions. As reliable numerical solutions of the complex-to-complex parabolic cylinder function are currently unavailable, asymptotic solution methods have been employed to investigate flow variations due to the Reynolds number R, the scaled phase speed U , and the scaled wavelength λ. Only sinusoidal shear-horizontal wall forcing has been considered. However, flows produced by more general wall forcing can be expressed by a linear combination of our solutions with each term multiplied by its Fourier coefficient.
Asymptotic methods, i.e. the boundary-layer and the WKBJ theories, have been utilized to study the flow. The underlying flow physics, revealed by the dominant balances in the governing equation, is gained by using these asymptotic theories and cannot be obtained through the exact analytical solution (10) or the numerical solutions computed in Sect. 4. While the boundary-layer method and the WKBJ approach both produce excellent approximations to the flow, there are particular advantages for each method. The simplicity of the boundary-layer solutions compared to the WKBJ solution is noticeable, aided by the outer solutions being identically zero. The boundarylayer method readily enables the determination of the boundary-layer thickness, which is not available if the WKBJ method is employed. Furthermore, the link between the different physical effects is elucidated better utilizing the boundary-layer approach. There are also advantages of the WKBJ approach over the boundary-layer method. Without further scaling (introducing a second boundary layer close to y = 2), the WKBJ method readily generates solutions which are valid across the full channel.
The WKBJ solution also identifies the turning-point behaviour for 0 ≤ U ≤ 3/2, R −1/2 λ R, which is not explained by the boundary-layer method. As the standard WKBJ solution (25) is unbounded near the turning point, solutions have been found by a WKBJ composite expansion and the Langer method. While the Langer solution is simpler, the composite WKBJ expansion has the benefits of explicitly determining the turning-point location and of quantifying the thickness of this thin layer. This is important because the physics changes there. The WKBJ theory shows that, when the streamwise diffusion effects are negligible, in the turning-point layer the Poiseuille convection and the convection due to the waves cancel out so that the wall-normal viscous stresses are constant. The high accuracy of the asymptotic solutions is quantified by comparing them with the numerical profiles in Appendix D.
We have finally presented ideas on how to generate the travelling waves for microfluidic flow mixing via surface acoustic forcing and electro-osmotic actuation. In the latter case, matched asymptotic expansion has been useful to obtain the Helmholtz-Smoluchowski velocity that drives the bulk spanwise flow.
linearly independent solutions of the Airy equation, Ai(y) and Ai e 2iπ/3 y , are utilized to facilitate the matching with the outer WKBJ behaviour above and below the turning-point region.
The turning-point behaviour given by (28b) is not valid when the quadratic term in (40) is as large as the linear term and hence cannot be neglected. This occurs when y 0 is close to the centreline because P (1) = 0. Asŷ = O(1), the linear and the quadratic terms balance when y 0 − 1 = O(ε 1/3 ), which is found by use of the definition of P(y). Inside this region, the quadratic term in the expansion (40) cannot be neglected and, rather than the Airy equation, the local behaviour is governed by an ordinary differential equation whose solutions are parabolic cylinder functions. Consistent with our earlier approach of not developing an algorithm for complex-to-complex parabolic cylinder functions, the particular case with the turning point close to the centreline will not be discussed further.

Appendix B: Matching for WKBJ turning-point composite expansion
When the solution about the turning point is given by (28b), W 0 and W m are matched over a region ε 1/3 < y − y 0 < ε 1/5 , which gives rise to the connection formulae: Similarly, W w and W 0 are matched over a region −ε 1/5 < y − y 0 < −ε 1/3 , which leads to a w = a 0 ε 1/12 iP (y 0 ) The derivation of the connection formulae is usually given for cases where iP(y) is a strictly real valued function and is found in Bender and Orszag [38]. The derivation is extended now to the problem of a strictly complex-value iP(y). For y > y 0 , the Taylor expansion of P(y) about y = y 0 is used in the integral for κ(y), and hence Consequently for y → ∞, the local solution about the turning point has the form Comparing the coefficients of the exponential terms in (43) and (45) gives rise to the connection formula (41), while similar matching between (28c) and (28b) gives rise to the connection formula (42).

Appendix C: WKBJ turning-point composite solutions
In this Appendix, the possible flow configurations that can arise due to the relationship between the matching regions, the wall and the centreline, are discussed. As illustrated from left to right in Fig. 12(right), one finds: -Case 1 (y 0 < ε 1/5 and y 0 > 1 − ε 1/5 ): W 0 is valid across the full half channel, with W(y) = W 0 (y), where A 0 , B 0 , A 1 , and B 1 are defined in (46). -Case 2 (y 0 < ε 1/5 and y 0 < 1 − ε 1/3 ): W 0 is valid near the wall and the turning point, while W m is valid near the centreline. Here W(y) = W m (y) where A 0 and B 0 are defined in (46). -Case 3 (y 0 > ε 1/3 and y 0 > 1 − ε 1/5 ): W 0 is valid near the turning point and the centreline, while W w is valid near the wall. Here W(y) = W w (y) + W 0 (y) − W c 0w (y), with where A 1 and B 1 are defined in (46). -Case 4 (y 0 > ε 1/3 and y 0 < 1 − ε 1/3 ): W 0 is valid near the turning point, while W w is valid near the wall and W m is valid near the centreline. Here W(y) is given by (27), with .
When y 0 < ε 1/5 and y 0 > ε 1/3 the wall is contained within the matching region between W w (y) and W 0 (y), and when y 0 < 1 − ε 1/3 and y 0 > 1 − ε 1/5 the centreline is contained in a matching region between W m (y) and W 0 (y). Of the four cases to be described, the agreement between the WKBJ turning-point solution and the numerical solution is least good in case 1, where the imaginary parts of the composite turning-point solution and the numerical solution show a relative error of 26% at the turning point. This is because the local turning-point solution W 0 in (28b) has to be employed to describe the behaviour close to the wall, the turning point, and near the centreline. It is expedient to investigate this case even if it is less accurate than the subsequent cases, because the coefficients (46) are also present in the more complicated expansions. For cases 2, 3, and 4, the composite expansions are in excellent agreement with the numerical solutions.

Appendix D: Errors in numerical and asymptotic solutions
In order to assess the relative accuracy of the asymptotic approximations, the error measure is calculated, where W is the approximate solution generated by each approach and W num is the numerical solution calculated utilizing the method outlined in Sect. 4. Profiles are compared on a uniform mesh with 8000 uniformly spaced points over the range 0 ≤ y ≤ 1. In addition to determining the errors associated with each asymptotic method, the error associated with the numerical solution is estimated via a comparison with a second numerical solution calculated on a grid with a mesh spacing of exactly double that of the first numerical solution. Figure 14 -The error in the numerical solution increases as the value of ε decreases. This is because smaller values of ε produce higher wall-normal velocity gradients close to the channel wall and, as ε decreases, the approximation of these gradients on a uniform grid becomes less accurate. -For ε 10 −2 , the WKBJ method outperforms the boundary-layer approximation. In this region the errors associated with the two methods decay at comparable rates as the value of ε decreases. The decay of E as ε decreases is fully expected as ε is the small asymptotic parameter. -For ε 10 −2 , the apparent errors associated with the asymptotic solutions start to increase as the value of ε decreases. This is because the asymptotic solutions are compared with the numerical solution. The asymptotic solutions tend to the exact solution as ε → 0, whereas the numerical error increases, as discussed in the first bullet point.
-For very small ε, the errors associated with the WKBJ method and the boundary-layer method are very similar, which suggests that the two asymptotic methods are in close agreement both with each other and with the exact solution. This is to be expected for small ε.
To investigate the relative accuracy of the WKBJ turning-point expansion and the Langer method when a turning point is present, the error associated with each method is calculated for U = 1.125 (y 0 = 0.5) and a range of values of ε, with the errors shown in Fig. 14(right). In terms of the hierarchy of flow solutions with turning points illustrated in Fig. 12(right), the solution about the turning point is valid across the full half channel for ε 2 −5 (Sect. 6.1, case 1), while the full five-term composite expansion (27) is valid for ε 2 −3 (Sect. 6.1, case 4). The limits of validity for case 1 and case 4 are denoted by vertical lines in Fig. 14(right).
-The Langer method outperforms the turning-point expansion (case 4) across a wide range of values of ε. The decay rate of both methods is similar as the value of ε decreases. -Throughout the region of common validity (2 −5 ε 2 −3 ), the error associated with the five-term composite expansion (case 4) is smaller than that observed with just the turning-point solution (case 1). This is consistent with the behaviour shown for ε = 0.05 in Sect. 6.1. -For ε 10 −2 , the asymptotic methods are more accurate than the numerical method, as shown in Fig. 14. -For very small ε, good agreement between the turning point expansion and the Langer solution is obtained as both tend to the exact solution. -For ε ≈ 1, the errors associated with the turning point expansion and the Langer method are O(1), and hence these asymptotic methods are of no use in this regime. This is in marked contrast to the earlier analysis for U = 10, where the errors associated with the boundary-layer method (E = 4.0 × 10 −3 ) and the WKBJ approximation (E = 2.4 × 10 −4 ), when ε = 1 indicate that these methods may still be of use.