2.1. Governing equations

The continuity equation and Navier–Stokes equation for a viscous compressible fluid can be written as

(2.1a) \begin{equation} \begin{cases}\dfrac{\partial \rho }{\partial t}+\boldsymbol{\nabla }\boldsymbol{\cdot }\left(\rho \boldsymbol{v}\right)=0,\\[11pt] \rho \dfrac{\partial \boldsymbol{v}}{\partial t}+\rho \left(\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\right)\boldsymbol{v}=-\boldsymbol{\nabla }\!p+\mu \boldsymbol{\nabla} ^{2}\boldsymbol{v}+\left(\mu _{b}+\dfrac{\mu }{3}\right)\boldsymbol{\nabla }\left(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}\right)\!, \end{cases} \end{equation}

where $\rho , \boldsymbol{v}=(u,v,w), p, \mu$ and $\mu _{b}$ are the mass density, fluid velocity, fluid pressure, shear viscosity and bulk viscosity, respectively. For the acoustic field in the present SAW-driven cavity, the acoustic Mach number remains much smaller than unity, and the entropy variation induced by viscous dissipation and heat transfer is negligible (Li et al. Reference Li, Liang, Kabla, Zhang and Yang2023; Marques & Tomaz Da Silva Reference Marques and Tomaz Da Silva2025).Therefore, the thermodynamic relation between pressure and density is characterised by

(2.1b) \begin{equation} c^{2}=\left(\frac{\partial p}{\partial \rho }\right)_{s}\!\!,\end{equation}

where $c$ is the isentropic speed of sound in the liquid.

2.2. Perturbation method

The difficulty of solving the nonlinear terms in (2.1a ) lies in the large time scale difference between the applied periodic ultrasonic oscillations with frequency, ranging from 1 to 100 MHz, and the induced nonlinear streaming flow with characteristic time spanning from tenths of a second to several minutes (Nama et al. Reference Nama, Barnkob, Mao, Kähler, Costanzo and Huang2015). Since the MHz acoustic excitation introduces only small perturbations to the fluid, the physical quantities are expanded in terms of a small non-dimensional parameter $\epsilon$ , which characterises the relative acoustic amplitude. For the present SAW-driven cavity, $\epsilon$ can be estimated by the acoustic Mach number based on the amplitude of the imposed velocity at the boundary, i.e.

(2.2) \begin{equation}\epsilon \sim \textit{Ma}=W/c_{0},\end{equation}

where $\textit{Ma}$ is the acoustic Mach number, W is the characteristic boundary velocity amplitude induced by the SAW excitation, as defined in (2.5), and $c_{0}$ is the speed of sound in water. For the present operating condition, $W\approx 0.01$ m s⁻¹, giving $Ma\approx 6\times 10^{-6}\ll 1$ . Therefore, following a small-amplitude perturbation approach (Köster Reference Köster2007), the physical quantities are decomposed into the zeroth-order component (quiescent steady flow), the first-order component (harmonic oscillatory response), the second-order component (giving rise to the time-averaged acoustic streaming when averaged over an acoustic cycle) and higher-order small terms

(2.3) \begin{equation} \begin{cases} \boldsymbol{v}=\boldsymbol{v}_{0}+\epsilon \boldsymbol{v}_{1}+\epsilon ^{2}\boldsymbol{v}_{2}+O(\epsilon ^{3}),\\[3pt] p=p_{0}+\epsilon p_{1}+\epsilon ^{2}p_{2}+O(\epsilon ^{3}),\\[3pt] \rho =\rho _{0}+\epsilon \rho _{1}+{\epsilon ^{2}}\rho _{2}+O(\epsilon ^{3}), \end{cases} \end{equation}

where the subscripts 0, 1, 2 refer to the zeroth-order, first-order and the second-order components; $O(\epsilon ^{3})$ represent third-order small terms. Contributions beyond second order are assumed to be negligible and thus discarded.

Under the isentropic assumption in (2.1b ), substituting (2.3) in to (2.1) and collecting only the first-order terms, we arrive at the linear governing equations

(2.4) \begin{equation} \begin{cases} \dfrac{\partial \rho _{1}}{\partial t}+\rho _{0}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}_{1}=0,\\[11pt] \rho _{0}\dfrac{\partial \boldsymbol{v}_{1}}{\partial t}=-\boldsymbol{\nabla }\!p_{1}+\mu \boldsymbol{\nabla} ^{2}\boldsymbol{v}_{1}+\left(\mu _{b}+\dfrac{\mu }{3}\right)\boldsymbol{\nabla }\left(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}_{1}\right)\!,\\[11pt] p_{1}=c^{2}\rho _{1}. \end{cases} \end{equation}

Figure 1.

Schematic of the studied acoustofluidic models with different SAW configurations. (a) A microfluidic cavity whose bottom and top walls (or only the bottom wall) are piezoelectric substrates, each configured with a pair of opposing IDTs. The SAWs generated by these IDTs couple into the liquid and superpose. (b) A fluid chamber is configured on a substrate where two pairs of IDTs are orthogonally arranged to generate a 2-D SSAW field inside the fluid. (c) The fluid region is abstracted as a cubic domain. Red arrows indicate the direction in which periodic boundary conditions are applied, and blue arrows mark the positive direction of velocity boundary conditions.

Figure 1. Long description

Panel (a) shows a fluid cavity between parallel piezoelectric substrates. Pairs of opposing interdigital transducers on the substrate surfaces generate counter-propagating surface acoustic waves that leak into and superpose within the liquid; the upper substrate may be omitted for single-sided actuation. Panel (b) shows a chamber on one substrate surrounded by two orthogonal transducer pairs, producing a two-dimensional standing surface acoustic wave. Panel (c) abstracts the liquid as a cuboid of square base length L and height H. The x and y axes lie in the substrate plane and z is vertical. Red arrows mark periodic boundary conditions on the four side faces; blue arrows mark velocity boundary conditions on the top and bottom faces.

2.3. Boundary conditions

Typical SAW-based microfluidic systems consist of piezoelectric substrates, IDT pairs and microchannels. By applying an alternating current signal to the IDTs, SAWs are actuated on the surface of the piezoelectric substrate and propagate perpendicular to the orientation of the IDT fingers. The SAWs can be introduced into both continuous flows (figure 1a ) and quiescent fluids (figure 1b ), and the IDTs can be arranged in various configurations to generate controllable acoustofluidic field within the fluid domain. These devices can be abstracted to model as illustrated in figure 1(c), where the fluid domain is constrained in a cubic cavity with a height H and a square base of edge length L. The Eulerian reference frame $O(x,y,z)$ is placed at the centre of the bottom boundary of the liquid. Both the x and y axes are parallel to the plane of the substrates, and the axis z is perpendicular to the substrate plane. The upper and bottom boundaries are in direct contact with the piezoelectric substrates, through which leaky SAWs couple into the fluid. When a pair of IDTs is configured in parallel on opposite sides of the fluid cavity, counterpropagating SAWs form a 1-D SSAW. When two pairs of IDTs (four IDTs in total) are positioned orthogonally on a substrate (as shown in figure 1(b)), two orthogonal 1-D SSAWs interfere to form a 2-D SSAW field. Building on these traditional set-ups, we develop a model, for the first time, considering the SSAW fields in three dimensions, where 2-D SSAWs are formed simultaneously on both the upper and lower boundaries of the fluid region. With this framework, we systematically investigate the resulting 3-D acoustofluidic fields induced within the fluid domain by dual-sided 2-D SSAWs.

As the SAWs are spatially periodic along their propagation paths on the substrate (i.e. in the x- and y- directions), the resulting acoustofluidic field inherits this spatial periodicity when the substrate and fluid are assumed isotropic and the system operates in steady state. To utilise this property, we choose the edge lengths L of the square base of the domain to be an integer multiple of the SAW wavelength $\lambda _{\textit{SAW}}$ , thereby any physical quantity $Q(x,y,z,t)$ , such as the velocity $\boldsymbol{v}=(u,v,w)$ , pressure $p$ and density $\rho$ , takes the same value on opposing boundaries of the fluid domain in the x- and y- directions. Thus, periodic boundary conditions can be applied to four side boundaries of the fluid domain. Notably, setting the edge lengths of the base face of the cubic fluid domain to the same value L is not a requirement. This choice is made merely for convenience when expressing the periodic boundary conditions

(2.5) \begin{equation} \begin{cases} Q\left(-\dfrac{L}{2},y,z,t\right)=Q\left(+\dfrac{L}{2},y,z,t\right)\!,\\[15pt] Q\left(x,-\dfrac{L}{2},z,t\right)=Q\left(x,+\dfrac{L}{2},z,t\right)\!. \end{cases} \end{equation}

Because the top and bottom faces of the fluid domain (the surfaces normal to the z-direction) are in direct contact with the piezoelectric substrates, continuity requires that the substrates and the adjacent fluid share the same velocity. Therefore, as illustrated in figure 1(c), both the upper and lower surfaces of the fluid region are described by the velocity boundary conditions. Since the SAW excitation frequency is typically in the MHz range, which is several orders of magnitude higher than the characteristic time scales of particle migration and acoustic streaming, substrate vibration is treated as a steady, time-harmonic pattern that does not vary on the slow time scale of particle movement. For a single-frequency harmonic oscillation, the substrate surface displacement $\boldsymbol{u}_{\textit{sub}}$ and velocity $\boldsymbol{v}$ are related by $\boldsymbol{v}=i\omega \boldsymbol{u}_{\textit{sub}}\textit{Re}(\exp (i\omega t))$ . Therefore, imposing either a displacement or a velocity boundary condition is equivalent.

The velocity boundary condition that used to prescribe SAW propagation on the substrate has been reported in Li et al. (Reference Li, Liang, Kabla, Zhang and Yang2023). The out-of-plane (along the z-direction) vibration generated by orthogonally aligned IDTs can be expressed as

(2.6a) \begin{equation} \begin{cases} w_{+x,b}\left(x,y,0,t\right)=W_{+x,b}\boldsymbol{\cdot }\textit{Re}(\exp (-C_{\textit{d},x}x-i(k_{\textit{SAW},x}x-\omega t+\varphi _{+x,b}))),\\[3pt] w_{-x,b}\left(x,y,0,t\right)=W_{-x,b}\boldsymbol{\cdot }\textit{Re}(\exp (C_{\textit{d},x}x-i(-k_{\textit{SAW},x}x-\omega t+\varphi _{-x,b}))),\\[3pt] w_{+y,b}\left(x,y,0,t\right)=W_{+y,b}\boldsymbol{\cdot }\textit{Re}(\exp (-C_{\textit{d},y}y-i(k_{\textit{SAW},y}y-\omega t+\varphi _{+y,b}))),\\[3pt] w_{-y,b}\left(x,y,0,t\right)=W_{-y,b}\boldsymbol{\cdot }\textit{Re}(\exp (C_{\textit{d},y}y-i(-k_{\textit{SAW},y}y-\omega t+\varphi _{-y,b}))), \end{cases} \end{equation}

and

(2.6b) \begin{equation} \begin{cases} w_{+x,t}\left(x,y,H,t\right)=W_{+x,t}\boldsymbol{\cdot }\textit{Re}(\exp (-C_{\textit{d},x}x-i(k_{\textit{SAW},x}x-\omega t+\varphi _{+x,t}))),\\[3pt] w_{-x,t}\left(x,y,H,t\right)=W_{-x,t}\boldsymbol{\cdot }\textit{Re}(\exp (C_{\textit{d},x}x-i(-k_{\textit{SAW},x}x-\omega t+\varphi _{-x,t}))),\\[3pt] w_{+y,t}\left(x,y,H,t\right)=W_{+y,t}\boldsymbol{\cdot }\textit{Re}(\exp (-C_{\textit{d},y}y-i(k_{\textit{SAW},y}y-\omega t+\varphi _{+y,t}))),\\[3pt] w_{-y,t}\left(x,y,H,t\right)=W_{-y,t}\boldsymbol{\cdot }\textit{Re}(\exp (C_{\textit{d},y}y-i(-k_{\textit{SAW},y}y-\omega t+\varphi _{-y,t}))), \end{cases} \end{equation}

where $k_{\textit{SAW},x}$ and $k_{SAW,y}$ are the SAW wavenumbers along the $x$ - and $y$ -directions on the substrate, respectively; $C_{\textit{d},x}$ and $C_{\textit{d},y}$ are the corresponding attenuation coefficients of the SAWs along the $x$ - and $y$ -directions; $\varphi$ and $W$ denote the initial phase and velocity amplitude of each SAW component, respectively. The subscripts ‘b’ and ‘t’ denote the bottom and top substrates, and the subscripts ‘+’ and ‘–’ identify waves travelling in the positive and negative directions of the corresponding axis, respectively. For instance, $\varphi _{+x,b}$ is the initial phase of the wave travelling in the +x direction on the bottom substrate; $i$ is the imaginary unit, i.e. $i=\sqrt{-1} ; \textit{Re}(\boldsymbol{\cdot })$ denotes the real-part operator; $\omega$ is the angular frequency of the substrate vibration. It should be noted that (2.6) prescribes only the out-of-plane substrate vibration velocity $w$ . The in-plane vibration components $u$ and $v$ associated with the SAWs are neglected in the main text, since they are generally much smaller than the out-of-plane component in typical acoustofluidic systems (Jannesar & Hamzehpour Reference Jannesar and Hamzehpour2021; Li et al. Reference Li, Liang, Kabla, Zhang and Yang2023).

For the analytical development in the main text, two further simplifications are adopted. First, for isotropic substrates, or for specially arranged actuation conditions under which the SAWs propagating along different in-plane directions exhibit the same propagation characteristics (Shen et al. Reference Shen2024), one may take $k_{\textit{SAW},x}=k_{\textit{SAW},y}=k_{\textit{SAW}}$ and $C_{d,x}=C_{d,y}$ . Second, the attenuation coefficients $C_{\textit{d},x}$ and $C_{\textit{d},y}$ are neglected in the main text analysis. This approximation is commonly used in microfluidic applications, for which SAW attenuation over the characteristic propagation distance remains relatively weak. Specifically, with the attenuation coefficient reported by Nama et al. (Reference Nama, Barnkob, Mao, Kähler, Costanzo and Huang2015), $C_{d}=116\,\mathrm{m}^{-1}$ , the wave amplitude still retains approximately $70\mathrm{\,\% }$ of its initial value after propagating over $3000\,{\unicode[Arial]{x03BC}}\mathrm{m}$ . Accordingly, to keep the main derivation focused on the principal 3-D control mechanism, the attenuation terms are omitted in the main text.

Imposing velocity boundary conditions on both the top and bottom surfaces of the fluid domain is an approach of generalisation, rather than a specialisation. In practice, in addition to devices where the fluid channel is sandwiched between two piezoelectric substrates, there are also devices in which only one side of the channel is in contact with a vibrating substrate, while the opposite side interfaces a non-vibrating material. For such single-side-actuated devices, if the non-vibrating wall is acoustically rigid, we simply set the z-direction velocity amplitude in (2.6b ) to zero, i.e. $W_{j,t}=0$ or $W_{j,b}=0$ ( $j=+x,-x,+y,-y$ ), to enforce total reflection. If the non-vibrating wall is an acoustically soft material, the impedance boundary condition can be prescribed (Nama et al. Reference Nama, Barnkob, Mao, Kähler, Costanzo and Huang2015).

In (2.6), the reference phase is chosen as the mean of the initial phases of the two opposing travelling SAWs on the bottom substrate along the x-axis. To achieve independent control of the acoustic field in the x, y and z directions, SSAWs on the top and bottom substrates need to share the same vibration mode. Under this condition, the PNs and ANs of the 2-D SSAW fields on the top and bottom substrates are vertically aligned in the z-direction. Therefore, the eight initial phases of SAWs from the eight IDTs on the two substrates reduce to only four independent parameters

(2.7a) \begin{equation} \begin{cases} \varphi _{+x,b}=-\varphi _{x},\\[3pt] \varphi _{-x,b}=\varphi _{x},\\[1pt] \varphi _{+y,b}=-\varphi _{y}+\varphi _{\varDelta},\\[1pt] \varphi _{-y,b}=\varphi _{y}+\varphi _{\varDelta}\end{cases} \end{equation}

and

(2.7b) \begin{equation} \begin{cases} \varphi _{+x,t}=-\varphi _{x}-\varphi _{z},\\[1pt] \varphi _{-x,t}=\varphi _{x}-\varphi _{z},\\[1pt] \varphi _{+y,t}=-\varphi _{y}+\varphi _{\varDelta}-\varphi _{z},\\[1pt] \varphi _{-y,t}=\varphi _{y}+\varphi _{\varDelta} -\varphi _{z}\end{cases} \end{equation}

Under this parameterisation, the phase difference between counter-propagating waves along the x-direction on both substrates is $2\varphi _{x}$ ; the phase difference between counter-propagating waves along the y-direction is $2\varphi _{y}$ ; the phase difference between SSAWs propagating along x- and y- directions is $\varphi _{\varDelta}$ ; and the phase difference between the SSAW fields on the two substrates is $\varphi _{\textit{z}}$ . While we operate reduction on initial phases, one still can get the analytical solutions of the linear governing equations in (2.4) using the method detailed in § 2.5 if every IDT is allowed its own independent phase (eight distinct values). Such fully arbitrary phase sets are, however, rare in practical device designs and lead to considerable mathematical complexity. Therefore, they are not further pursued in the main text.

To maintain a well-defined 2-D SSAW and avoid undesirable asymmetric wave fields, excitation of the four IDTs on each substrate with equal amplitude is assumed, i.e. $2W_{j,b}=W_{b}$ and $2W_{j,t}=W_{t}$ for $j=+x,-x,+y,-y$ . The vibrations of the two substrates in (2.6) can be further condensed into (Bernard et al. Reference Bernard, Doinikov, Marmottant, Rabaud, Poulain and Thibault2017)

(2.8a) \begin{equation} \begin{cases} w_{b}(x,y,0,t)&\!\!=\,\sum _{j=+x,-x,+y,-y}w_{j,b}(x,y,0,t)\\[1pt] &\!\! =\,W_{b}\boldsymbol{\cdot }\textit{Re}(\exp (i\omega t)\boldsymbol{\cdot }(\cos (k_{\textit{SAW}}x-\varphi _{x})+\cos (k_{\textit{SAW}}y-\varphi _{y})\boldsymbol{\cdot }\exp (-i\varphi _{\varDelta}))),\\[1pt] w_{t}(x,y,H,t)&\!\!=\,\sum _{j=+x,-x,+y,-y}w_{j,t}(x,y,H,t)\\[1pt] &\!\! =\,W_{t}\boldsymbol{\cdot }\textit{Re}(\exp (i\omega t+i\varphi _{z})\boldsymbol{\cdot }(\cos (k_{\textit{SAW}}x-\varphi _{x})+\cos (k_{\textit{SAW}}y-\varphi _{y})\boldsymbol{\cdot }\exp (-i\varphi _{\varDelta}))) ,\end{cases} \end{equation}

and

(2.8b) \begin{equation} \begin{cases} u_{b}(x,y,0,t)=u_{t}(x,y,H,t)=0,\\[1pt] v_{b}(x,y,0,t)=v_{t}(x,y,H,t)=0 ,\end{cases} \end{equation}

where $u_{b}(x,y,0,t)$ and $u_{t}(x,y,H,t)$ are velocities along the x-direction; $v_{b}(x,y,0,t)$ and $v_{t}(x,y,H,t)$ are velocities along the y-direction. Although an analytical solution for the acoustic field can still be derived when the IDTs operate at unequal amplitudes, this situation is rare in practice, and the resulting field is no longer a pure standing wave. Consequently, it is not discussed in the main text.

The above assumptions, however, are not intended to be universally valid. In particular, they may become inadequate when substrate anisotropy, elliptical motion or directional attenuation is significant. For this reason, a general form of boundary conditions and the corresponding analytical solutions without the above simplifications of material properties and parameterisation of initial phases are provided in the Supplementary Material where is available at https://doi.org/10.1017/jfm.2026.11970. Since the general formulation in the Supplementary Material leads to the same conclusions for the applications considered in § 3, the simplified form of the analytical solution under the approximations is adopted in this section.

2.4. Reduction of dimensions

From (2.8), the total wave field on each substrate can be describe by the superposition of two SSAWs in the x- and y- directions, dependent on $(x,z,t)$ or $(y,z,t)$ , respectively. Since (2.4) is linear, the 3-D problem can be reduced to the sum of two 2-D solutions. By this reduction of dimension, the original 2-D velocity boundary condition (2.8a ) is decomposed to the sum of 1-D components as follows:

(2.9) \begin{equation} \begin{cases} w_{b}\left(x,y,0,t\right)=w_{x,b}\left(x,t\right)+w_{y,b}\left(y,t\right)\!,\\[3pt] w_{t}\left(x,y,H,t\right)=w_{x,t}\left(x,t\right)+w_{y,t}\left(y,t\right)\!, \end{cases} \end{equation}

where

\begin{equation*} \begin{cases} w_{x,b}(x,t)=W_{b}\boldsymbol{\cdot }\textit{Re}(\exp (i\omega t)\boldsymbol{\cdot }\cos (k_{\textit{SAW}}x-\varphi _{x})), \qquad \qquad \qquad \qquad \qquad \quad \,\kern0.5pt (2.10\textit{a})\\[3pt] w_{y,b}(y,t)=W_{b}\boldsymbol{\cdot }\textit{Re}(\exp (i\omega t-i\varphi _{\varDelta})\boldsymbol{\cdot }\cos (k_{\textit{SAW}}y-\varphi _{y})), \qquad \qquad \qquad \qquad \, \,\kern1pt (2.10\textit{b})\\[3pt] \begin{array}{c}\!\! w_{x,t}(x,t)=W_{t}\boldsymbol{\cdot }\textit{Re}(\exp (i\omega t+i\varphi _{z})\boldsymbol{\cdot }\cos (k_{\textit{SAW}}x-\varphi _{x})), \qquad \qquad \qquad \qquad \,\,\kern6pt (2.10\textit{c})\\[3pt] \!\!w_{y,t}(y,t)=W_{t}\boldsymbol{\cdot }\textit{Re}(\exp (i\omega t-i\varphi _{\varDelta}+i\varphi _{z})\boldsymbol{\cdot }\cos (k_{\textit{SAW}}y-\varphi _{y})). \qquad \qquad \qquad\kern3pt (2.10\textit{d}) \end{array} \end{cases} \end{equation*}

This decomposition yields two independent 2-D subproblems: one expressed with the independent variables $(x,z,t)$ and the other with $(y,z,t)$ . Each can be solved separately, and the full 3-D solution to (2.4) is then obtained by superposing the two independent 2-D solutions. For example, considering the $(x,z,t)$ subproblem, the boundary conditions (2.10a ) and (2.10c ) are applied, and the original 3-D governing equations reduce to the following 2-D form of (2.4):

(2.11) \begin{equation} \begin{cases} \dfrac{\partial \rho _{1}}{\partial t}+\rho _{0}\left(\dfrac{\partial u_{1}}{\partial x}+\dfrac{\partial w_{1}}{\partial z}\right)=0,\\[11pt] \rho _{0}\dfrac{\partial u_{1}}{\partial t}=-\dfrac{\partial p_{1}}{\partial x}+\left(\dfrac{4\mu }{3}+\mu _{b}\right)\dfrac{\partial ^{2}u_{1}}{\partial x^{2}}+\mu \dfrac{\partial ^{2}u_{1}}{\partial z^{2}}+\left(\dfrac{\mu }{3}+\mu _{b}\right)\dfrac{\partial ^{2}w_{1}}{\partial x\partial z},\\[11pt] \rho _{0}\dfrac{\partial w_{1}}{\partial t}=-\dfrac{\partial p_{1}}{\partial z}+\left(\dfrac{4\mu }{3}+\mu _{b}\right)\dfrac{\partial ^{2}w_{1}}{\partial z^{2}}+\mu \dfrac{\partial ^{2}w_{1}}{\partial x^{2}}+\left(\dfrac{\mu }{3}+\mu _{b}\right)\dfrac{\partial ^{2}u_{1}}{\partial x\partial z},\\[11pt] p_{1}=c^{2}\rho _{1} .\end{cases} \end{equation}

Figure 2 provides a geometric illustration of this dimensionality reduction. Taking the top substrate as an example, $w_{x,t}(x,t)$ corresponds to a SSAW that varies only along the x-direction and remains uniform along the y-direction; $w_{y,t}(y,t)$ corresponds to a SSAW that varies only along the y-direction. When these two components are superposed to form $w_{t}(x,y,H,t)$ , a 2-D SSAW field is formed. The resulting 2-D SSAW patterns are dependent on the phase difference $\varphi _{\varDelta}$ (Bernard et al. Reference Bernard, Doinikov, Marmottant, Rabaud, Poulain and Thibault2017). Yet, not every value of $\varphi _{\varDelta}$ produces a qualitatively distinct pattern. Because $\exp (0)=1$ and $\exp (i\pi /2)=i$ , any phase-shifted wave with an arbitrary $\varphi _{\varDelta}$ can be expressed as a linear combination of the two canonical cases shown in figure 2 ( $\varphi _{\varDelta}=0$ and $\varphi _{\varDelta}=\pi /2$ ). The different patterns illustrated in figure 2 are explained analytically later in § 3.2. The same property holds for the 2-D SSAW on the bottom substrate.

Figure 2.

Geometric interpretation of reducing a 2-D SAW velocity boundary into a superposition of 1-D SAWs.

Figure 2. Long description

The left surface varies periodically only along the x direction, and the middle surface varies periodically only along the y direction. Adding the two orthogonal one-dimensional standing-wave components produces a two-dimensional surface field. When the relative phase difference phi Delta is zero, the result has a broad central maximum surrounded by a connected low region and corner maxima. When phi Delta equals pi over 2, the result is a checkerboard-like array of alternating maxima and minima. The diagram illustrates the dimensional reduction used in the analysis: the full two-dimensional boundary can be represented as the sum of two independent one-dimensional components.

2.5. First-order field for viscous fluid

As discussed in § 2.3, the fluid domain is driven by boundary conditions that oscillate at a single frequency, thus the steady-state solution of (2.4) can be sought directly in the frequency domain. Furthermore, § 2.4 introduced the decomposition of the 3-D problem to the sum of 2-D subproblems. Without loss of generality, we first examine the $(x,z,t)$ -dependent 2-D subproblem

(2.12) \begin{align} \begin{cases} \boldsymbol{v}_{1}=\textit{Re}\left(\dfrac{\widetilde{\boldsymbol{V}}(z)}{2}\boldsymbol{\cdot }\left(\exp \left(i\omega t-ik_{\textit{SAW}}x+i\varphi _{x}\right)+\exp \left(i\omega t+ik_{\textit{SAW}}x-i\varphi _{x}\right)\right)\right)\!,\\[11pt] p_{1}=\textit{Re}\left(\dfrac{\tilde{P}(z)}{2}\boldsymbol{\cdot }\left(\exp \left(i\omega t-ik_{\textit{SAW}}x+i\varphi _{x}\right)+\exp \left(i\omega t+ik_{\textit{SAW}}x-i\varphi _{x}\right)\right)\right)\!,\\[11pt] \rho _{1}=\textit{Re}\left(\dfrac{\tilde{\rho }(z)}{2}\boldsymbol{\cdot }\left(\exp \left(i\omega t-ik_{\textit{SAW}}x+i\varphi _{x}\right)+\exp \left(i\omega t+ik_{\textit{SAW}}x-i\varphi _{x}\right)\right)\right)\!. \end{cases} \end{align}

Here, the tilde ‘ $\tilde{\boldsymbol{\cdot }} $ ’ denotes the complex amplitude independent of time, and $\widetilde{\boldsymbol{V}}(z)=(\tilde{U}(z),\tilde{V}(z),\tilde{W}(z))$ . For every variable defined in (2.12), we can solve separately for the components containing $\exp (i\omega t-ik_{\textit{SAW}}x+i\varphi _{x})$ and $\exp (i\omega t+ik_{\textit{SAW}}x-i\varphi _{x})$ . Any first-order quantity $Q_{1}$ , e.g. $\boldsymbol{v}_{1},p_{1},\rho _{1}$ , can therefore be decomposed as

\begin{align*} & \hspace{7.5pc} \begin{cases} Q_{1}=\textit{Re}\left(Q_{1,+}+Q_{1,-}\right)\!, & \qquad \qquad \quad (2.13\textit{a})\\[5pt] Q_{1,+}=\dfrac{\tilde{Q}(z)}{2}\boldsymbol{\cdot }\exp \left(i\omega t+ik_{\textit{SAW}}x-i\varphi _{x}\right), & \qquad \qquad \quad (2.13\textit{b})\\[5pt] Q_{1,-}=\dfrac{\tilde{Q}(z)}{2}\boldsymbol{\cdot }\exp \left(i\omega t-ik_{\textit{SAW}}x+i\varphi _{x}\right)\!, & \qquad \qquad \quad (2.13\textit{c}) \end{cases} \end{align*}

where $Q_{1,+}$ represents any first-order quantity carrying component of $\exp (i\omega t+ik_{\textit{SAW}}x-i\varphi _{x})$ , and the first-order quantities $Q_{1,+}$ and $Q_{1,-}$ are formally symmetric. We therefore first solve for the component $Q_{1,+}$ . The decomposition in (2.13) shows that $Q_{1,+}$ satisfies

(2.14) \begin{equation} \begin{cases} \dfrac{\partial Q_{1,+}}{\partial x}=ik_{\textit{SAW}}Q_{1,+},\\[11pt] \dfrac{\partial Q_{1,+}}{\partial t}=i\omega Q_{1,+} .\end{cases} \end{equation}

Substituting the decomposed $Q_{1,+}$ into (2.11), we obtain the following ordinary differential equations (ODEs) in the variable z:

\begin{align*} \begin{cases} i\omega \rho _{1,+}+\rho _{0}\left({ik_{\textit{SAW}}}u_{1,+}+\dfrac{\mathrm{d}w_{1,+}}{\mathrm{d}z}\right)=0, & \quad (2.15\textit{a})\\[9pt] i\omega \rho _{0}u_{1,+}=-ik_{\textit{SAW}}p_{1,+}-\left(\dfrac{4\mu }{3}+\mu _{b}\right){k}_{\textit{SAW}}^{2}u_{1,+}+\mu \dfrac{\mathrm{d}^{2}u_{1,+}}{\mathrm{d}z^{2}}+ik_{\textit{SAW}}\left(\dfrac{\mu }{3}+\mu _{b}\right)\dfrac{\mathrm{d}w_{1,+}}{\mathrm{d}z}, &\quad (2.15\textit{b})\\[9pt] i\omega \rho _{0}w_{1,+}=-\dfrac{\mathrm{d}p_{1,+}}{\mathrm{d}z}+\left(\dfrac{4\mu }{3}+\mu _{b}\right)\dfrac{\mathrm{d}^{2}w_{1,+}}{\mathrm{d}z^{2}}-\mu {k}_{\textit{SAW}}^{2}w_{1,+}+ik_{\textit{SAW}}\left(\dfrac{\mu }{3}+\mu _{b}\right)\dfrac{\mathrm{d}u_{1,+}}{\mathrm{d}z}, &\quad (2.15\textit{c})\\[9pt] p_{1,+}=c^{2}\rho _{1,+}. &\quad (2.15\textit{d}) \end{cases} \end{align*}

Eliminating $\rho _{1,+}$ by substituting (2.15d ) into (2.15a ) yields

(2.16) \begin{equation}i\omega p_{1,+}+\rho _{0}c^{2}\left({ik_{\textit{SAW}}}u_{1,+}+\frac{\mathrm{d}w_{1,+}}{\mathrm{d}z}\right)=0.\end{equation}

Substituting (2.16) into (2.15b ) and (2.15c ) gives

(2.17) \begin{align} \begin{cases} i\omega \rho _{0}u_{1,+}=k_{\textit{SAW}}\dfrac{\rho _{0}c^{2}}{\omega }\bigg({ik_{\textit{SAW}}}u_{1,+}+\dfrac{\mathrm{d}w_{1,+}}{\mathrm{d}z}\bigg)-\bigg(\dfrac{4\mu }{3}+\mu _{b}\bigg){k}_{\textit{SAW}}^{2}u_{1,+}+\mu \dfrac{\mathrm{d}^{2}u_{1,+}}{\mathrm{d}z^{2}}+ik_{\textit{SAW}} \bigg(\dfrac{\mu }{3}+\mu _{b}\bigg)\dfrac{\mathrm{d}w_{1,+}}{\mathrm{d}z},\\[13pt] i\omega \rho _{0}w_{1,+}=-i\dfrac{\rho _{0}c^{2}}{\omega }\bigg(ik_{\textit{SAW}}\dfrac{\mathrm{d}u_{1,+}}{\mathrm{d}z}+\dfrac{\mathrm{d}^{2}w_{1,+}}{\mathrm{d}z^{2}}\bigg)+\bigg(\dfrac{4\mu }{3}+\mu _{b}\bigg)\dfrac{\mathrm{d}^{2}w_{1,+}}{\mathrm{d}z^{2}}-\mu {k}_{\textit{SAW}}^{2}w_{1,+}+ik_{\textit{SAW}}\bigg(\dfrac{\mu }{3}+\mu _{b}\bigg)\dfrac{\mathrm{d}u_{1,+}}{\mathrm{d}z} .\end{cases} \end{align}

Equation (2.17) includes a set of coupled, second-order, linear, homogeneous ODEs. For a viscous fluid they can be solved by an eigen-mode (characteristic-mode) method. Setting $\boldsymbol{\Gamma}_{+}=[\!\begin{array}{cc} u_{1,+} & w_{1,+} \!\end{array}]^{T}$ , then (2.17) can be rewritten as

(2.18) \begin{equation}\boldsymbol{M}_{\mathbf{2}}\frac{\mathrm{d}^{2}\boldsymbol{\Gamma}_{+}}{\mathrm{d}z^{2}}+\boldsymbol{M}_{\mathbf{1}}\frac{\mathrm{d}\boldsymbol{\Gamma}_{+}}{\mathrm{d}z}+\boldsymbol{M}_{\mathbf{0}}\boldsymbol{\Gamma}_{+}=0,\end{equation}

where

(2.19) \begin{align} \begin{cases} \boldsymbol{M}_{\mathbf{2}}=\left[\begin{array}{cc} \mu & 0\\[5pt] 0 & \left(\dfrac{4\mu }{3}+\mu _{b}\right)-i\dfrac{\rho _{0}c^{2}}{\omega } \end{array}\right]\!,\\[15pt] \boldsymbol{M}_{\mathbf{1}}=\left[\begin{array}{cc} 0 & k_{\textit{SAW}}\dfrac{\rho _{0}c^{2}}{\omega }+ik_{\textit{SAW}}\left(\dfrac{\mu }{3}+\mu _{b}\right)\\[8pt] k_{\textit{SAW}}\dfrac{\rho _{0}c^{2}}{\omega }+ik_{\textit{SAW}}\left(\dfrac{\mu }{3}+\mu _{b}\right) & 0 \end{array}\right]\!,\\[15pt] \boldsymbol{M}_{\mathbf{0}}=\left[\begin{array}{cc} -i\omega \rho _{0}+{ik}_{\textit{SAW}}^{2}\dfrac{\rho _{0}c^{2}}{\omega }-\left(\dfrac{4\mu }{3}+\mu _{b}\right){k}_{\textit{SAW}}^{2} & 0\\[5pt] 0 & -i\omega \rho _{0}-\mu {k}_{\textit{SAW}}^{2} \end{array}\right]\!. \end{cases} \end{align}

Owing to the linearity of (2.17), its solution can be sought in the form

(2.20) \begin{equation}\boldsymbol{\Gamma}_{+}=\sum _{j=1,2,3,4}{C_{+,j}}\boldsymbol{E}_{+,j}\exp \left(\lambda _{+,j}z\right)\boldsymbol{\cdot }\exp \left(i\omega t+ik_{\textit{SAW}}x-i\varphi _{x}\right),\end{equation}

where the coefficients $C_{+,j}$ are constants determined by the boundary conditions; $\boldsymbol{E}_{+,j}$ is the generalised eigenvector (independent of $x,z$ and $t$ ) and $\lambda _{+,j}$ is the corresponding generalised eigenvalue.

Figure 3.

Eigenvalue spectra for liquids of different viscosities.

Figure 3. Long description

The left plot shows two eigenvalue pairs symmetric about the origin in the complex plane. Engine-oil eigenvalues are red, water eigenvalues are blue, and inviscid-water eigenvalues are green. Inviscid water has only the nearly vertical bulk-wave pair, with lambda 1 equal to 52.4 i times 10 cubed per metre. Water has a similar bulk pair, 0.004 plus 52.4 i, and a much larger viscous pair, 6428 plus 6428 i, aligned approximately 45 degrees from the real axis. Engine oil has two pairs of more comparable magnitude, 7.36 plus 80 i and 216.1 plus 217.1 i. The accompanying table lists shear viscosity, bulk viscosity, and these eigenvalues for the three fluids. Large real parts represent rapid viscous attenuation near boundaries; nearly imaginary eigenvalues represent bulk acoustic-wave propagation.

Substituting (2.20) into (2.18) yields

(2.21) \begin{equation}\left[\begin{array}{cc} {\lambda }_{+,j}^{2}\mu -\left(i\omega \rho _{0}+{k}_{\textit{SAW}}^{2}R\right) & \lambda _{+,j}ik_{\textit{SAW}}(R-\mu )\\[7pt] \lambda _{+,j}ik_{\textit{SAW}}(R-\mu ) & {\lambda }_{+,j}^{2}R-i\omega \rho _{0}-\mu {k}_{\textit{SAW}}^{2} \end{array}\right]\boldsymbol{E}_{+,j}=0,\end{equation}

where $R=((({4\mu })/{3})+\mu _{b})-i(({\rho _{0}c^{2}})/{\omega })$ .

Taking the determinant of the matrix in (2.21) and setting it to zero gives a quartic algebraic equation for $\lambda _{+,j}$

(2.22) \begin{align} \mu R{\lambda }_{+,j}^{4}&-\big(Ri\omega \rho _{0}+\mu \big(i\omega \rho _{0}+2Rk_{\textit{SAW}}^{2}\big)\big){\lambda }_{+,j}^{2} \nonumber\\[3pt] &+\big(i\omega \rho _{0}+{k}_{\textit{SAW}}^{2}R\big)\big(i\omega \rho _{0}+\mu {k}_{\textit{SAW}}^{2}\big)=0. \end{align}

Solving this quartic equation yields four roots for $\lambda _{+,j}$ , which correspond to the four generalised eigenvalues of (2.21). Since (2.22) contains only even powers of $\lambda _{+,j}$ , the eigenvalues can be written as

(2.23) \begin{align}\lambda _{+,j}=& \pm \left(2\mu R\right)^{-\frac{1}{2}}\left(i\omega \rho _{0}R+\mu \left(i\omega \rho _{0}+2R {k}_{\textit{SAW}}^{2}\right)\pm \left(\left(i\omega \rho _{0}R+\mu \left(i\omega \rho _{0}+2Rk_{\textit{SAW}}^{2}\right)\right)^{2} \right. \right. \notag \\&\quad \left. \left. -4\mu R\left(i\omega \rho _{0}+{k}_{\textit{SAW}}^{2}R\right)\left(i\omega \rho _{0}+\mu {k}_{\textit{SAW}}^{2}\right)\right)^{\frac{1}{2}}\right)^{\frac{1}{2}}.\end{align}

Figure 3 shows the generalised eigenvalues $\lambda _{+,j}$ computed from (2.23) for different liquids. For an ideally inviscid liquid, only two finite eigenvalues exist and both are purely imaginary. For low-viscosity liquids such as water, two additional eigenvalues with large magnitudes appear at approximately 45° and −135° in the complex plane, while the remaining eigenvalues cluster near the imaginary axis. For high-viscosity liquids such as engine oil, four eigenvalues of comparable magnitude emerge, symmetrically distributed across the first and third quadrants. Both the highly viscous engine oil and the less viscous water have four eigenvalues, where $\lambda _{+,3}$ and $\lambda _{+,4}$ exhibit a markedly larger real part than $\lambda _{+,1}$ and $\lambda _{+,2}$ . These real components $\lambda _{+,3}$ and $\lambda _{+,4}$ account for viscous dissipation near the boundary layer, which will be explained later. In contrast, $\lambda _{+,1}$ and $\lambda _{+,2}$ are almost purely imaginary, representing a bulk acoustic wave (BAW) outside the boundary layer. In § 3.4 we will compare how the boundary-layer thickness varies with the real parts of $\lambda _{+,3}$ and $\lambda _{+,4}$ for fluids of different viscosities.

From (2.23) we obtain the four distinct generalised eigenvalues $\lambda _{+,j}(\mu , \mu _{b}, \omega , \rho _{0}, c, k_{\textit{SAW}})$ , $j=1,2,3,4$ , of the $\boldsymbol{\Gamma}_{+}$ field. To obtain the generalised eigenvalues $\lambda _{-,j}$ for the companion field $\boldsymbol{\Gamma}_{-}$ , one can reverse the sign of every $k_{\textit{SAW}}$ appearing in the matrix of (2.21) and re-evaluate the eigenvalues of matrix in the following equation:

(2.24) \begin{equation}\left[\begin{array}{cc} {\lambda }_{-,j}^{2}\mu -\left(i\omega \rho _{0}+{k}_{\textit{SAW}}^{2}R\right) & -\lambda _{-,j}ik_{\textit{SAW}}(R-\mu )\\[5pt] -\lambda _{-,j}ik_{\textit{SAW}}(R-\mu ) & {\lambda }_{-,j}^{2}R-i\omega \rho _{0}-\mu {k}_{\textit{SAW}}^{2} \end{array}\right]\boldsymbol{E}_{-,j}=0.\end{equation}

From the structure of (2.24) and (2.21), one can find that $\lambda _{+,j}=\lambda _{-,j}$ . Henceforth, we simply use $\lambda _{j}$ for either set.

Inserting the generalised eigenvalues $\lambda _{j}$ into (2.21) and (2.24) gives the generalised eigenvectors $\boldsymbol{E}_{+,j}$ and $\boldsymbol{E}_{-,j}$ . Writing them in terms of velocity components as $\boldsymbol{E}_{+,j}=[{E}_{+,j}^{u} \ \ {E}_{+,j}^{w}]^{T}$ and $\boldsymbol{E}_{-,j}=[{E}_{-,j}^{u} \ \ {E}_{-,j}^{w}]^{T}$ obtains

(2.25a) \begin{align} \begin{cases} {E}_{+,j}^{u}=i\dfrac{{\lambda }_{j}^{2}R-i\omega \rho _{0}-\mu {k}_{\textit{SAW}}^{2}}{\lambda _{j}k_{\textit{SAW}}(R-\mu )},\\[7pt] {E}_{+,j}^{w}=1 ,\end{cases} \end{align}

and

(2.25b) \begin{align} \begin{cases} {E}_{-,j}^{u}=-{E}_{+,j}^{u},\\[5pt] {E}_{-,j}^{w}={E}_{+,j}^{w}. \end{cases} \end{align}

Because (2.4) is linear, the velocity amplitudes $W_{b}$ and $W_{t}$ prescribed by the boundary conditions (2.10) can be separated from the constants $C_{+,j}$ in (2.20)

(2.26a) \begin{equation}C_{+,j}=C_{+,b,j}\frac{W_{b}}{2}+C_{+,t,j}\frac{W_{t}\exp \left(i\varphi _{z}\right)}{2}.\end{equation}

Symmetrically,

(2.26b) \begin{equation}C_{-,j}=C_{-,b,j}\frac{W_{b}}{2}+C_{-,t,j}\frac{W_{t}\exp \left(i\varphi _{z}\right)}{2}.\end{equation}

According to the definitions of $C_{+,b,j}$ and $C_{+,t,j}$ in (2.26), the constants can be found in

(2.27a) \begin{equation}\left[\begin{array}{cccc} {E}_{+,1}^{u} & {E}_{+,2}^{u} & {E}_{+,3}^{u} & {E}_{+,4}^{u}\\[5pt] {E}_{+,1}^{w} & {E}_{+,2}^{w} & {E}_{+,3}^{w} & {E}_{+,4}^{w}\\[5pt] {E}_{+,1}^{u}\exp (\lambda _{1}H) & {E}_{+,2}^{u}\exp (\lambda _{2}H) & {E}_{+,3}^{u}\exp (\lambda _{3}H) & {E}_{+,4}^{u}\exp (\lambda _{4}H)\\[5pt] {E}_{+,1}^{w}\exp (\lambda _{1}H) & {E}_{+,2}^{w}\exp (\lambda _{2}H) & {E}_{+,3}^{w}\exp (\lambda _{3}H) & {E}_{+,4}^{w}\exp (\lambda _{4}H) \end{array}\right]\left[\begin{array}{c} C_{+,b,1}\\[5pt] C_{+,b,2}\\[5pt] C_{+,b,3}\\[5pt] C_{+,b,4} \end{array}\right]=\left[\begin{array}{c} 0\\[5pt] 1\\[5pt] 0\\[5pt] 0 \end{array}\right]\!,\end{equation}

and

(2.27b) \begin{equation}\left[\begin{array}{cccc} {E}_{+,1}^{u} & {E}_{+,2}^{u} & {E}_{+,3}^{u} & {E}_{+,4}^{u}\\[5pt] {E}_{+,1}^{w} & {E}_{+,2}^{w} & {E}_{+,3}^{w} & {E}_{+,4}^{w}\\[5pt] {E}_{+,1}^{u}\exp (\lambda _{1}H) & {E}_{+,2}^{u}\exp (\lambda _{2}H) & {E}_{+,3}^{u}\exp (\lambda _{3}H) & {E}_{+,4}^{u}\exp (\lambda _{4}H)\\[5pt] {E}_{+,1}^{w}\exp (\lambda _{1}H) & {E}_{+,2}^{w}\exp (\lambda _{2}H) & {E}_{+,3}^{w}\exp (\lambda _{3}H) & {E}_{+,4}^{w}\exp (\lambda _{4}H) \end{array}\right]\left[\begin{array}{c} C_{+,t,1}\\[5pt] C_{+,t,2}\\[5pt] C_{+,t,3}\\[5pt] C_{+,t,4} \end{array}\right]=\left[\begin{array}{c} 0\\[5pt] 0\\[5pt] 0\\[5pt] 1 \end{array}\right]\!.\end{equation}

Replacing every subscript ‘+’ with ‘–’ in (2.27) gives the analogous expressions for $C_{-,b,j}$ and $C_{-,t,j}$ that satisfy the boundary conditions. Inserting the relation between $\boldsymbol{E}_{+,j}$ and $\boldsymbol{E}_{-,j}$ from (2.25) into (2.27) yields

(2.28) \begin{align} \begin{cases} C_{+,b,j}=C_{-,b,j},\\[3pt] C_{+,t,j}=C_{-,t,j} .\end{cases} \end{align}

Consequently, we no longer distinguish between $C_{+,b,j}$ and $C_{-,b,j}$ , and denote both as $C_{b,j}$ . Likewise, $C_{+,t,j}$ and $C_{-,t,j}$ are both written as $C_{t,j}$ .

Inserting the expressions for $C_{b,j}$ and $C_{t,j}$ obtained from (2.27), together with the generalised eigenvalues in (2.23) and eigenvectors in (2.25), into (2.20) gives

(2.29a) \begin{equation}\boldsymbol{\Gamma}_{+}= \!\sum _{j=1,2,3,4} \!\!\left(C_{b,j}\frac{W_{b}}{2}+C_{t,j}\frac{W_{t}\exp \left(i\varphi _{z}\right)}{2}\right)\boldsymbol{E}_{+,j}\exp (\lambda _{j}z )\boldsymbol{\cdot }\exp (i\omega t+ik_{\textit{SAW}}x-i\varphi _{x}).\end{equation}

Symmetrically,

(2.29b) \begin{equation}\boldsymbol{\Gamma}_{-}= \! \sum _{j=1,2,3,4} \! \! \left(C_{b,j}\frac{W_{b}}{2}+C_{t,j}\frac{W_{t}\exp \left(i\varphi _{z}\right)}{2}\right)\boldsymbol{E}_{-,j}\exp (\lambda _{j}z)\boldsymbol{\cdot }\exp (i\omega t-ik_{\textit{SAW}}x+i\varphi _{x}).\end{equation}

Combining (2.29a ) and (2.29b ) for the $(x,z,t)$ subproblem and simplifying the equations yields

(2.30) \begin{align}&\left[\!\begin{array}{c} u_{1}\\ w_{1} \end{array}\!\right]=\boldsymbol{\Gamma}=\textit{Re}\big(\boldsymbol{\Gamma}_{+}+\boldsymbol{\Gamma}_{-}\big) \notag \\[3pt] &\, =\textit{Re}\Bigg(\!\exp (i\omega t)\boldsymbol{\cdot }\sum _{j=1,2,3,4}\left[\begin{array}{c} -\dfrac{{\lambda }_{j}^{2}R-i\omega \rho _{0}-\mu {k}_{\textit{SAW}}^{2}}{\lambda _{j}k_{\textit{SAW}}(R-\mu )}\sin (k_{\textit{SAW}}x-\varphi _{x})\\[11pt] \cos (k_{\textit{SAW}}x-\varphi _{x}) \end{array}\right]\big(C_{b,j}W_{b}+C_{t,j}W_{t}\exp (i\varphi _{z})\big)\exp (\lambda _{j}z)\!\Bigg).\end{align}

Applying the same procedure outlined in (2.12)–(2.30) to the $(y,t)$ -dependent boundary conditions given by (2.10b ) and (2.10d ), we obtain the corresponding $(y,z,t)$ subproblem. Superposing these 2-D solutions gives the complete 3-D first-order velocity field

(2.31) \begin{align} \begin{cases} u_{1}=\textit{Re}\Bigg(\!\!-\exp (i\omega t)\sin (k_{\textit{SAW}}x-\varphi _{x})\sum _{j=1,2,3,4}\dfrac{{\lambda }_{j}^{2}R-i\omega \rho _{0}-\mu {k}_{\textit{SAW}}^{2}}{\lambda _{j}k_{\textit{SAW}}(R-\mu )}\big(C_{b,j}W_{b}+C_{t,j}W_{t}\exp (i\varphi _{z})\big)\exp (\lambda _{j}z)\!\Bigg),\\[13pt] v_{1}=\textit{Re}\Bigg(\!\!-\exp (i\omega t-i\varphi _{\varDelta})\sin (k_{\textit{SAW}}y-\varphi _{y})\sum _{j=1,2,3,4}\dfrac{{\lambda }_{j}^{2}R-i\omega \rho _{0}-\mu {k}_{\textit{SAW}}^{2}}{\lambda _{j}k_{\textit{SAW}}(R-\mu )}\big(C_{b,j}W_{b}+C_{t,j}W_{t}\exp (i\varphi _{z})\big)\exp (\lambda _{j}z)\!\Bigg),\\[13pt] w_{1}=\textit{Re}\big(\!\exp (i\omega t)\sum _{j=1,2,3,4}\big(C_{b,j}W_{b}+C_{t,j}W_{t}\exp (i\varphi _{z})\big)\exp (\lambda _{j}z )\\[6pt] \qquad\qquad\qquad\qquad\qquad\qquad\kern-1.5pt\boldsymbol{\cdot }\,\big(\cos (k_{\textit{SAW}}x-\varphi _{x})+\exp (-i\varphi _{\varDelta})\cos (k_{\textit{SAW}}y-\varphi _{y})\big)\big), \end{cases} \end{align}

and the first-order pressure and density fields

(2.32) \begin{align} \begin{cases} p_{1}=c^{2}\rho _{1},\\[11pt] \rho _{1}=\textit{Re}(i\frac{\rho _{0}}{\omega }\exp (i\omega t) \sum _{j=1,2,3,4}\Bigg( \frac{i\omega \rho _{0}+\mu {k}_{\textit{SAW}}^{2}-{\lambda }_{j}^{2}\mu }{\lambda _{j}(R-\mu )} \Bigg)(C_{b,j}W_{b}+C_{t,j}W_{t}\exp (i\varphi _{z}))\exp (\lambda _{j}z)\\[11pt] \quad\quad\boldsymbol{\cdot }\,(\cos (k_{\textit{SAW}}x-\varphi _{x})+\exp (-i\varphi _{\varDelta})\cos (k_{\textit{SAW}}y-\varphi _{y}))) ,\end{cases} \end{align}

where $C_{b,j}$ and $C_{t,j}$ , defined in (2.27), are functions of $\mu , \mu _{b}, \omega , \rho _{0}, c, k_{\textit{SAW}},H$ . Equations (2.31) and (2.32) give the complete first-order solution for the viscous case. These first-order velocity, pressure and density fields also provide the basis for constructing the time-averaged second-order acoustic streaming field, which will be formulated in § 2.7. Before introducing the second-order extension, we first complete the discussion of the first-order problem by considering the inviscid limit in § 2.6.

2.6. First-order field for inviscid fluid

As seen in figure 3, for the case of water, the eigenvalues $\lambda _{3}$ and $\lambda _{4}$ (which represent viscous dissipation) are much larger in magnitude than $\lambda _{1}$ and $\lambda _{2}$ (which describe BAW propagation). This implies that the characteristic dissipation length in water is far shorter than the BAW wavelength. Consequently, for microfluidic applications that use water or other low-viscosity fluids, viscous dissipation can be neglected in almost the entire interested region, yielding a compact closed-form solution. Dropping viscosity by removing the highest-order coefficient in (2.21) yields

(2.33) \begin{equation}-c^{2}{\lambda }_{+,j}^{2}-\omega ^{2}+{k}_{\textit{SAW}}^{2}c^{2}=0,\end{equation}

where we can solve the two roots for $\lambda _{+}$ , as follows:

(2.34) \begin{align} \begin{cases} \lambda _{+,1}=ik_{\bot },\\[3pt] \lambda _{+,2}=-ik_{\bot }, \end{cases} \end{align}

and

(2.35) \begin{equation}k_{\bot }=\sqrt{{k}_{\textit{f}}^{2}-{k}_{\textit{SAW}}^{2}},\end{equation}

where $k_{\textit{f}}=({\omega }/{c})$ is the wavenumber in water, and the physical meaning of $k_{\bot }$ will be discussed further in § 3.3.

Similarly, (2.21) is simplified under the inviscid assumption

(2.36) \begin{equation}\left[\!\!\begin{array}{cc} -i\omega +{ik}_{\textit{SAW}}^{2}\frac{c^{2}}{\omega } & \lambda _{+}\left(k_{\textit{SAW}}\frac{c^{2}}{\omega }\right)\\[7pt] \lambda _{+}\left(k_{\textit{SAW}}\frac{c^{2}}{\omega }\right) & -i\frac{c^{2}}{\omega }{\lambda }_{+}^{2}-i\omega \end{array}\!\!\right]\boldsymbol{E}_{+,j}=0.\end{equation}

Inserting the eigenvalues obtained from (2.34) into (2.36) yields the eigenvectors

(2.37) \begin{align} \begin{cases} E_{+,1}=\big[\begin{array}{cc} k_{\textit{SAW}}/k_{\bot } & 1 \end{array}\big]^{T},\\[7pt] E_{+,2}=\big[\begin{array}{cc} -k_{\textit{SAW}}/k_{\bot } & 1 \end{array}\big]^{T} .\end{cases} \end{align}

Substituting the eigenvectors in (2.37) and the eigenvalues in (2.34) into (2.20), and imposing the boundary conditions (2.10a ) and (2.10c ) in the form of (2.11), $C_{1}$ and $C_{2}$ can be solved by

(2.38) \begin{align} \begin{cases} w_{1,+} (x,z=0,t)=\dfrac{W_{b}}{2}\boldsymbol{\cdot }\exp (i\omega t+ik_{\textit{SAW}}x-i\varphi _{x})= (C_{1}+C_{2})\exp (i\omega t+ik_{\textit{SAW}}x-i\varphi _{x}),\\[13pt] \begin{array}{ll} & w_{1,+} (x,z=H,t)=\dfrac{W_{t}\exp (i\varphi _{z})}{2}\boldsymbol{\cdot }\exp (i\omega t+ik_{\textit{SAW}}x-i\varphi _{x}) \end{array}\\[13pt] \quad\qquad\qquad\qquad\qquad\kern-2pt =\big(C_{1}\exp (H\lambda _{+,1})+C_{2}\exp (H\lambda _{+,2})\big)\exp (i\omega t+ik_{\textit{SAW}}x-i\varphi _{x}), \end{cases} \end{align}

i.e.

(2.39) \begin{align} \begin{cases} \dfrac{W_{b}}{2}=C_{1}+C_{2},\\[9pt] \begin{array}{ll} & \dfrac{W_{t}\exp \left(i\varphi _{z}\right)}{2} \end{array}=\left(C_{1}\exp \left(iHk_{\bot }\right)+C_{2}\exp \left(-iHk_{\bot }\right)\right)\!. \end{cases} \end{align}

Solving (2.39) for $C_{1}$ and $C_{2}$ gives

(2.40) \begin{align} \begin{cases} C_{1}=\dfrac{W_{t}\exp \left(i\varphi _{z}\right)-W_{b}\exp \left(-iHk_{\bot }\right)}{4i\sin \left(Hk_{\bot }\right)},\\[13pt] C_{2}=\dfrac{W_{b}}{2}-\dfrac{W_{t}\exp \left(i\varphi _{z}\right)-W_{b}\exp \left(-iHk_{\bot }\right)}{4i\sin \left(Hk_{\bot }\right)} .\end{cases} \end{align}

Following the same procedure outlined in (2.14)–(2.22) and (2.33)–(2.40), we obtain the complementary components $\boldsymbol{v}_{1,-},p_{1,-}$ and $\rho _{1,-}$ defined by the decomposition (2.12). Inserting (2.40) into (2.20) and then summing $\boldsymbol{v}_{1,+}$ and $\boldsymbol{v}_{1,-}$ according to (2.13a ) yields the first-order velocity field that satisfies the boundary conditions (2.9a ) and (2.9c )

(2.41) \begin{align} \begin{cases} u_{1}=\textit{Re}\Bigg(\!\exp (i\omega t )\dfrac{k_{\textit{SAW}}}{k_{\bot }}\Bigg(\!\dfrac{W_{t}\exp (i\varphi _{z})-W_{b}\exp (-iHk_{\bot })}{\sin (Hk_{\bot })}\cos (zk_{\bot })-iW_{b}\exp (-izk_{\bot })\!\Bigg)\sin (k_{\textit{SAW}}x-\varphi _{x})\!\Bigg),\\[13pt] w_{1}=\textit{Re}\Bigg(\!\exp (i\omega t)\Bigg(\!\dfrac{W_{t}\exp (i\varphi _{z})-W_{b}\exp (-iHk_{\bot })}{\sin (Hk_{\bot })}\sin (zk_{\bot })+W_{b}\exp (-izk_{\bot })\!\Bigg)\cos (k_{\textit{SAW}}x-\varphi _{x})\!\Bigg). \end{cases} \end{align}

Likewise, the $(y,z,t)$ subproblem described by boundary conditions in (2.10b ) and (2.10d ) can be solved in the same way. Superposing the 2-D sub-solutions gives the complete 3-D first-order velocity field

(2.42) \begin{align} \begin{cases} u_{1}=\textit{Re} \Bigg(\!\exp (i\omega t)\dfrac{k_{\textit{SAW}}}{k_{\bot }}\Bigg(\!\dfrac{W_{t}\exp (i\varphi _{z})-W_{b}\exp (-iHk_{\bot })}{\sin (Hk_{\bot })}\cos (zk_{\bot })-iW_{b}\exp (-izk_{\bot })\!\Bigg)\sin (k_{\textit{SAW}}x-\varphi _{x})\!\Bigg),\\[13pt] v_{1}=\textit{Re}\Bigg(\!\exp (i\omega t-i\varphi _{\varDelta})\dfrac{k_{\textit{SAW}}}{k_{\bot }}\Bigg(\dfrac{W_{t}\exp (i\varphi _{z})-W_{b}\exp (-iHk_{\bot })}{\sin (Hk_{\bot })}\cos (zk_{\bot })-iW_{b}\exp (-izk_{\bot })\!\Bigg)\!\Bigg)\\[13pt] \quad\quad\boldsymbol{\cdot }\,\sin (k_{\textit{SAW}}y-\varphi _{y}),\\[13pt] w_{1}=\textit{Re}(\exp (i\omega t)\Bigg(\dfrac{W_{t}\exp (i\varphi _{z})-W_{b}\exp (-iHk_{\bot })}{\sin (Hk_{\bot })}\sin (zk_{\bot })+W_{b}\exp (-izk_{\bot })\Bigg)\\[13pt] \quad\quad\boldsymbol{\cdot }\,(\cos (k_{\textit{SAW}}x-\varphi _{x})+\exp (-i\varphi _{\varDelta})\cos (k_{\textit{SAW}}y-\varphi _{y}))) .\end{cases} \end{align}

The first-order pressure and density fields are obtained by substituting (2.42) into (2.4), as follows:

(2.43) \begin{align} \begin{cases} p_{1}=c^{2}\rho _{1},\\[13pt] \rho _{1}=\textit{Re}\Big(i\frac{\rho _{0}}{\omega }\exp (i\omega t)\Big(\frac{{k}_{\textit{SAW}}^{2}}{k_{\bot }}+k_{\bot }\Big)\Big(\frac{W_{t}\exp (i\varphi _{z})-W_{b}\exp (-iHk_{\bot })}{\sin (Hk_{\bot })}\cos (zk_{\bot })-iW_{b}\exp (-izk_{\bot })\Big)\\[13pt] \quad\quad\boldsymbol{\cdot }\,(\cos (k_{\textit{SAW}}x-\varphi _{x})+\exp (-i\varphi _{\varDelta})\cos (k_{\textit{SAW}}y-\varphi _{y}))\Big). \end{cases} \end{align}

2.7. Second-order streaming field

Following the perturbation method in § 2.2, substituting (2.3) into (2.1) and collecting the second-order terms yields the governing equations for the second-order response. In the present section, we focus on the quasi-steady acoustic streaming field, whose characteristic evolution time is much longer than the acoustic period of the SAW excitation. Therefore, the second-order equations are averaged over one acoustic period, and the variables $\boldsymbol{v}_{2}$ , $p_{2}$ and $\rho _{2}$ hereafter denote the time-averaged second-order velocity, pressure and density fields, respectively. The governing equations are then obtained as

(2.44) \begin{align} \begin{cases} \boldsymbol{\nabla }\!p_{2}-\mu \boldsymbol{\nabla} ^{2}\boldsymbol{v}_{2}-\left(\mu _{b}+\dfrac{\mu }{3}\right)\boldsymbol{\nabla }\left(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}_{2}\right)=\left\langle -\rho _{1}\dfrac{\partial v_{1}}{\partial t}-\rho _{0}\left(\boldsymbol{v}_{1}\boldsymbol{\cdot }\boldsymbol{\nabla }\right)\boldsymbol{v}_{1}\right\rangle, \\[11pt] \rho _{0}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}_{2}=\left\langle -\boldsymbol{\nabla }\boldsymbol{\cdot }\left(\rho _{1}\boldsymbol{v}_{1}\right)\right\rangle \!,\\[11pt] p_{2}=c^{2}\rho _{2}, \end{cases} \end{align}

where the operator $\langle Q(t)\rangle$ denotes the time average of the variable $Q(t)$ , i.e. $\langle Q(t)\rangle := ({1}/{T}){\int }_{t}^{t+T}Q(t)\,\text{d}t$ . The boundary conditions are

(2.45) \begin{align} \begin{cases} \boldsymbol{v}_{2}\left(x,y,0\right)=-\left\langle \left(\boldsymbol{\xi }_{1}\left(x,y,0\right)\boldsymbol{\cdot }\boldsymbol{\nabla }\right)\boldsymbol{v}_{1}\left(x,y,0\right)\right\rangle \!, \\[11pt] \boldsymbol{v}_{2}\left(x,y,H\right)=-\left\langle \left(\boldsymbol{\xi }_{1}\left(x,y,H\right)\boldsymbol{\cdot }\boldsymbol{\nabla }\right)\boldsymbol{v}_{1}\left(x,y,H\right)\right\rangle \!,\end{cases} \end{align}

where $\boldsymbol{\xi }_{1}(x,y,z)=\int \boldsymbol{v}_{1}(x,y,z)\mathrm{d}t$ denotes the acoustic displacement. Following the method introduced in § 2.4, (2.44) can be decomposed into a set of 2-D subproblems. Here, we present the representative $(x,z)$ subproblem associated with the $x$ -directional first-order field:

(2.46) \begin{align} \begin{cases} \dfrac{\partial p_{2}}{\partial x}-\mu \left(\dfrac{\partial ^{2}u_{2}}{\partial x^{2}}+\dfrac{\partial ^{2}u_{2}}{\partial z^{2}}\right)-\left(\mu _{b}+\dfrac{\mu }{3}\right)\dfrac{\partial }{\partial x}\left(\dfrac{\partial u_{2}}{\partial x}+\dfrac{\partial w_{2}}{\partial z}\right)=\left\langle -\rho _{1}\dfrac{\partial u_{1}}{\partial t}-\rho _{0}\left(u_{1}\dfrac{\partial u_{1}}{\partial x}+w_{1}\dfrac{\partial u_{1}}{\partial z}\right)\right\rangle \!,\\[13pt]\dfrac{\partial p_{2}}{\partial z}-\mu \left(\dfrac{\partial ^{2}w_{2}}{\partial x^{2}}+\dfrac{\partial ^{2}w_{2}}{\partial z^{2}}\right)-\left(\mu _{b}+\dfrac{\mu }{3}\right)\dfrac{\partial }{\partial z}\left(\dfrac{\partial u_{2}}{\partial x}+\dfrac{\partial w_{2}}{\partial z}\right)=\left\langle -\rho _{1}\dfrac{\partial w_{1}}{\partial t}-\rho _{0}\left(u_{1}\dfrac{\partial w_{1}}{\partial x}+w_{1}\dfrac{\partial w_{1}}{\partial z} \! \right) \! \right\rangle \!,\\[13pt] \rho _{0}\left(\dfrac{\partial u_{2}}{\partial x}+\dfrac{\partial w_{2}}{\partial z}\right)=\left\langle -\dfrac{\partial \left(\rho _{1}u_{1}\right)}{\partial x}-\dfrac{\partial \left(\rho _{1}w_{1}\right)}{\partial z}\right\rangle \!.\end{cases} \end{align}

We first consider the inhomogeneous terms in (2.46). According to (2.31) and (2.32)

(2.47a) \begin{align} \begin{cases} u_{1}=\textit{Re}\left(\tilde{u}_{1}\right)=\dfrac{\tilde{u}_{1}+{\tilde{u}}_{1}^{*}}{2},\\[9pt] w_{1}=\textit{Re}\left(\tilde{w}_{1}\right)=\dfrac{\tilde{w}_{1}+{\tilde{w}}_{1}^{*}}{2},\\[9pt] \rho _{1}=\textit{Re}\left(\tilde{\rho }_{1}\right)=\dfrac{\tilde{\rho }_{1}+{\tilde{\rho }}_{1}^{*}}{2} ,\end{cases} \end{align}

where $Q^{*}$ denotes the complex conjugate of $Q$

(2.47b) \begin{align} \begin{cases} \tilde{u}_{1}=\exp (i\omega t)\boldsymbol{\cdot }\sin (k_{\textit{SAW}}x-\varphi _{x})\boldsymbol{\cdot }\sum _{j=1,2,3,4}{G}_{j}^{u}\exp (\lambda _{j}z),\\[7pt] \tilde{w}_{1}=\exp (i\omega t)\boldsymbol{\cdot }\cos (k_{\textit{SAW}}x-\varphi _{x})\boldsymbol{\cdot }\sum _{j=1,2,3,4}{G}_{j}^{w}\exp (\lambda _{j}z),\\[7pt] \tilde{\rho }_{1}=\exp (i\omega t)\boldsymbol{\cdot }\cos (k_{\textit{SAW}}x-\varphi _{x})\boldsymbol{\cdot }\sum _{j=1,2,3,4}{G}_{j}^{\rho }\exp (\lambda _{j}z), \end{cases} \end{align}

where ${G}_{j}^{u}$ , ${G}_{j}^{w}$ and ${G}_{j}^{\rho }$ are constants independent of $(x,y,z,t)$ . Similarly, all quantities denoted by ${G}_{j}^{q}$ in § 2.7 are constants obtained from (2.31) and (2.32) and are independent of $(x,y,z,t)$ . Noting that $\langle \exp (i\omega t)\exp (i\omega t)\rangle =0$ and $\langle \exp (i\omega t)\exp (-i\omega t)\rangle =1$ , the inhomogeneous terms in (2.46) can be reorganised into products of functions depending only on x and z, respectively, i.e.

(2.48) \begin{align} \begin{cases} \left\langle -\rho _{1}\dfrac{\partial u_{1}}{\partial t}-\rho _{0}\left(u_{1}\dfrac{\partial u_{1}}{\partial x}+w_{1}\dfrac{\partial u_{1}}{\partial z}\right)\right\rangle =\dfrac{1}{8}\sin \left(2k_{\textit{SAW}}x-2\varphi _{x}\right){F}_{1,{s}}^{\textit{inh}}(z),\\[13pt] \left\langle -\rho _{1}\dfrac{\partial w_{1}}{\partial t}-\rho _{0}\left(u_{1}\dfrac{\partial w_{1}}{\partial x}+w_{1}\dfrac{\partial w_{1}}{\partial z}\right)\right\rangle =\dfrac{1}{8}\cos \left(2k_{\textit{SAW}}x-2\varphi _{x}\right){F}_{2,\mathrm{c}}^{\textit{inh}}(z)+\dfrac{1}{8}{F}_{2,0}^{\textit{inh}}(z),\\[13pt] \left\langle -\dfrac{\partial \left(\rho _{1}u_{1}\right)}{\partial x}-\dfrac{\partial \left(\rho _{1}w_{1}\right)}{\partial z}\right\rangle =\dfrac{1}{8}\cos \left(2k_{\textit{SAW}}x-2\varphi _{x}\right){F}_{3,\mathrm{c}}^{\textit{inh}}(z)+\dfrac{1}{8}{F}_{3,0}^{\textit{inh}}(z), \end{cases} \end{align}

where ${F}_{j,\mathrm{s}}^{\textit{inh}}(z)$ , ${F}_{j,\mathrm{c}}^{\textit{inh}}(z)$ and ${F}_{j,0}^{\textit{inh}}(z)$ are real functions of z and can be written as sums of exponential functions. The subscripts ‘s’, ‘c’ and ‘0’ indicate that the corresponding x-dependent factors are sine, cosine and constant terms, respectively. Similarly, all quantities denoted by ${F}_{j}^{q}(z)$ in this section are functions of z obtained from (2.31) and (2.32). The boundary conditions in (2.45) can be written in the same form as (2.48), i.e.

(2.49) \begin{align} \begin{cases} u_{2}\left(x,0\right)=\dfrac{1}{8}\sin \left(2k_{\textit{SAW}}x-2\varphi _{x}\right){F}_{1,{s}}^{\textit{bond}}(0),\\[13pt] u_{2}\left(x,H\right)=\dfrac{1}{8}\sin \left(2k_{\textit{SAW}}x-2\varphi _{x}\right){F}_{1,{s}}^{\textit{bond}}(H),\\[13pt] w_{2}\left(x,0\right)=\dfrac{1}{8}\left(\cos \left(2k_{\textit{SAW}}x-2\varphi _{x}\right){F}_{1,\mathrm{c}}^{\textit{bond}}(0)+{F}_{1,0}^{\textit{bond}}(0)\right)\!,\\[13pt] w_{2}\left(x,H\right)=\dfrac{1}{8}\left(\cos \left(2k_{\textit{SAW}}x-2\varphi _{x}\right){F}_{1,\mathrm{c}}^{\textit{bond}}(H)+{F}_{1,0}^{\textit{bond}}(H)\right)\!. \end{cases} \end{align}

Since both the inhomogeneous terms in (2.48) and the boundary conditions in (2.49) are separable, the same method as that used in § 2.5 can be adopted, i.e.

(2.50) \begin{align} \begin{cases} \tilde{u}_{2}=\tilde{U}_{2,+}(z)\exp \left(i\left(2k_{\textit{SAW}}x-2\varphi _{x}\right)\right)+\tilde{U}_{2,-}(z)\exp \left(i\left(-2k_{\textit{SAW}}x+2\varphi _{x}\right)\right)+\tilde{U}_{2,0}, \\[7pt] \tilde{w}_{2}=\tilde{W}_{2,+}(z)\exp \left(i\left(2k_{\textit{SAW}}x-2\varphi _{x}\right)\right)+\tilde{W}_{2,-}(z)\exp \left(i\left(-2k_{\textit{SAW}}x+2\varphi _{x}\right)\right)+\tilde{W}_{2,0},\\[7pt] \tilde{p}_{2}=\tilde{P}_{2,+}(z)\exp \left(i\left(2k_{\textit{SAW}}x-2\varphi _{x}\right)\right)+\tilde{P}_{2,-}(z)\exp \left(i\left(-2k_{\textit{SAW}}x+2\varphi _{x}\right)\right)+\tilde{P}_{2,0} .\end{cases} \end{align}

Taking the component in (2.50) that contains the factor $\exp (i(2k_{\textit{SAW}}x-2\varphi _{x}))$ , the differential operators in (2.46) are applied in their frequency-domain form, with $({\partial }/{\partial x})=2ik_{\textit{SAW}}$ and $({\partial }/{\partial z})=({\mathrm{d}}/{\mathrm{d}z})$ , respectively. Therefore, (2.46) can be reduced to a set of linear inhomogeneous ODEs in z

\begin{align*} \begin{cases} 2ik_{\textit{SAW}}\tilde{P}_{2,+}(z)-\mu \dfrac{\mathrm{d}^{2}\tilde{U}_{2,+}(z)}{\mathrm{d}z^{2}}+\left(\mu _{b}+\dfrac{4\mu }{3}\right)4{k}_{\textit{SAW}}^{2}\tilde{U}_{2,+}(z)-2i\left(\mu _{b}+\dfrac{\mu }{3}\right)k_{\textit{SAW}}\dfrac{\mathrm{d}\tilde{W}_{2,+}(z)}{\mathrm{d}z}={F}_{1,+}^{\textit{inh}}(z), & \,\,\, (2.51\textit{a}) \\[11pt] \dfrac{\mathrm{d}\tilde{P}_{2,+}(z)}{\mathrm{d}z}+4\mu {k}_{\textit{SAW}}^{2}\tilde{W}_{2,+}(z)-2i\left(\mu _{b}+\dfrac{\mu }{3}\right)k_{\textit{SAW}}\dfrac{\mathrm{d}\tilde{U}_{2,+}(z)}{\mathrm{d}z}-\left(\mu _{b}+\dfrac{4\mu }{3}\right)\dfrac{\mathrm{d}^{2}\tilde{W}_{2,+}(z)}{\mathrm{d}z^{2}}={F}_{2,+}^{\textit{inh}}(z), & \,\,\, (2.51\textit{b})\\[11pt] \rho _{0}\left(2ik_{\textit{SAW}}\tilde{U}_{2,+}(z)+\dfrac{\mathrm{d}\tilde{W}_{2,+}(z)}{\mathrm{d}z}\right)={F}_{3,+}^{\textit{inh}}(z), & \,\,\, (2.51\textit{c}) \end{cases} \end{align*}

where

(2.51d) \begin{align} \begin{cases} {F}_{1,+}^{\textit{inh}}(z)=-\dfrac{i}{16}\sum _{\begin{array}{c} m=1,2,3,4\\ n=1,2,3,4 \end{array}}{G}_{mn,1}^{\textit{inh}}\exp \left(\left({\lambda }_{m}^{*}+\lambda _{n}\right)z\right)\!,\\[11pt] {F}_{2,+}^{\textit{inh}}(z)=\dfrac{1}{16}\sum _{\begin{array}{c} m=1,2,3,4\\ n=1,2,3,4 \end{array}}{G}_{mn,2}^{\textit{inh}}\exp \left(\left({\lambda }_{m}^{*}+\lambda _{n}\right)z\right)\!,\\[11pt] {F}_{3,+}^{\textit{inh}}(z)=\dfrac{1}{16}\sum _{\begin{array}{c} m=1,2,3,4\\ n=1,2,3,4 \end{array}}{G}_{mn,3}^{\textit{inh}}\exp \left(\left({\lambda }_{m}^{*}+\lambda _{n}\right)z\right)\!. \end{cases} \end{align}

The corresponding boundary conditions are

(2.52) \begin{align} \begin{cases} \tilde{U}_{2,+}(0)=-\dfrac{i}{16}{F}_{1,{s}}^{\textit{bond}}(0),\\[11pt] \tilde{U}_{2,+}(H)=-\dfrac{i}{16}{F}_{1,{s}}^{\textit{bond}}(H), \\[11pt]\tilde{W}_{2,+}(0)=\dfrac{1}{16}{F}_{1,\mathrm{c}}^{\textit{bond}}(0),\\[11pt] \tilde{W}_{2,+}(H)=\dfrac{1}{16}{F}_{1,\mathrm{c}}^{\textit{bond}}(H). \end{cases} \end{align}

From (2.51c ), we obtain

(2.53) \begin{equation}\tilde{U}_{2,+}(z)=\frac{1}{2ik_{\textit{SAW}}}\left(\frac{{F}_{3,+}^{\textit{inh}}(z)}{\rho _{0}}-\frac{\mathrm{d}\tilde{W}_{2,+}(z)}{\mathrm{d}z}\right)\!. \end{equation}

Substituting (2.53) into (2.51a ) gives

(2.54) \begin{align}\tilde{P}_{2,+}(z)&=\frac{\mu }{4{k}_{\textit{SAW}}^{2}}\frac{d^{3}\tilde{W}_{2,+}(z)}{dz^{3}}-\mu \frac{\mathrm{d}\tilde{W}_{2,+}(z)}{\mathrm{d}z}+\frac{{F}_{1,+}^{\textit{inh}}(z)}{2ik_{\textit{SAW}}} \notag\\ &\quad -\frac{\mu }{4\rho _{0}{k}_{\textit{SAW}}^{2}}\frac{d^{2}{F}_{3,+}^{inh}(z)}{dz^{2}}+\left(\mu _{b}+\frac{4\mu }{3}\right)\frac{{F}_{3,+}^{inh}(z)}{\rho _{0}}. \end{align}

Substituting (2.53) and (2.54) into (2.51b ) gives

(2.55a) \begin{equation}\left(\frac{\mathrm{d}^{2}}{\mathrm{d}z^{2}}-4{k}_{\textit{SAW}}^{2}\right)^{2}\tilde{W}_{2,+}(z)=S(z),\end{equation}

where the inhomogeneous term is

(2.55b) \begin{equation}S(z)=\frac{2ik_{\textit{SAW}}}{\mu }\frac{\mathrm{d}{F}_{1,+}^{inh}(z)}{\mathrm{d}z}+\frac{4{k}_{\textit{SAW}}^{2}}{\mu }{F}_{2,+}^{\textit{inh}}(z)+\frac{1}{\rho _{0}}\frac{\mathrm{d}^{3}{F}_{3,+}^{inh}(z)}{\mathrm{d}z^{3}}-\frac{4{k}_{\textit{SAW}}^{2}}{\rho _{0}}\frac{\mathrm{d}{F}_{3,+}^{inh}(z)}{\mathrm{d}z}.\end{equation}

The characteristic equation of the homogeneous part of (2.55a ) has two double roots. Thus, the homogeneous solution is

(2.56) \begin{equation}{\tilde{W}}_{2,+}^{\textit{hom}}(z)=\big({C}_{+,1}^{\textit{hom}}+z{C}_{+,2}^{\textit{hom}}\big)\exp \left(2k_{\textit{SAW}}z\right)+\big({C}_{+,3}^{\textit{hom}}+z{C}_{+,4}^{\textit{hom}}\big)\exp \left(-2k_{\textit{SAW}}z\right).\end{equation}

According to (2.51d ), the inhomogeneous term in (2.55b ) can be written as

(2.57a) \begin{equation}S(z)=\sum _{\begin{array}{c} m=1,2,3,4\\ n=1,2,3,4 \end{array}}S_{mn}\exp (s_{mn}z),\end{equation}

where

(2.57b) \begin{align} \begin{cases} S_{mn}=\dfrac{1}{16}\left(\dfrac{2k_{\textit{SAW}}s_{mn}}{\mu }{G}_{mn,1}^{\textit{inh}}+\dfrac{4{k}_{\textit{SAW}}^{2}}{\mu }{G}_{mn,2}^{\textit{inh}}+\dfrac{{s}_{mn}^{3}-4{k}_{\textit{SAW}}^{2}s_{mn}}{\rho _{0}}{G}_{mn,3}^{\textit{inh}}\right)\\[11pt] s_{mn}={\lambda }_{m}^{*}+\lambda _{n}. \end{cases} \end{align}

The corresponding particular solution for inhomogeneous part is

(2.58) \begin{equation}{\tilde{W}}_{2,+}^{\textit{inh}}(z)=\sum _{\begin{array}{c} m=1,2,3,4\\ n=1,2,3,4 \end{array}}\frac{S_{mn}}{\left({s}_{mn}^{2}-4{k}_{\textit{SAW}}^{2}\right)^{2}}\exp \left(s_{mn}z\right).\end{equation}

Therefore, the general solution of $\tilde{W}_{2,+}(z)$ is

(2.59) \begin{align}\tilde{W}_{2,+}(z)&={\tilde{W}}_{2,+}^{\textit{hom}}(z)+{\tilde{W}}_{2,+}^{\textit{inh}}(z)=\big({C}_{+,1}^{\textit{hom}}+z{C}_{+,2}^{\textit{hom}}\big)\exp \left(2k_{\textit{SAW}}z\right) \notag \\& +\big({C}_{+,3}^{\textit{hom}}+z{C}_{+,4}^{\textit{hom}}\big)\exp (-2k_{\textit{SAW}}z)+ \!\! \sum _{\begin{array}{c} m=1,2,3,4\\ n=1,2,3,4 \end{array}}\frac{S_{mn}}{\left({s}_{mn}^{2}-4{k}_{\textit{SAW}}^{2}\right)^{2}}\exp (s_{mn}z ).\end{align}

From (2.53) and (2.54), $\tilde{U}_{2,+}(z)$ and $\tilde{P}_{2,+}(z)$ can be obtained. To determine the homogeneous coefficients ${C}_{+,j}^{\textit{hom}}$ , substituting (2.53) into the boundary conditions in (2.52) yields

(2.60) \begin{align} \begin{cases} \dfrac{\mathrm{d}\tilde{W}_{2,+}(0)}{\mathrm{d}z}=\dfrac{{F}_{3,+}^{\textit{inh}}(0)}{\rho _{0}}-\dfrac{k_{\textit{SAW}}}{8}{F}_{1,{s}}^{\textit{bond}}(0),\\[13pt] \dfrac{\mathrm{d}\tilde{W}_{2,+}(H)}{\mathrm{d}z}=\dfrac{{F}_{3,+}^{\textit{inh}}(H)}{\rho _{0}}-\dfrac{k_{\textit{SAW}}}{8}{F}_{1,{s}}^{\textit{bond}}(H),\\[13pt] \tilde{W}_{2,+}(0)=\dfrac{1}{16}{F}_{1,\mathrm{c}}^{\textit{bond}}(0),\\[13pt]\tilde{W}_{2,+}(H)=\dfrac{1}{16}{F}_{1,\mathrm{c}}^{\textit{bond}}(H). \end{cases} \end{align}

Substituting (2.59) into (2.60) gives

(2.61a) \begin{equation}\boldsymbol{M}_{\text{bond}}\left[\!\begin{array}{c} {C}_{+,1}^{\textit{hom}}\\[11pt] {C}_{+,2}^{\textit{hom}}\\[11pt] {C}_{+,3}^{\textit{hom}}\\[11pt] {C}_{+,4}^{\textit{hom}} \end{array}\!\right]=\left[\!\begin{array}{c} \dfrac{{F}_{3,+}^{inh}(0)}{\rho _{0}}-\dfrac{k_{\textit{SAW}}}{8}{F}_{1,s}^{bond}(0)-\sum _{\begin{array}{c} m=1,2,3,4\\ n=1,2,3,4 \end{array}}\dfrac{s_{mn}S_{mn}}{\left({s}_{mn}^{2}-4{k}_{\textit{SAW}}^{2}\right)^{2}}\\[11pt] \dfrac{{F}_{3,+}^{inh}(H)}{\rho _{0}}-\dfrac{k_{\textit{SAW}}}{8}{F}_{1,s}^{bond}(H)-\sum _{\begin{array}{c} m=1,2,3,4\\ n=1,2,3,4 \end{array}}\dfrac{s_{mn}S_{mn}}{\left({s}_{mn}^{2}-4{k}_{\textit{SAW}}^{2}\right)^{2}}\exp \left(s_{mn}H\right)\\[11pt] \dfrac{1}{16}{F}_{1,\mathrm{c}}^{\textit{bond}}(0)-\sum _{\begin{array}{c} m=1,2,3,4\\ n=1,2,3,4 \end{array}}\dfrac{S_{mn}}{\left({s}_{mn}^{2}-4{k}_{\textit{SAW}}^{2}\right)^{2}}\\[11pt] \dfrac{1}{16}{F}_{1,\mathrm{c}}^{\textit{bond}}(H)-\sum _{\begin{array}{c} m=1,2,3,4\\ n=1,2,3,4 \end{array}}\dfrac{S_{mn}}{\left({s}_{mn}^{2}-4{k}_{\textit{SAW}}^{2}\right)^{2}}\exp \left(s_{mn}H\right) \end{array}\!\right]\!,\end{equation}

where

(2.61b) \begin{equation} \boldsymbol{M}_{\textit{bond}}= \left[\!\begin{array}{cccc} 2k_{\textit{SAW}} & 1 & -2k_{\textit{SAW}} & 1\\[9pt] 2k_{\textit{SAW}}\exp \left(2k_{\textit{SAW}}H\right) & \left(1+2Hk_{\textit{SAW}}\right)\exp \left(2k_{\textit{SAW}}H\right) & -2k_{\textit{SAW}}\exp \left(-2k_{\textit{SAW}}H\right) & \left(1-2Hk_{\textit{SAW}}\right)\exp \left(-2k_{\textit{SAW}}H\right)\\[9pt] 1 & 0 & 1 & 0\\[9pt] \exp \left(2k_{\textit{SAW}}H\right) & H\exp \left(2k_{\textit{SAW}}H\right) & \exp \left(-2k_{\textit{SAW}}H\right) & H\exp \left(-2k_{\textit{SAW}}H\right) \end{array}\!\right]\!. \end{equation}

For the component in (2.50) containing the factor $\exp (-i(2k_{\textit{SAW}}x-2\varphi _{x}))$ , repeating (2.51)–(2.61) gives

(2.62) \begin{align} \begin{cases} \tilde{U}_{2,+}(z)=-\tilde{U}_{2,-}(z),\\[9pt] \tilde{W}_{2,+}(z)=\tilde{W}_{2,-}(z) \!. \end{cases} \end{align}

Similarly, for the spatially uniform components $\tilde{U}_{2,0}(z)$ , $\tilde{W}_{2,0}(z)$ and $\tilde{P}_{2,0}(z)$ , the operators become $({\partial }/{\partial x})=0$ and $({\partial }/{\partial z})=({\mathrm{d}}/{\mathrm{d}z})$ . Therefore, the following linear inhomogeneous ODEs in the z-direction are obtained:

\begin{align*} &\hspace{7pc} \begin{cases} -\mu \dfrac{\mathrm{d}^{2}\tilde{U}_{2,0}(z)}{\mathrm{d}z^{2}}=0, &\qquad \qquad \, (2.63\textit{a})\\[9pt] \dfrac{\mathrm{d}\tilde{P}_{2,0}(z)}{\mathrm{d}z}-\left(\mu _{b}+\dfrac{4\mu }{3}\right)\dfrac{\mathrm{d}^{2}\tilde{W}_{2,0}(z)}{\mathrm{d}z^{2}}={F}_{2,0}^{\textit{inh}}(z), &\qquad \qquad \, (2.63\textit{b})\\[9pt] \rho _{0}\dfrac{\mathrm{d}\tilde{W}_{2,0}(z)}{\mathrm{d}z}={F}_{3,0}^{\textit{inh}}(z), &\qquad \qquad \, (2.63\textit{c}) \end{cases} \end{align*}

where

(2.63d) \begin{align} \begin{cases}\!\!\! & {F}_{2,0}^{\textit{inh}}(z)=\dfrac{1}{8}\sum _{\begin{array}{c} m=1,2,3,4\\ n=1,2,3,4 \end{array}}{G}_{mn,2,0}^{\textit{inh}}\exp \left(\left({\lambda }_{m}^{*}+\lambda _{n}\right)z\right)\\ & {F}_{3,0}^{\textit{inh}}(z)=\dfrac{1}{8}\sum _{\begin{array}{c} m=1,2,3,4\\ n=1,2,3,4 \end{array}}{G}_{mn,3,0}^{\textit{inh}}\exp \left(\left({\lambda }_{m}^{*}+\lambda _{n}\right)z\right)\!. \end{cases} \end{align}

Solving (2.63) with the corresponding boundary conditions gives

(2.64) \begin{align} \begin{cases} & \tilde{U}_{2,0}(z)=0\\ & \tilde{W}_{2,0}(z)=\dfrac{1}{8}{F}_{1,0}^{\textit{bond}}(0)+\dfrac{1}{8}\sum _{\begin{array}{c} m=1,2,3,4\\ n=1,2,3,4 \end{array}}\dfrac{{G}_{mn,3,0}^{\textit{inh}}}{\rho _{0}\left({\lambda }_{m}^{*}+\lambda _{n}\right)}\left(\exp \left(\left({\lambda }_{m}^{*}+\lambda _{n}\right)z\right)-1\right)\!. \end{cases} \end{align}

Combining the results in (2.50) with (2.59), (2.62) and (2.64), we obtain the second-order solution for the $(x, z)$ subproblem

(2.65a) \begin{align}\tilde{w}_{2}&=2\left(\cos \left(2k_{\textit{SAW}}x-2\varphi _{x}\right)\right)\Big(\Big({C}_{+,1}^{\textit{hom}} +z{C}_{+,2}^{\textit{hom}}\Big) \exp \left(2k_{\textit{SAW}}z\right) \notag \\&\quad +\left({C}_{+,3}^{\textit{hom}}+z{C}_{+,4}^{\textit{hom}}\right)\exp \left(-2k_{\textit{SAW}}z\right) \notag \\&\quad +\sum _{\begin{array}{c} m=1,2,3,4\\ n=1,2,3,4 \end{array}}\dfrac{S_{mn}}{\left({s}_{mn}^{2}-4{k}_{\textit{SAW}}^{2}\right)^{2}}\exp \left(s_{mn}z\right)\Big)+\dfrac{1}{8}{F}_{1,0}^{\textit{bond}}(0) \notag \\&\quad +\dfrac{1}{8}\sum _{\begin{array}{c} m=1,2,3,4\\ n=1,2,3,4 \end{array}}\frac{{G}_{mn,3,0}^{\textit{inh}}}{\left({\lambda }_{m}^{*}+\lambda _{n}\right)\rho _{0}}\left(\exp \left(\left({\lambda }_{m}^{*}+\lambda _{n}\right)z\right)-1\right)\!,\end{align}

and

(2.65b) \begin{align}\tilde{u}_{2}&=\dfrac{\sin \left(2k_{\textit{SAW}}x-2\varphi _{x}\right)}{k_{\textit{SAW}}}\Bigg(\dfrac{{F}_{3,+}^{\textit{inh}}(z)}{\rho _{0}}-\Big(2k_{\textit{SAW}}{C}_{+,1}^{\textit{hom}}+{C}_{+,2}^{\textit{hom}} +2zk_{\textit{SAW}}{C}_{+,2}^{\textit{hom}}\Big) \notag \\&\quad\exp \left(2k_{\textit{SAW}}z\right)-\Big(\!-2k_{\textit{SAW}}{C}_{+,3}^{\textit{hom}}+{C}_{+,4}^{\textit{hom}} -2zk_{\textit{SAW}}{C}_{+,4}^{\textit{hom}}\Big)\exp \left(-2k_{\textit{SAW}}z\right) \notag \\&\quad-\sum _{\begin{array}{c} m=1,2,3,4\\ n=1,2,3,4 \end{array}}\frac{s_{mn}S_{mn}}{\left({s}_{mn}^{2}-4{k}_{\textit{SAW}}^{2}\right)^{2}}\exp \left(s_{mn}z\right)\!\Bigg).\end{align}

The physical time-averaged second-order velocity field is obtained by taking the real part of (2.65).

3.1. First-order velocity field

In § 3.1, the computational domain is defined as a cubic fluid region with a square base length L of $\lambda _{\textit{SAW}}=300 \,{\unicode[Arial]{x03BC}}\mathrm{m}$ and a height H of 50 $\,{\unicode[Arial]{x03BC}}\mathrm{m}$ , resulting in 62 450 mesh nodes in the region of interest. The acoustic velocity and pressure fields within the fluid are simulated using thermoviscous acoustics, a frequency-domain interface and the total number of degrees of freedom (DOFs) in the 3-D frequency-domain simulations is 484,315. A higher DOF generally indicates increased computational cost. To ensure generality, the acoustic field depicted in figure 4 was chosen with parameters $({W_{b}}/{\omega })={}1$ nm, $({W_{t}}/{\omega })=-0.7$ nm, $\varphi _{x}=({\pi }/{2})+1$ and $\varphi _{y}=({\pi }/{2})$ . Applying these parameters, figures 4(a) and 4(b) show the structure of representative acoustofluidic fields with $\varphi _{\varDelta }=\pi /2$ and $\varphi _{\varDelta }=0$ , respectively.

Figure 4.

Characteristic first-order acoustic velocity fields for (a) $\varphi _{\varDelta }=\pi /2$ and (b) $\varphi _{\varDelta }=0$ . Comparison between theory and finite-element simulation is conducted along three characteristic cutlines in the principal coordinate directions. Red crosses mark the simulated values at computational grid nodes along these cutlines, and the black solid lines represent the analytical prediction.

Figure 4. Long description

Panel (a), for phi Delta equal to pi over 2, shows a periodic velocity-magnitude field with alternating peaks and valleys in both in-plane directions. Panel (b), for phi Delta equal to zero, shows a different pattern dominated by a central high-velocity region and larger in-plane cells. In each panel, the central three-dimensional colour map is intersected by x-, y-, and z-directed cutlines. The surrounding graphs compare analytical predictions, drawn as black solid curves, with finite-element values, shown as red crosses. The curves and markers nearly overlap along every cutline, with coefficients of determination above 0.99. The displayed velocity ranges reach approximately 0.12 metres per second in panel (a) and 0.17 metres per second in panel (b).

Due to the low viscosity of liquid water, the viscous effect is limited to a very narrow region within the boundary layers, which is detailed in § 3.4. The viscous effect outside of the boundary layers, however, is negligible. Therefore, in most region, the theoretical calculations using (2.31) closely matched the results obtained by the inviscid solution in (2.42). Therefore, we hereby employ the inviscid formula (2.42) to calculate the acoustic field within water, and we present the analytical results of the z-direction velocity component $w_{1}$ in figure 4. To validate the theoretical results from our 3-D numerical model, we specifically calculate the magnitude of $w_{1}$ along three cutlines in the x-, y- and z- directions, as shown in figure 4. The black lines in figure 4 represent the theoretically calculated $w_{1}$ magnitude along these intersecting cutlines, while the red cross markers denote the simulation data at each node position.

To assess the agreement between theoretical and numerical results, we employed the coefficient of determination $R^{2}$ , which quantifies the degree to which a theoretical framework accounts for the variability in the numerical data. An $R^{2}$ value of $1$ indicates a perfect fit between the theoretical model and the numerical data, while $R^{2}=0$ implies that the theoretical model performs no better than simply predicting the mean of the numerical data. The coefficient $R^{2}$ is computed using the Numpy 2.0.0 lib in Python 3.10.15. For the two characteristic acoustic fields in figure 4, the $R^{2}$ value consistently exceeds 0.99. The two cases with lowest $R^{2}$ values appear in the profiles of $w_{1}$ along the z-axis, where $R^{2}=0.998$ . The slight discrepancy observed near the minimum $w_{1}$ magnitude at $z=40 \,{\unicode[Arial]{x03BC}}\mathrm{m}$ is attributed to the relatively coarse computational mesh resolution in this region, which causes the abrupt changes in the $w_{1}$ magnitude not captured as accurately as in other cases. Nevertheless, an assessment of $R^{2}=0.998$ still indicates high reliability of the theoretical framework. All other comparisons along the x- and y- directions yield $R^{2}$ values greater than 0.999, demonstrating that the theoretical calculations accurately represent the numerical simulation results and significantly reduce computational cost for similar scenarios. Additional comparisons for the pressure field, the other velocity components and the viscous analytical formulation against thermoviscous simulations are provided in Appendix A, where similarly good agreement is observed.

It is noteworthy that the velocity distributions shown in figure 4 exhibit similar profiles along the z direction yet display distinct periodicities along the x and y directions. For a phase difference $\varphi _{\varDelta }=\pi /2$ (figure 4a ), the velocity field within the square domain of side length $L=\lambda _{\textit{SAW}}$ exhibits two pairs of peaks and valleys along both the x and y axes, indicating a spatial periodicity of $\lambda _{\textit{SAW}}/2$ in these directions. In contrast, for $\varphi _{\varDelta }=0$ (figure 4b ), the velocity distribution shows only one peak–trough pair along x axis and two peaks of different magnitudes along the y axis, corresponding to a periodicity of $\lambda _{\textit{SAW}}$ along both axes. These contrasting acoustic field patterns and the associated arrangements of PNs enable distinct microfluidic functionalities (Ding et al. Reference Ding, Lin, Kiraly, Yue, Li, Chiang, Shi, Benkovic and Huang2012; Melde et al. Reference Melde, Mark, Qiu and Fischer2016). Accordingly, § 3.2 presents a quantitative and qualitative analysis of the acoustic field shape and nodal positions for $\varphi _{\varDelta }=\pi /2$ and $\varphi _{\varDelta }=0$ .

3.2. Geometry of nodal surface

In typical applications of acoustofluidic particle manipulation, particles suspended in the liquid migrate towards either PNs or ANs based on their acoustic contrast factor (Baudoin & Thomas Reference Baudoin and Thomas2020). In this section, we analyse the spatial locations of PNs within the 3-D standing-wave field using the inviscid analytical solutions given in (2.43). To determine the position of PNs, we identify the locations where the complex amplitude of the first-order pressure field vanishes in (2.43), which yields

(3.1) \begin{align} \begin{cases} 0=f_{1}\left(x,y\right)\boldsymbol{\cdot }f_{2}(z),\\[9pt] f_{1}\left(x,y\right)=\cos \left(k_{\textit{SAW}}x-\varphi _{x}\right)+\exp \left(-i\varphi _{\varDelta}\right)\cos \left(k_{\textit{SAW}}y-\varphi _{y}\right)\!,\\[9pt] f_{2}(z)=\dfrac{W_{t}\exp \left(i\varphi _{z}\right)-W_{b}\exp \left(-iHk_{\bot }\right)}{\sin \left(Hk_{\bot }\right)}\cos \left(zk_{\bot }\right)-iW_{b}\exp \left(-izk_{\bot }\right)\!. \end{cases} \end{align}

Equation (3.1) describes the pressure nodal surfaces within the 3-D acoustic field. The equation factorises into two independent components: a lateral term $f_{1}(x,y)$ depending only on the $x$ - and $y$ -coordinates, and a vertical term $f_{2}(z)$ depending only on the $z$ -coordinate. This decomposition indicates that the full 3-D pressure nodal surface can be understood as the combination of two independent nodal surface systems. Setting $f_{1}(x,y)$ or $f_{2}(z)$ in (3.1) to zero yields the conditions for the two pressure nodal surface systems. Specifically, the condition $f_{1}(x,y)=0$ yields pressure nodal surfaces that are parallel to the z-axis, while the surface defined by $f_{2}(z)=0$ forms planar pressure nodal surfaces perpendicular to the z-axis. Considering the distinctive behaviour of the pressure nodal surface in the z-direction for a dual-side-actuated acoustofluidic system, a detailed analysis of the planar nodal positions defined by $f_{2}(z)=0$ will be presented in § 3.3.

The pressure nodal surface system parallel to the z-axis exhibits a more complex structure. When $\varphi _{\varDelta }=\pi /2$ , the PN condition $f_{1}(x,y)=0$ becomes

(3.2) \begin{align} \begin{array}{cc} & \\ & \cos \left(k_{\textit{SAW}}x-\varphi _{x}\right)+i\cos \left(k_{\textit{SAW}}y-\varphi _{y}\right)=0,\\ & \end{array} \end{align}

which holds only when $\cos (k_{\textit{SAW}}x-\varphi _{x})=0$ and $\cos (k_{\textit{SAW}}y-\varphi _{y})=0$ . This yields the family of solutions

(3.3) \begin{align} \begin{cases} x=\dfrac{\left(n+\dfrac{1}{2}\right)\pi +\varphi _{x}}{k_{\textit{SAW}}},\\[11pt] y=\dfrac{\left(m+\dfrac{1}{2}\right)\pi +\varphi _{y}}{k_{\textit{SAW}}}, \end{cases} n,m\in \mathbb{Z}, \end{align}

indicating that, for $\varphi _{\varDelta }=\pi /2$ , the PNs form a lattice of columns parallel to the z-axis, each separated by a distance of $\pi /k_{\textit{SAW}}=\lambda _{\textit{SAW}}/2$ . This conclusion is consistent with the results in figure 5(a), where the three columns from left to right show the low-pressure isosurface computed by (2.43), the pressure distribution obtained from numerical simulations and the experimentally observed flow field, respectively. As depicted in the left column of figure 5(a), these PN columns are denoted by the orange cylinders, while planar pressure nodal surfaces can also be observed perpendicular to the z-direction. The simulations likewise reveal that the nodal surfaces parallel to the z-axis form uniformly spaced columns with a periodicity of $\lambda _{\textit{SAW}}/2$ . Although the experimental image shows slightly distorted and vague nodal patterns near the sidewalls, the PNs at the central region of the fluid cavity retains equidistant, and the overall acoustic field preserves the same lattice structure as in analytical and numerical results. Notably, because this acoustic field configuration does not produce continuous nodal lines along both the $x$ - and $y$ -directions, it is particularly well suited for in-plane acoustofluidic particle manipulation, where particles can be rapidly and reliably trapped at these nodal columns to enable precise in-plane (x- and y- directional) control (Shen et al. Reference Shen2024).

Figure 5.

Characteristic first-order pressure fields for (a) $\varphi _{\varDelta }=\pi /2$ and (b) $\varphi _{\varDelta }=0$ . Each column shows, from left to right: the analytical pressure isosurface ( $| p_{1}| =0.1| p_{1}| _{max }$ ), the numerical pressure field and the top-view phase-contrast micrograph of the corresponding flow field reproduced from Bernard et al. (Reference Bernard, Doinikov, Marmottant, Rabaud, Poulain and Thibault2017). The analytical and computational domains have a lateral dimension of $L=2\lambda _{\textit{SAW}}$ and a height of $H=\lambda _{\textit{SAW}}/3$ . In Bernard et al. (Reference Bernard, Doinikov, Marmottant, Rabaud, Poulain and Thibault2017), the experimental set-up uses $L=12\lambda _{\textit{SAW}}$ . Yellow cylinders and dots highlight the vertically oriented pressure nodal columns parallel to the z-axis. Semi-transparent yellow planes and solid grid lines indicate the nodal planes parallel to the z axis.

Figure 5. Long description

Two rows compare phase differences phi Delta equal to pi over 2 and zero. In each row, the left image is the analytical low-pressure isosurface, the centre image is the numerical pressure-magnitude field, and the right image is a top-view phase-contrast micrograph reproduced from Bernard et al. For phi Delta equal to pi over 2, the low-pressure surfaces form a regular lattice of vertical pressure-node pillars parallel to the z axis, with spacing equal to half a surface-acoustic-wave wavelength; the numerical and experimental fields show the corresponding dense lattice. For phi Delta equal to zero, the analytical field contains extended pressure-node planes and a larger periodic grid. The numerical field and experimental micrograph show the same coarser, diagonally oriented lattice pattern.

When $\varphi _{\varDelta }=0$ , (3.1) generates a different family of pressure nodal surfaces parallel to the z-axis, expressed as

(3.4) \begin{align} \cos \left(k_{\textit{SAW}}x-\varphi _{x}\right)+\cos \left(k_{\textit{SAW}}y-\varphi _{y}\right)=0 . \end{align}

The sufficient and necessary condition for (3.4) to hold is

(3.5) \begin{align} \pm \left(y-\dfrac{\varphi _{y}}{k_{\textit{SAW}}}\right)=x-\dfrac{\varphi _{x}}{k_{\textit{SAW}}}+\dfrac{\left(2 n+1\right)\pi }{k_{\textit{SAW}}}, \quad n\in \mathbb{Z}. \end{align}

As shown in figure 5(b), for $\varphi _{\varDelta }=0$ , the resulting pressure field features two orthogonal sets of PN planes parallel to the z-axis, and each is oriented at an angle of 45 degrees to the x- and y- axes. The same pair of orthogonally intersecting surface systems are also observed in both simulations and experiments. Since the pressure gradient vanishes within each pressure nodal surface, for the applications of acoustofluidic particle manipulation, particles under this acoustic field are attracted to these nodal surfaces and cease to move. Therefore, precise trapping of particles at a specific point within the XY plane becomes challenging, and this acoustic field is generally avoided for in-plane (x- and y- directional) particle manipulation.

Equations (3.3) and (3.5) further demonstrate that the in-plane position $(x,y)$ of the PNs can be independently tuned by adjusting $(\varphi _{x},\varphi _{y})$ . Specifically, small modifications $({\Delta} \varphi _{x},{\Delta} \varphi _{y})$ to the IDT actuation phases lead to corresponding shifts in nodal position given by $({\Delta} x,{\Delta} y)=({\Delta} \varphi _{x}/k_{\textit{SAW}},{\Delta} \varphi _{y}/k_{\textit{SAW}} )$ . This provides the theoretical basis for x- and y- directional acoustic manipulation, as demonstrated in Shen et al. (Reference Shen2024). Notably, since only $f_{1}(x,y)$ in (3.1) depends on $(\varphi _{x},\varphi _{y})$ , tuning the XY-plane nodal positions via phase $(\varphi _{x},\varphi _{y})$ adjustments does not affect the nodal surfaces in the $z$ -direction. Conversely, adjustments to the z-axis nodal-plane positions based on $f_{2}(z)$ do not affect the PN positions in the XY plane. Consequently, the pressure nodal surfaces in this acoustofluidic system can be independently controlled along all three spatial directions. In practical device design, this controllability is most reliably realised in the central working region of a sufficiently large cavity, where the influence of the lateral boundaries is weak and the periodic approximation remains valid.

It should be emphasised that the comparison with the experimental data reproduced from Bernard et al. (Reference Bernard, Doinikov, Marmottant, Rabaud, Poulain and Thibault2017) is intended to validate the predicted nodal structures on the XY plane only. Since the reported experimental observations are obtained from the top-view imaging of a conventional SSAW configuration, they provide direct information about the in-plane particle distribution patterns and thus enable quantitative comparison of the planar nodal geometries. However, they do not constitute a direct experimental validation of the dual-substrate 3-D configuration considered in the present work. Instead, the out-of-plane nodal structure along the z-direction is examined separately in § 3.3 through comparison between the theoretical prediction and numerical simulation.

3.3. The PN locations along the z-direction

As discussed in the previous section, a planar pressure nodal plane perpendicular to the z-axis is defined by $f_{2}(z)=0$ , written as

(3.6) \begin{align} \dfrac{W_{t}\exp \left(i\varphi _{z}\right)-W_{b}\exp \left(-iHk_{\bot }\right)}{\sin \left(Hk_{\bot }\right)}\cos \left(zk_{\bot }\right)-iW_{b}\exp \left(-izk_{\bot }\right)=0. \end{align}

Let the imaginary and real parts of (3.6) be equal to zero

\begin{align*} & \hspace{8.5pc} \begin{cases} 0=W_{t}\cos (\varphi _{z})\cos (zk_{\bot })-W_{b}\cos ((H-z)k_{\bot }), & \qquad \, (3.7\textit{a})\\[7pt] 0=W_{t}\sin (\varphi _{z}). & \qquad \, (3.7\textit{b}) \end{cases} \end{align*}

Figure 6.

Influence of the fluid cavity height H on the out-of-plane PN position l. The amplitude ratio is set to $\xi =-0.7$ . Red crosses denote the nodal positions extracted from simulations, and the black solid curves show the analytical prediction of l as a function of H. The region with light-blue background indicates the physically meaningful fluid domain, and the theoretical solutions fall outside this region are not considered.

Figure 6. Long description

The horizontal axis is cavity height H and the vertical axis is the out-of-plane pressure-node position l, both in micrometres, for amplitude ratio xi equal to minus 0.7. The light-blue triangular region denotes physically allowable positions between the bottom and top of the cavity. Black curves give the analytical solutions and red crosses give finite-element results; the crosses closely follow the curves. Multiple branches appear as the cavity height increases, indicating additional vertically separated pressure nodes. Insets show pressure fields for H equal to 50, 100, 150, 200, and 250 micrometres. The branches repeat with increasing height, consistent with the predicted vertical periodicity.

The condition for (3.7b ) to hold is $\varphi _{z}=n\pi$ , $n\in \mathbb{Z}$ . Since $\varphi _{z}=(2n+1)\pi$ and $\varphi _{z}=2n\pi$ simply change the sign of $\cos (\varphi _{z})$ in (3.7a ), both sets of values satisfy the equation. Therefore, we take $\varphi _{z}=0$ for simplicity and define a non-dimensional amplitude ratio between SAWs on the top and bottom substrate $\xi =W_{t}/W_{b}$ , which allows us to rewrite (3.7a ) as

(3.8a) \begin{equation}l=z| _{{p_{1}}=0}=\frac{1}{k_{\bot }}\arctan \left(\frac{\xi -\cos \left(Hk_{\bot }\right)}{\sin \left(Hk_{\bot }\right)}\right)+\frac{n\pi }{k_{\bot }},\quad n\in \mathbb{Z} ,\end{equation}

where $l$ is the position where nodal plane intersects z-axis. Inverting (3.8a ) with respect to the SAW amplitude ratio $\xi$ yields

(3.8b) \begin{equation}\xi =\frac{W_{t}}{W_{b}}=\frac{\cos \left(\left(H-l\right)k_{\bot }\right)}{\cos \left(lk_{\bot }\right)}.\end{equation}

Equation (3.8a ) enables direct calculation of the vertical nodal-plane position $l$ for acoustofluidic devices with arbitrary channel height H and SAW amplitude ratio $\xi$ . Conversely, if one aims to manipulate the location l of the vertical nodal plane, (3.8b ) can be applied to determine the corresponding amplitude ratio $\xi$ .

Firstly, the influence of channel height H on the nodal-plane position l is investigated in figure 6. By using a representative amplitude ratio $\xi =-0.7$ , black curves in figure 6 illustrate the theoretical nodal-plane position l as a function of channel height H. The numerical simulation results, denoted by red crosses, show a maximum deviation of less than $0.5 \,{\unicode[Arial]{x03BC}}\mathrm{m}$ against the theoretical calculation, demonstrating good agreement. The term $n\pi /k_{\bot }$ in (3.8a ) indicates a series of equidistant curves, consistent with numerical observations of multiple PNs spaced by $\pi /k_{\bot }=60 \,{\unicode[Arial]{x03BC}}\mathrm{m}$ in channels with $H\gt 100 \,{\unicode[Arial]{x03BC}}\mathrm{m}$ ( $k_{\bot }=52395\,\mathrm{m}^{-1}$ is derived from (2.35)). Furthermore, the dependence on $\cos (Hk_{\bot })$ and $\sin (Hk_{\bot })$ implies a periodicity of $Hk_{\bot }=2\pi$ , and therefore, channels with a height difference of ${\Delta} H=2\pi /k_{\bot }=120 \,{\unicode[Arial]{x03BC}}\mathrm{m}$ will exhibit vertical nodal planes at identical positions. With this theoretical framework, (3.8a ) offers rigorous design guidelines for acoustofluidic devices that require precise out-of-plane acoustic PN manipulation.

In addition to the geometric effects, the amplitude ratio between SAWs on the upper and lower substrates $\xi$ also plays a critical role in controlling the location of the out-of-plane pressure nodal surface. Figure 7 illustrates the qualitative and quantitative results of the calculated nodal-plane positions l for a range of amplitude ratios $\xi$ . The channel height H is fixed at $50 \,{\unicode[Arial]{x03BC}}\mathrm{m}$ for clarity, as only a single nodal plane is presented at this height. The solid lines represent the theoretically calculated nodal-plane positions, while the red crosses indicate the simulation results. Due to the significant variation in nodal position l when $\xi$ spans from −2 to 1, seven, equally spaced amplitude ratios (with an interval of 0.5) were selected for numerical simulation across this range. The maximum deviation between theory and simulation is less than $0.5 \,{\unicode[Arial]{x03BC}}\mathrm{m}$ , demonstrating excellent agreement.

Figure 7.

Influence of the amplitude ratio $\xi$ on the out-of-plane PN position l. Panels $ \unicode{x2460}$ – $ \unicode{x2464}$ in (a) show cross-sectional views of representative simulated acoustic pressure fields under different amplitude ratios $ \xi $ . In (b), the solid black curve gives the analytical prediction for a water-filled cavity with sinusoidally vibrating SAWs on substrates, and the red crosses mark the corresponding finite-element simulation results. The dashed black curve shows the idealised limit in which the substrates vibrate rigidly and $k_{\bot }=k_{\textit{f}}$ , and the blue circles indicate the numerical results obtained under this limit.

Figure 7. Long description

The left column shows five cross-sectional pressure-magnitude fields for representative amplitude ratios xi, with the pressure-node distance l measured from the bottom wall. As xi changes, the nodal plane moves through the 50-micrometre-high cavity; near xi equal to minus one, antiphase actuation produces high pressure at the mid-plane and nodes near the walls. The right graph plots l against the amplitude ratio W sub t divided by W sub b from minus 3 to 3. The solid black curve is the analytical result for sinusoidally vibrating substrates and the red crosses are finite-element results, agreeing to within about 0.5 micrometres. A dashed black curve and blue circles show the rigid-substrate limit k perpendicular equal to k f; this approximation differs most strongly around xi equal to minus one.

A particular interesting case occurs when $\xi =-1$ , corresponding to antiphase vibration of the upper and lower substrates. Under this condition, the central region of the fluid experiences simultaneous compression and extension due to the influences from both the upper and lower substrates, resulting in a higher pressure in the channel mid-plane. This behaviour is consistent with the formation of ANs in the central region and PNs near the bottom and the ceiling of channel, as demonstrated in figure 7(a).

In the previous derivations of acoustic field for inviscid fluids, the parameter $k_{\bot }$ is particularly defined in (2.35) and adopted throughout our work. Its physical significance has been previously discussed in detail by Li et al. (Reference Li, Liang, Kabla, Zhang and Yang2023), and here we briefly clarify its relation to the Rayleigh leakage angle. According to the Huygens–Fresnel principle (Devendran et al. Reference Devendran, Collins, Ai and Neild2017), a SAW propagating on a solid substrate experiences refraction when it enters the fluid. The refraction angle, namely Rayleigh angle, is determined by the wavenumber in water $k_{\textit{f}}=({\omega }/{c})$ and the wavenumber on the substrate $k_{\textit{SAW}}$ . Following the notations in Li et al. (Reference Li, Liang, Kabla, Zhang and Yang2023), the Rayleigh angle is expressed as

(3.9) \begin{equation}\theta _{R}=\arcsin \left(\frac{k_{\textit{SAW}}}{k_{\textit{f}}}\right)\!.\end{equation}

When a pair of counter-propagating travelling surface acoustic waves interferes on the same substrate, both wavefronts refract as they enter the fluid. The wavefronts of the leaked SAWs in the fluid then interfere to form a standing acoustic field, whose effective wavenumber is

(3.10) \begin{equation}k_{\bot }=\frac{k_{\textit{SAW}}}{\tan \left(\theta _{R}\right)}=\frac{k_{\textit{SAW}}}{\tan \left(\arcsin (k_{\textit{SAW}}/k_{\textit{f}})\right)}=\sqrt{{k}_{\textit{f}}^{2}-{k}_{\textit{SAW}}^{2}}.\end{equation}

The result is equivalent to the definition of $k_{\bot }$ in (2.35), thereby confirming the consistency between the present formulation and the previous result reported in Li et al. (Reference Li, Liang, Kabla, Zhang and Yang2023).

Equation (3.9) demonstrates that a smaller $k_{\textit{SAW}}$ results in a smaller Rayleigh angle. Therefore, in the limiting case when $k_{\textit{SAW}}\rightarrow 0$ , the SAW propagation along the substrate undergoes negligible refraction when it couples into the fluid. Physically, this corresponds to $\lambda _{\textit{SAW}}\rightarrow \infty$ , such that the substrates can be interpreted as rigid ‘pistons’ oscillating with high frequency along the z-axis (normal to the surface), and the leaked waves behave as a BAW. Figure 7 shows the calculated nodal positions under this limit, demonstrating a high degree of agreement between the theory and simulation.

3.4. Viscous effects

While (2.42) and (2.43) can accurately describe acoustofluidic phenomena for most fluids, they may be insufficient for scenarios involving highly viscous fluid or demanding high computational accuracy. For instance, in petroleum processing during oil fractionation (Luo et al. Reference Luo, Gong, He, Zhang and He2021; Fasano et al. Reference Fasano, Li, Diez, D’Addesa, Manor, Cummings and Kondic2025) and biomedical engineering for cell sorting in viscous fluids such as blood plasma (Lim et al. Reference Lim, Hong, Cho, Han, Kim, Kim and Nam2025), the incorporation of viscosity via (2.31) and (2.32) provides a more precise description of the underlying physics. Specifically, to examine the boundary-layer effect in water, figure 8 presents a direct comparison of the first-order velocity component $u_{1}$ under viscous and inviscid conditions. As shown in figure 8(a), the overall standing-wave structure and the order of magnitude of $u_{1}$ remain similar in the two cases, indicating that viscosity does not fundamentally alter the overall acoustic field. However, clear differences appear in the immediate vicinity of the vibrating substrates, as highlighted in figures 8(b) and 8(c), where the viscous solution satisfies the no-slip condition and develops a rapid near-wall transition that is absent in the inviscid solution. The penetration depth (Marques & Tomaz Da Silva Reference Marques and Tomaz Da Silva2025) $\delta$ can be estimated to be

(3.11) \begin{equation}\delta =\sqrt{2\mu /\rho _{0}\omega }.\end{equation}

For liquid water, $\delta =0.16 \,{\unicode[Arial]{x03BC}}\mathrm{m}$ , which indicates that the boundary-layer effects are confined within a very thin region adjacent to the substrate. As detailed by Doinikov et al. (Reference Doinikov, Thibault and Marmottant2017), the pronounced velocity overshoot near the edge of the boundary layer arises from the transition between the thin viscous layer and the bulk region, which is evident in figures 8(b) and 8(c). Therefore, the viscosity of water does not change the global standing-wave pattern, but it modifies the near-wall velocity distribution and thereby influences the nodal positions in the $z$ -direction.

Figure 8.

Boundary-layer velocity profiles for water under viscous and inviscid conditions. Solid lines denote the viscous solution, and dotted lines denote the inviscid solution. (a) First-order velocity component $u_{1}$ along the z-axis for different amplitude ratios ξ. (b) Magnified view near the bottom substrate. (c) Magnified view near the top substrate.

Figure 8. Long description

Panel (a) plots the first-order velocity component u1 against z from 0 to 50 micrometres for amplitude ratios xi ranging from minus 2 to 1. Colours identify xi; solid curves are viscous solutions and dotted curves are inviscid solutions. Across most of the cavity the paired curves overlap, while visible deviations occur only close to z equal to 0 and 50 micrometres. Panels (b) and (c) magnify the bottom and top one-micrometre regions. Shaded strips and vertical markers indicate the viscous penetration depth delta. The viscous solution satisfies the wall velocity and rapidly approaches the inviscid profile outside each thin boundary layer, demonstrating that viscosity is localised near the fluid-substrate interfaces in water.

This behaviour can also be interpreted from the eigenvalue analysis in (2.23). As discussed in § 2.5, the eigenvalue $\lambda _{4}=(6428+6428i)\,\mathrm{mm}^{-1}$ is associated with the viscous effect and contains a large real component, indicating rapid spatial attenuation. In contrast, $\lambda _{1}=(0.004+52.4i)\,\mathrm{mm}^{-1}$ , which describes the inviscid effect, has a negligible real component. Owing to the large real component of $\lambda _{4}$ , viscous effects decay exponentially to almost zero in water, whereas the waves in the inviscid bulk region represented by $\lambda _{1}$ can propagate over considerable distances without significant attenuation. This explains why the viscous and inviscid solutions remain close in the bulk region while exhibiting substantial differences near the fluid–substrate interfaces.

For most liquids commonly encountered in acoustofluidic applications, the influence of viscosity can be quantified by the viscous penetration depth $\delta$ and the associated acoustic attenuation. Under typical MHz-frequency conditions, $\delta$ remains much smaller than the channel height and the acoustic wavelength for low-to-moderately viscous fluids. In this situation, viscous effects are mainly confined to the thin boundary layers adjacent to the vibrating substrates, while the acoustic field in the central part of the cavity remains virtually inviscid. As viscosity increases further, however, the penetration depth grows and viscous attenuation becomes more pronounced, so that the effect of viscosity will no longer be restricted to the immediate vicinity of the boundaries but modifies the flow over the whole domain. To illustrate this transition from the boundary-layer-confined viscous effects to more substantial cavity-scale influence, we next consider engine oil as a representative highly viscous fluid.

To further investigate how strong viscosity influences the nodal-plane positions in the acoustic field, we considered a highly viscous material, namely engine oil. The material properties of the engine oil employed in this study were $\rho _{0}=887.59\,\mathrm{kg}\,\mathrm{m}^{-3}$ , $c=1000$ m s⁻¹ (Rooke et al. Reference Rooke, Brunskill, Li, Taghizadeh, Hunter, He, Lu and Dwyer-Joyce2023), $\mu =0.79017\,\mathrm{Pa}\,\mathrm{s}$ (Alhouli et al. Reference Alhouli, AlShammri, Abed and Alali2015; Tripathi & Vinu Reference Tripathi and Vinu2015; Barcenas & Miller Reference Barcenas and Miller2025) and $\mu _{B}=0.79\,\mathrm{Pa}\boldsymbol{\cdot }\mathrm{s}$ . Since precise measurements of $\mu _{B}$ are unavailable, it was approximated by referring to $\mu$ . To assess the sensitivity of the results to this approximation, additional calculations were performed with $\mu _{B}=0.4$ and $1.2\,\mathrm{ Pa}\boldsymbol{\cdot }\mathrm{s}$ , i.e. $\pm 50\mathrm{\,\% }$ around the adopted value. The resulting changes in $\lambda _{1}$ and $\lambda _{2}$ were less than $1\mathrm{\,\% }$ in magnitude and less than $1^{\circ }$ in inclination angle, while the changes in $\lambda _{3}$ and $\lambda _{4}$ were below $10^{-4}$ in magnitude, with negligible angular variation. These results indicate that the conclusions of the present analysis are insensitive to the assumed bulk viscosity.

The decomposition in (3.1) can also be applied to (2.32), yielding

(3.12) \begin{align} \begin{cases} 0=f_{1,\mathrm{vis}}\left(x,y\right)\boldsymbol{\cdot }f_{2,\mathrm{vis}}(z),\\[11pt] f_{1,\mathrm{vis}}\left(x,y\right)=\cos \left(k_{\textit{SAW}}x-\varphi _{x}\right)+\exp \left(-i\varphi _{\varDelta}\right)\cos \left(k_{\textit{SAW}}y-\varphi _{y}\right)\!,\\[11pt] f_{2,\mathrm{vis}}(z)=\sum _{j=1,2,3,4}\left(\dfrac{i\omega \rho _{0}+\mu {k}_{\textit{SAW}}^{2}-\mu \lambda _{j}^{2}}{\lambda _{j}(R-\mu )}\right)\\[11pt]\qquad \left(C_{b,j}W_{b}+C_{t,j}W_{t}\exp \left(i\varphi _{z}\right)\right)\exp \left(\lambda _{j}z\right)\!. \end{cases} \end{align}

Equation (3.12) shows that, for highly viscous fluids, the full 3-D pressure nodal surface can be decomposed into two independent nodal surface systems. From (3.1) and (3.12), $f_{1,\mathrm{vis}}(x,y)=f_{1}(x,y)$ , indicating that the shape of nodal surface parallel to the z-axis is identical for inviscid and viscous conditions. However, the presence of real components in the eigenvalues $\lambda _{1}$ and $\lambda _{2} $ under highly viscous cases (as shown in the table in figure 3) directly cause the nodal system $f_{2,\mathrm{vis}}(z)$ values to deviate from those in inviscid fluid (described by $f_{2}(z)$ in (3.1)). Therefore, the introduction of viscosity mainly affects the PN locations along the z-direction. In figure 9, the acoustic pressure amplitude along the z-direction is captured under different cases. The black solid lines indicate the inviscid condition obtained by (2.43), while the blue and red dashed lines denote the viscous scenarios calculated from (2.32) with and without a phase correction, respectively. From figure 9, substantial discrepancies near the nodal positions can be observed between the inviscid and viscous curves. This indicates that the inviscid formulation in (2.43) is no longer sufficient to accurately capture the acoustic field under highly viscous conditions.

Figure 9.

Acoustic pressure profiles along the z-direction for a high-viscosity fluid under different amplitude ratios $\xi$ . The black solid curves are the inviscid prediction from (2.43), the blue dashed curves show the viscous results from (2.32) and the red dotted curves are viscous prediction obtained from (2.32) with additional phase compensation $\varphi _{z}$ to reduce the viscosity-induced residual pressure at the nodal positions.

Figure 9. Long description

Six stacked panels plot the acoustic-pressure magnitude along z from 0 to 50 micrometres for xi equal to minus 2, minus 1.5, minus 1, minus 0.5, 0.5, and 1. Black solid curves are inviscid predictions, blue dashed curves are viscous predictions without phase correction, and red dotted curves are viscous predictions with phase compensation phi z. Viscous attenuation shifts the minima and leaves non-zero residual pressure at the nominal nodal positions, especially for negative amplitude ratios. The phase-compensated red curves restore minima close to zero and generally follow the inviscid nodal locations, while retaining viscosity-related changes in amplitude.

As depicted in figure 9, the blue dashed lines highlight another challenge posed by high-viscosity fluids: the absence of locations where the acoustic pressure reaches an exact zero. In an inviscid fluid, the two counter-propagating waves along the out-of-plane direction do not undergo amplitude attenuation, and their superposition always produces exact zero-pressure locations regardless of the amplitude ratio $\xi$ , which can also be concluded from (3.8). For viscous fluid, however, solely tunning amplitude ratio $\xi$ cannot guarantee a zero pressure at the nodes, thereby impacting the efficiency and accuracy of acoustofluidic manipulation. To restore a zero pressure on the nodal plane in viscous environment, one can adjust $\varphi _{z}$ corresponding to the amplitude ratio $\xi$ . Specifically, setting $f_{2,\mathrm{vis}}(z)$ in (3.12) to zero yields

\begin{align*} & \begin{cases} \exp \left(i\varphi _{z}\right)=-\dfrac{1}{ \xi } \dfrac{\sum _{j=1,2,3,4}\left(\dfrac{i\omega \rho _{0}+\mu {k}_{\textit{SAW}}^{2}-\mu \lambda _{j}^{2}}{\lambda _{j}(R-\mu )}\right)C_{b,j}\exp \left(\lambda _{j}l\right)}{\sum _{j=1,2,3,4}\left(\dfrac{i\omega \rho _{0}+\mu {k}_{\textit{SAW}}^{2}-\mu \lambda _{j}^{2}}{\lambda _{j}(R-\mu )}\right)C_{t,j}\exp \left(\lambda _{j}l\right)}, &\quad (3.13\textit{a})\\[17pt] 1=\dfrac{1}{ |\xi| } \left|\left|\dfrac{\sum _{j=1,2,3,4}\left(\dfrac{i\omega \rho _{0}+\mu {k}_{\textit{SAW}}^{2}-\mu \lambda _{j}^{2}}{\lambda _{j}(R-\mu )}\right)C_{b,j}\exp \left(\lambda _{j}l\right)}{\sum _{j=1,2,3,4}\left(\dfrac{i\omega \rho _{0}+\mu {k}_{\textit{SAW}}^{2}-\mu \lambda _{j}^{2}}{\lambda _{j}(R-\mu )}\right)C_{t,j}\exp \left(\lambda _{j}l\right)}\right|\right|\!. &\quad (3.13\textit{b}) \end{cases} \end{align*}

Equation (3.13b ) determines whether a zero-pressure nodal plane can exist and where plane is located, while (3.13a ) gives the required phase compensation $\varphi _{z}$ to achieve a zero pressure at the nodal plane. Since (3.13b ) is a function of only $l$ , the Newton–Raphson method is implemented to determine the nodal-plane position $l$ . Then, the nodal-plane position $l$ acquired from (3.13b ) is substituted into (3.13a ) to obtain the corresponding compensation phase, $\varphi _{z}$ .

By adjusting the compensation phase $\varphi _{z}$ , as illustrated by the red dashed line in figure 9, the sharp zero-pressure wells can be restored. This confirms that amplitude–phase modulation is essential for achieving zero acoustic pressure nodal planes in highly viscous fluids. Furthermore, it is worth noting that adjusting $\varphi _{z}$ does not provides a global correction, instead it only minimises acoustic pressure at a specific targeted PN. For example, when ${\unicode[Arial]{x03BE}} =-1$ , two PNs appear at $z=6 $ and $z=44 \,{\unicode[Arial]{x03BC}}\mathrm{m}$ . Appling compensation to the PN at $z=6 \,{\unicode[Arial]{x03BC}}\mathrm{m}$ inevitably increases the pressure at $z=44 \,{\unicode[Arial]{x03BC}}\mathrm{m}$ . This selective compensation strategy for specific nodal positions suggests a potential approach for targeted acoustofluidic particle manipulation in viscous environments.

3.5. Second-order velocity field

The analytical formulation developed in § 2.7 is further validated by comparing the predicted second-order acoustic streaming field with numerical simulations. The same fluid domain and material parameters as those used in the first-order validation are adopted. The first-order acoustic fields for the viscous fluid obtained from (2.31) and (2.32) are used to construct the nonlinear forcing terms in (2.44), and the resulting acoustic streaming field is calculated analytically from the solution given in (2.65). A representative two-dimensional $(x,z)$ subproblem is considered to validate the second-order streaming formulation by comparing the analytical results with the numerical simulations. Figure 10 shows the second-order velocity component $w_{2}$ together with the corresponding streamlines for two representative amplitude ratios, $\xi =1$ and $\xi =-0.7$ .

Figure 10.

Validation of the analytical solution for the second-order acoustic streaming field. Contours show the second-order velocity component ( $w_{2}$ ), with streamlines indicating the corresponding quasi-steady streaming patterns. Red crosses denote the numerical results, and black solid lines denote the analytical solutions along the selected cutlines. (a) Results for $\xi =1$ . (b) Results for $\xi =-0.7$ .

Figure 10. Long description

Panel (a), for xi equal to 1, shows a symmetric second-order vertical-velocity field w2 with four counter-rotating streaming cells and a colour range of approximately minus 42 to 42 micrometres per second. Panel (b), for xi equal to minus 0.7, shows a different, lower-magnitude two-cell pattern with a range of approximately minus 12 to 12 micrometres per second. For each case, two graphs to the right plot w2 along selected vertical cutlines. Red crosses are numerical results and black solid curves are analytical solutions. The profiles nearly overlap; the reported coefficients of determination are about 0.98 for both cutlines in panel (a), and 0.997 and 0.98 in panel (b).

As shown in figure 10(a), when $\xi =1$ , the streaming field forms a periodic array of counter-rotating vortices. The magnitude of $w_{2}$ reaches approximately $42 \,{\unicode[Arial]{x03BC}}\mathrm{m}$ s⁻¹, which is several orders of magnitude smaller than the first-order acoustic velocity, consistent with its second-order origin. The analytical profiles along the selected cutlines agree well with the numerical results, with ( $R^{2}=0.98$ ) for both comparisons. For $\xi =-0.7$ , as shown in figure 10(b), the streaming pattern is modified by the change in the relative vibration amplitude of the two substrates. The maximum magnitude of $w_{2}$ decreases to approximately $12 \,{\unicode[Arial]{x03BC}}\mathrm{m}$ s⁻¹, and the vortex structure becomes less symmetric compared with the case of $\xi =1$ . Nevertheless, the analytical solution still captures the numerical streaming profiles accurately, with $R^{2}=0.997$ and $R^{2}=0.98$ along the two selected cutlines.

These comparisons demonstrate that the frequency-domain perturbation method can accurately predict not only the first-order acoustic field but also the induced time-averaged second-order streaming flow. It should be noted that acoustic streaming is highly sensitive to the near-wall boundary-layer flow. Therefore, in the numerical simulations, the mesh was refined near the solid boundaries to better resolve the boundary-layer contribution to the streaming field, while a relatively coarser mesh was used in the central region of the channel to reduce the computational cost. This mesh distribution is also reflected in figure 10. Near the boundaries, the numerical results vary smoothly and agree closely with the analytical predictions. In contrast, in the central region, where the mesh is relatively sparse, the numerical profiles exhibit larger local fluctuations and deviations from the analytical curves. These discrepancies are therefore likely associated with the finite numerical resolution in the central region, while the overall agreement shown in figure 10 confirms the reliability of the proposed analytical framework for describing the representative 2-D acoustic streaming subproblem induced by SAW excitation.