2.1 Basic notation and kinematics

We begin by distinguishing two distinct time scales with corresponding time variables $t$ and $T$ . We will refer to $t$ as the fast time characterizing the acoustic wave time scale. We will refer to $T$ as the slow time characterizing the mean motion of the fluid. A distinction between these two times can be done by adopting the following relation (Xie & Vanneste Reference Xie and Vanneste2014b ):

where $\unicode[STIX]{x1D716}>0$ is a non-dimensional smallness parameter typically defined as a ratio of length scales (cf. Köster Reference Köster2006; Ding et al. Reference Ding, Li, Lin, Stratton, Nama, Guo, Slotcavage, Mao, Shi and Costanzo2013), and $\unicode[STIX]{x1D6FD}>1$ is a non-dimensional exponent, yet to be determined, whose value sets the meaning of time-scale separation in a rigorous way. We note that, in the approach followed by Xie & Vanneste (Reference Xie and Vanneste2014b ), the description of the physical model proceeds by introducing functions of $t$ only and letting functions of $T$ emerge through proper time averages. This is not the only way to present the effects of the time-scale separation. Another way to proceed is illustrated, for example, by Chini, Malecha & Dreeben (Reference Chini, Malecha and Dreeben2014). Here functions of the type $f(x,t)$ ( $x$ representing space) are re-defined as $f(x,t,T)$ where the relation between $t$ and $T$ is exactly that indicated above. If one adopts this representation, then the time-scale separation is also manifested through an additive representation of time derivatives. Specifically, we have the remapping $\unicode[STIX]{x2202}_{t}\rightarrow \unicode[STIX]{x2202}_{t}+\unicode[STIX]{x1D716}^{\unicode[STIX]{x1D6FD}}\unicode[STIX]{x2202}_{T}$ .

Next, referring to figure 1, we denote by $B_{0}\subset \mathscr{E}^{d}$ , $d=2,3$ , the reference configuration of the fluid, where $\mathscr{E}^{d}$ is the underlying $d$ -dimensional Euclidean point space in which the body moves. For simplicity $B_{0}$ will be taken as the initial configuration of the fluid (before acoustic excitation), assumed to be a quiescent state. We denote by $\mathscr{V}$ the translation (vector) space of $\mathscr{E}^{d}$ . As is typical in an Euclidean context, we define a second-order tensor to be an element of $\mathscr{L}(\mathscr{V},\mathscr{V})$ (Gurtin, Fried & Anand Reference Gurtin, Fried and Anand2010), i.e. a linear map from $\mathscr{V}$ into $\mathscr{V}$ .

Figure 1.

Schematic of the domains we consider. Here $B_{0}$ denotes the reference configuration of our system and $B(t)$ denotes the current configuration of the system. We view $B(t)$ as the outcome of a deformation $\unicode[STIX]{x1D74C}_{\unicode[STIX]{x1D709}}$ from an intermediate configuration $\mathscr{B}(T)$ , which we view as the configuration of the body whose motion is observed at the slow time scale.

As explained later, the system’s motion will be viewed as a sum of contributions of different orders. Given a quantity, say $\unicode[STIX]{x1D719}$ , the contribution of $n$ th order to this quantity will be denoted as $\unicode[STIX]{x1D719}^{(n)}$ . Finally, given a tensor $\unicode[STIX]{x1D63C}$ , its transpose and its inverse transpose will be denoted by $\unicode[STIX]{x1D63C}^{\text{T}}$ and $\unicode[STIX]{x1D63C}^{-\text{T}}$ , respectively.

We denote material fluid particles in $B_{0}$ by $\boldsymbol{X}$ . The current configuration of the fluid resulting from the acoustic excitation will be described as follows:

where, as is standard in continuum mechanics, for all $t>0$ , $\unicode[STIX]{x1D73B}(\boldsymbol{X},\bullet )$ is assumed to be a diffeomorphism between $B_{0}$ and $B(t)$ , and where, for all $\boldsymbol{X}\in B_{0}$ , $\unicode[STIX]{x1D73B}(\bullet ,t)\in C^{2}(\mathbb{R}^{+})$ . We assume that the motion $\unicode[STIX]{x1D73B}(\boldsymbol{X},t)$ is the composition of two motions, with smoothness properties in space and time so as to guarantee the assumed properties of $\unicode[STIX]{x1D73B}(\boldsymbol{X},t)$ . These motions are as follows:

and

such that, under the constraint in (2.1),

Here $\unicode[STIX]{x1D74C}(\boldsymbol{X},T)$ and $\unicode[STIX]{x1D74C}_{\unicode[STIX]{x1D709}}(\boldsymbol{x},t)$ will be referred to as the slow and the fast motion, respectively. Since we have defined $B_{0}$ , $\mathscr{B}(T)$ , and $B(t)$ as subsets of the same underlying Euclidean point space, we introduce corresponding displacement fields as follows:

where, as discussed later, $\boldsymbol{u}^{(0)}$ is the zeroth-order contribution to $\boldsymbol{u}$ . In view of (2.1) and (2.5), we have that, for $t>0$ and $\boldsymbol{X}\in B_{0}$ ,

The displacement field $\unicode[STIX]{x1D743}(\boldsymbol{x},t)$ is assumed to be a harmonic function of $t$ , with period equal to that of the harmonic excitation of the system.

To make the notation less cumbersome, in the remainder of the paper we will work under the conditions underlying the composition in (2.5). Furthermore, by the notation $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D719}}(\bullet )$ we denote the partial derivative of the quantity $(\bullet )$ with respect to the variable $\unicode[STIX]{x1D719}$ , where the latter can be scalar-, vector-, or tensor-valued. With the notation $\unicode[STIX]{x1D735}_{z}(\bullet )$ we denote the spatial gradient of $(\bullet )$ with respect to the position variable $\boldsymbol{z}$ , i.e.  $\unicode[STIX]{x1D735}_{z}(\bullet )=\unicode[STIX]{x2202}_{z}(\bullet )$ ; and by $\unicode[STIX]{x1D735}_{z}\boldsymbol{\cdot }(\bullet )$ we denote the corresponding divergence operator. Finally, the notation $\dot{\unicode[STIX]{x1D719}}$ will denote the material time derivative of the quantity $\unicode[STIX]{x1D719}$ , where the latter can be scalar-, vector-, or tensor-valued, and where by material time derivative we mean (cf. Gurtin et al. Reference Gurtin, Fried and Anand2010),

The deformation gradient of the slow motion is

Similarly, the deformation gradient of the fast (harmonic) motion is

where $\unicode[STIX]{x1D644}$ is the identity tensor and $\mathscr{L}^{+}(\mathscr{V},\mathscr{V})$ is the set of tensors with positive determinant. The positivity of $\det (\unicode[STIX]{x1D641}^{(0)})$ and $\det (\unicode[STIX]{x1D641}_{\unicode[STIX]{x1D709}})$ is a consequence of having assumed $\unicode[STIX]{x1D74C}_{\unicode[STIX]{x1D709}}$ to be a diffeomorphism for all $T>0$ and that $B_{0}\equiv B(0)$ . The deformation gradient of the overall motion $\unicode[STIX]{x1D73B}$ is then given as (Gurtin et al. Reference Gurtin, Fried and Anand2010)

For future use, we will denote $\det (\unicode[STIX]{x1D641}_{\unicode[STIX]{x1D709}})$ by $J_{\unicode[STIX]{x1D709}}$ , and by recalling the formula for the characteristic polynomial of a second-order tensor, we have the following (exact) representation formula for $J_{\unicode[STIX]{x1D709}}$ :

where $\mathscr{I}_{k}$ , $k=1,\ldots ,3$ , are the principal invariants of $\unicode[STIX]{x1D735}_{x}\unicode[STIX]{x1D743}$ , namely

Using (2.9), the material velocity field is

where, as implied by the first of (2.16), $\boldsymbol{v}$ is a field over $B(t)$ . For later use we define the following velocity field:

2.2 Governing equations

The motion of the fluid is governed by the balance laws for mass and momentum (Gurtin et al. Reference Gurtin, Fried and Anand2010) given by, respectively,

where $\unicode[STIX]{x1D70C}(\boldsymbol{y},t):B(t)\rightarrow \mathbb{R}$ is the mass density distribution, $\boldsymbol{b}(\boldsymbol{y},t):B(t)\rightarrow \mathscr{V}$ is the external body force density per unit mass, and $\unicode[STIX]{x1D64F}(\boldsymbol{y},t):B(t)\rightarrow \mathscr{L}(\mathscr{V},\mathscr{V})$ is the Cauchy stress.

We assume the fluid to be linear, viscous, and compressible with constitutive response function for the Cauchy stress given by

where $p$ is the fluid (static) pressure, $\unicode[STIX]{x1D707}$ and $\unicode[STIX]{x1D707}_{b}$ are constants representing the shear and bulk viscosities, respectively, and $p(\unicode[STIX]{x1D70C})$ is assumed to be the following linear relation

where $c_{0}$ and $\unicode[STIX]{x1D70C}_{0}$ are constants denoting the fluid’s speed of sound and mass density at rest, respectively.

2.3 Reformulation of the governing equations

We now reformulate the governing equations by mapping them onto body’s mean deformed configuration $\mathscr{B}(T)$ . A field defined over $B(t)$ mapped onto $\mathscr{B}(T)$ will be identified by an asterisk, i.e. given $\unicode[STIX]{x1D719}(\boldsymbol{y},t):B(t)\rightarrow \mathscr{W}$ , with $\unicode[STIX]{x1D719}$ a scalar-, vector-, or tensor-valued quantity, and $\mathscr{W}$ some appropriate vector space, $\unicode[STIX]{x1D719}^{\ast }(\boldsymbol{x},t):\mathscr{B}(T)\rightarrow \mathscr{W}$ will denote the corresponding mapped quantity according to the following definition:

A field defined over $B_{0}$ mapped onto $\mathscr{B}(T)$ will be denoted with the superscript symbol $\#$ , i.e. given $\unicode[STIX]{x1D719}(\boldsymbol{X},t):B_{0}\rightarrow \mathscr{W}$ , with $\unicode[STIX]{x1D719}$ a scalar-, vector-, or tensor-valued quantity, and $\mathscr{W}$ some appropriate vector space, $\unicode[STIX]{x1D719}^{\#}(\boldsymbol{x},t):\mathscr{B}(T)\rightarrow \mathscr{W}$ will denote the corresponding mapped quantity according to the following definition:

Using the above definition, equations (2.18) and (2.19) can be expressed on $\mathscr{B}(T)$ as (see appendix A for details)

where, with a slight abuse of notation, $\unicode[STIX]{x1D64B}^{\ast }$ is the Piola–Kirchhoff stress tensor on $\mathscr{B}(T)$ (cf. Gurtin et al. Reference Gurtin, Fried and Anand2010) defined as follows:

For later use, it is also convenient to consider the Lagrangian expression of the balance of mass. This reads $\unicode[STIX]{x1D70C}(\boldsymbol{y},t)\det \unicode[STIX]{x1D641}(\boldsymbol{X},t)=\unicode[STIX]{x1D70C}_{0}(\boldsymbol{X})$ with $\boldsymbol{y}=\unicode[STIX]{x1D73B}(\boldsymbol{X},t)$ (Gurtin et al. Reference Gurtin, Fried and Anand2010). Mapping this expression onto $\mathscr{B}(T)$ , and in view of (2.13) and (2.23), we have

2.4 Time average, time-scale separation, and orders of magnitude

As anticipated earlier, to model the behaviour of acoustofluidic devices, we assume that the displacement field $\unicode[STIX]{x1D743}(\boldsymbol{x},t)$ is the time-harmonic response to the time-harmonic excitation of the system. We will denote the period of this excitation by $\unicode[STIX]{x1D6F1}$ . As discussed later, the governing equations will be subjected to an asymptotic expansion using the smallness parameter $\unicode[STIX]{x1D716}$ . By grouping terms of equal order in $\unicode[STIX]{x1D716}$ one obtains a hierarchy of problems that will be shown to be linear in the underlying unknowns. The harmonic response we have just referred to is that applied through the boundary conditions of the first-order problem. The second-order problem is then derived through a time-averaging operation that necessitates the choice of a physically meaningful time period. For simplicity, here we present a treatment in which there is only one excitation frequency. Thus, the smallest time period over which time averaging makes sense coincides with the period of the harmonic forcing. If the harmonic forcing consists of the sum of contributions with different frequencies, due to the linearity of the first-order problem, one can obtain a first-order solution through superposition. If one further assumes that the frequencies in question are multiples of the smallest of such frequencies, then the averaging underlying the second-order problem would have to be carried out over the longest period, namely that corresponding to the smallest frequency. More complex scenarios in which the forcing is, say, quasi-periodic, require additional analysis that might be explored in future work.

We now introduce the notation to denote the time-average operation over the excitation period $\unicode[STIX]{x1D6F1}$ . For any generic quantity $\unicode[STIX]{x1D719}(\bullet ,t)$ (scalar-, vector-, or tensor-valued), the notation $\langle \unicode[STIX]{x1D719}\rangle$ will identify its time average over the period of harmonic excitation, i.e.

We note that, having assumed that $\unicode[STIX]{x1D743}(\boldsymbol{x},t)$ is a harmonic function of time with period $\unicode[STIX]{x1D6F1}$ , we have

A similar conclusion applies to any $\unicode[STIX]{x1D6F1}$ -periodic harmonic function. On the other hand, we observe that $\unicode[STIX]{x1D641}_{\unicode[STIX]{x1D709}}$ is merely $\unicode[STIX]{x1D6F1}$ -periodic and therefore $\langle \unicode[STIX]{x1D641}_{\unicode[STIX]{x1D709}}\rangle$ does not vanish. In fact, from the last of (2.12)

For future discussion, it is also important to be able to relate $\langle \unicode[STIX]{x1D719}\rangle$ for a generic function $\unicode[STIX]{x1D719}(\bullet ,t)$ to the smallness parameter $\unicode[STIX]{x1D716}$ . To this end, referring to (2.1), we set

The above expression is a simple re-expression of a quantity from the fast to the slow time scale. Next, assuming that $\unicode[STIX]{x1D719}(\bullet ,t)$ is differentiable in time, in the neighbourhood of a time instant $t=t_{0}$ , Taylor’s theorem allows us to write

Therefore, letting $T_{0}=\unicode[STIX]{x1D716}^{\unicode[STIX]{x1D6FD}}t_{0}$ , for such a $\unicode[STIX]{x1D719}(\bullet ,t)$ , we can express $\langle \unicode[STIX]{x1D719}\rangle (\bullet ,t)$ in a neighbourhood of $t_{0}$ as

Next, we observe that

Now, letting $\unicode[STIX]{x0394}T=O(1)$ be some meaningful time duration at the slow scale, in the study of microacoustofluidics it is reasonable to posit that the excitation period is a reference time interval for distinguishing fast motions from slow motions. Furthermore, in view of the definition in (2.28), it is reasonable to limit the context of the expansion in (2.33) to sub-intervals of the period $\unicode[STIX]{x1D6F1}$ . Therefore we posit that

With this stipulation, equations (2.33) and (2.34) imply that, at worst,

that is, the time-average operation is, at worst, equivalent to the remapping at the slow scale up to a remainder on the order of $\unicode[STIX]{x1D716}^{2\unicode[STIX]{x1D6FD}}$ .

As indicated in the Introduction, we will follow the approach in Nyborg (Reference Nyborg, Hamilton and Blackstock1998) (see also Köster Reference Köster2007; Vanneste & Bühler Reference Vanneste and Bühler2011; Muller et al. Reference Muller, Barnkob, Jensen and Bruus2012; Nama et al. Reference Nama, Barnkob, Mao, Kähler, Costanzo and Huang2015) where the response to a harmonic excitation is analysed via an asymptotic expansion of the velocity in contributions of different integer orders of the smallness parameter $\unicode[STIX]{x1D716}$ . Our approach differs from previous ones in that this expansion is done on the equations written over $\mathscr{B}(T)$ instead of $B(t)$ . Furthermore, as done in Xie & Vanneste (Reference Xie and Vanneste2014b ), we add the displacement to the fields to be expanded, with the understanding that the displacement field must be allowed to be on the order of the device’s largest linear dimension. Hence, we set the displacement field describing the mean motion of the system to be of order $O(\unicode[STIX]{x1D716}^{0})=O(1)$ , and we have already denoted this field by $\boldsymbol{u}^{(0)}$ . The remaining contribution to the displacement is then the harmonic contribution $\unicode[STIX]{x1D743}$ , which we assume to be $O(\unicode[STIX]{x1D716})$ . Also, as it is customary, we will assume that differentiation with respect to any of the position variables and with respect to the fast time $t$ does not alter the order of a quantity. Here, we remark that we are following in the footsteps of analyses like that in Köster (Reference Köster2007). In reality, this is not the only way to carry out asymptotic expansions, as shown in Xie & Vanneste (Reference Xie and Vanneste2014b ) in which very interesting alternative scenarios have been identified. The current work will explore such scenarios in future publications.

Summarizing, we have

With reference to (2.14), we also have

We will then follow previous approaches (Köster Reference Köster2007; Vanneste & Bühler Reference Vanneste and Bühler2011; Muller et al. Reference Muller, Barnkob, Jensen and Bruus2012; Nama et al. Reference Nama, Barnkob, Mao, Kähler, Costanzo and Huang2015) in assuming that:

1. (i) there is no order zero (background) velocity field;

2. (ii) the harmonic response in terms of velocity is of order $\unicode[STIX]{x1D716}$ ;

3. (iii) and the streaming flow is at most of order $\unicode[STIX]{x1D716}^{2}$ .

Using these assumptions, and enforcing consistency between the material time derivative of the displacement expansion and that of the velocity, yields a constraint equation for $\unicode[STIX]{x1D6FD}$ that determines the separation of time scales. Specifically, referring to (2.16), the mapped material velocity is given by

In view of (2.37), for the first-order contribution to be only $\boldsymbol{v}_{\unicode[STIX]{x1D709}}$ , the last term must necessarily be $O(\unicode[STIX]{x1D716}^{2})$ . This result matches the analysis of the ‘non-inertial regime’ described by Xie & Vanneste (Reference Xie and Vanneste2014b ). We then formally set

for which the first-order contribution to the material velocity is

whereas the second-order contribution is defined as follows:

As far as the mass density distribution is concerned, we simply postulate the existence of the zeroth-, first-, and the second-order contributions. With this said, we make the additional assumption that the time scale at which the second-order contribution is meaningful is the slow time scale. In other words, we posit that

Finally, we assume that the mass density distributions and external body force field are as follows:

2.5 Mean fields and Lagrangian-mean fields

Before presenting the hierarchy of boundary value problems (BVPs) generated by the assumptions stated in the previous section, we now revisit our kinematics and discuss the rationale for our assumptions.

In the work of Bradley (Reference Bradley1996, Reference Bradley2012), Nyborg (Reference Nyborg, Hamilton and Blackstock1998), Köster (Reference Köster2007), Vanneste & Bühler (Reference Vanneste and Bühler2011), Muller et al. (Reference Muller, Barnkob, Jensen and Bruus2012), Nama et al. (Reference Nama, Barnkob, Mao, Kähler, Costanzo and Huang2015) (and the many works that adopt the same kinematic framework) the domain of the equations governing the response of the fluid under a harmonic excitation is a time-independent domain (or control volume), which we will call $\unicode[STIX]{x1D6FA}$ . The fields defined over $\unicode[STIX]{x1D6FA}$ are expressed as functions of time $t$ and a position variable, which we will denote by $\boldsymbol{z}$ . The purpose of the present discussion is to compare and contrast the kinematics adopted in the paper with that of previous works.

We begin by comparing the meaning of the variable $\boldsymbol{z}$ with the three descriptors used here, namely, $\boldsymbol{X}$ , $\boldsymbol{x}$ , and $\boldsymbol{y}$ . Similar to what we have done, in the aforementioned references there is a slow and a fast time scale, and quantities at the slow time scale are averages of analogous quantities at the fast time scale. In the generalized Lagrangian mean (GLM) theory (Andrews & McIntyre Reference Andrews and Mcintyre1978; Bühler Reference Bühler2009) this averaging concept is also applied to the position $\boldsymbol{z}$ itself. Let us recall that $\boldsymbol{z}$ is viewed as the position at time $t$ of some material particle $P$ , with the stipulation that at time $t+\unicode[STIX]{x0394}t>t$ the position $\boldsymbol{z}$ will be occupied by a particle $Q$ such that, in general, $Q\neq P$ . If this understanding of $\boldsymbol{z}$ is the only one used, then $\boldsymbol{z}$ does not have a clear correspondence with any of the three descriptors introduced here. The situation is different if we interpret $\boldsymbol{z}$ as is done in the GLM theory, which posits that $\boldsymbol{z}$ is the mean position of $P$ in the sense that $\boldsymbol{z}$ is the average of all the positions that $P$ occupies over the period of excitation surrounding the time instant $t$ . It is this idea that is then developed into the notions like Stokes drift and Lagrangian-mean motion. To be explicit, as discussed in, say, Andrews & McIntyre (Reference Andrews and Mcintyre1978) and Bühler (Reference Bühler2009), a particle with nominal position $\boldsymbol{z}$ at time $t$ has true position at $\boldsymbol{z}+\boldsymbol{q}(\boldsymbol{z},t)$ . This is formalized as a map (cf. e.g. (2.2) in Andrews & McIntyre (Reference Andrews and Mcintyre1978)) as follows:

where $\boldsymbol{q}$ is referred to as the ‘disturbance-associated particle displacement field’ in Andrews & McIntyre (Reference Andrews and Mcintyre1978), and as the ‘lift field’ in Bühler (Reference Bühler2009), and where $\boldsymbol{q}(\boldsymbol{z},t)$ is required to have the following property:

This map is then used to define the concept of Lagrangian-mean average. In this paper, the Lagrangian mean of a field $\unicode[STIX]{x1D719}(\boldsymbol{z},t)$ will be denoted by $\unicode[STIX]{x1D719}^{L}(\boldsymbol{z},t)$ and defined as in Andrews & McIntyre (Reference Andrews and Mcintyre1978), Bühler (Reference Bühler2009), i.e.

Applying (2.46) and (2.47)

so that the nominal position $\boldsymbol{z}$ is (by definition) the Lagrangian mean of the true position. Focusing now on our kinematics, we recall it is the position $\boldsymbol{y}$ that identifies points in the current configuration. Hence, using (2.2) and (2.7), we have

Therefore, the kinematic framework we have adopted is simply that of the standard GLM theory, where $\boldsymbol{x}$ is the Lagrangian-mean position and $\unicode[STIX]{x1D743}(\boldsymbol{x},t)$ is the lift field. This is not surprising, in that we have adopted an ALE perspective and the ALE mapping from the current domain to a convenient (e.g. computationally) domain $\mathscr{B}(T)$ is exactly the map in (2.45) in the GLM theory. That is, at least when it comes to the kinematics of position, the GLM theory is an ALE scheme. At this point it is important to remark that there is no difficulty in reconciling that the domain $\mathscr{B}(T)$ is time-dependent whereas $\unicode[STIX]{x1D6FA}$ is not, even in the GLM theory. As already remarked, $\unicode[STIX]{x1D6FA}$ is time-independent because it is a control volume. By contrast $\mathscr{B}(T)$ is a configuration and as such it is time-dependent because material particles have moved to go from $B_{0}$ to $\mathscr{B}(T)$ . However, if the particles on the boundary of $\mathscr{B}(T)$ do not move, i.e. if $\boldsymbol{u}^{(0)}(\boldsymbol{X},T)=\mathbf{0}$ for $\boldsymbol{X}\in \unicode[STIX]{x2202}B_{0}$ , then $\mathscr{B}(T)$ will always occupy the same portion of $\mathscr{E}^{d}$ as $B_{0}$ , and it is this portion that coincides with $\unicode[STIX]{x1D6FA}$ . In other words, $\unicode[STIX]{x1D6FA}$ and $\mathscr{B}(T)$ coincide if they are both viewed as mere (geometric) point sets. We contend that, in the study of acoustofluidic devices, this boundary condition is a reasonable and convenient one to assume, namely, that the motion of the system, on average, does not displace the boundary of the fluid compartment of the device. Hence, in the remainder of the paper we formally assume that

We now turn to comparing the velocity descriptors. We start by observing that the velocity field $\boldsymbol{v}(\boldsymbol{y},t)$ does not find correspondence in the standard approaches because its domain is $B(t)$ , which is a time-dependent deformable domain, as opposed to the control volume $\unicode[STIX]{x1D6FA}$ . However, as we have already remarked that our kinematic framework is that of the GLM theory, we should be able to find correspondences between said theory and our approach. To this end, let us recall that

Hence, from (2.16), and recalling that $\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D743}$ is harmonic and $\unicode[STIX]{x1D716}^{\unicode[STIX]{x1D6FD}}\unicode[STIX]{x2202}_{T}\boldsymbol{u}^{(0)}=\unicode[STIX]{x2202}_{t}\boldsymbol{u}^{(0)}$ , we have

Referring to (2.42), we also see that

Both (2.52) and (2.53) are results worth discussing. Looking ahead at the numerical implementation of our scheme, we will select $\langle \boldsymbol{v}^{\ast (2)}\rangle$ as one of the primary unknowns of our time-averaged second-order problem. Hence, equation (2.53) implies that the Lagrangian-mean velocity is $O(\unicode[STIX]{x1D716}^{2})$ , and it is precisely one of the fields that our solution scheme delivers directly.

Another issue concerning our kinematic choices concerns the field $\boldsymbol{u}^{(0)}$ . In a fluids problem in which there is no need to determine the displacement field, one can simply not solve for it. However, being able to solve for the field $\boldsymbol{u}^{(0)}$ , whether as a field with domain $B_{0}$ or as a field $\boldsymbol{u}^{(0)\#}$ defined over $\mathscr{B}(T)$ (see definition in (2.23)), is essential when interested in modelling the streaming flows of complex fluids that have an elastic component to their constitutive response and/or to problems in which there are deformable solid objects immersed in the fluid. Hence, we discuss here how one could pose the problem of determining such a field. Clearly, any such discussion needs to start from considering carefully (2.52). If the displacement field is described through the function $\boldsymbol{u}^{(0)\#}$ , then (2.52) must be rewritten as

Whether in using (2.52) or (2.54), the presence of a time-averaging operator makes the relationship between the primary unknown $\boldsymbol{v}^{L}$ and the displacement a nonlinear integro-differential equation. In those cases where considerations based on the smallness parameter $\unicode[STIX]{x1D716}$ allow it, it is possible to consider an approximation that reduces to a partial differential equation. To illustrate this point, we will focus on (2.52). From the discussion of (2.42), we know that $\unicode[STIX]{x2202}_{t}\boldsymbol{u}^{(0)}=O(\unicode[STIX]{x1D716}^{\unicode[STIX]{x1D6FD}})$ . Since $\unicode[STIX]{x1D641}_{\unicode[STIX]{x1D709}}\unicode[STIX]{x2202}_{t}\boldsymbol{u}^{(0)}=\unicode[STIX]{x2202}_{t}\boldsymbol{u}^{(0)}+\unicode[STIX]{x1D735}_{x}\unicode[STIX]{x1D743}\unicode[STIX]{x2202}_{t}\boldsymbol{u}^{(0)}$ , in view of the third of (2.37), we can then write

This result implies that if an approximation up to order $\unicode[STIX]{x1D6FD}+1$ is acceptable, one can reconstruct the displacement field directly from the above equation. We note that, from a numerical viewpoint, FSI schemes called immersed methods (Peskin Reference Peskin2002; Wang & Liu Reference Wang and Liu2004; Heltai & Costanzo Reference Heltai and Costanzo2012) are able to reconstruct the displacement field from an equation such as (2.55). Higher-order corrections that take into account the time-averaging operations are possible and, in principle, can also be handled via immersed methods.

Before going to the next section, we observe that (2.55), along with (2.50) implies that

3.1 Remarks on the balance of mass

The character of the boundary value problems induced by the selected orders of magnitude is somewhat different, to the best of the authors’ knowledge, from what is currently found in the literature. To elucidate this point, we begin with considerations based on the balance of mass in Lagrangian form reported in (2.27). Using (2.14) and (2.43), we can rewrite the balance of mass as

Expanding, using (2.15), and matching orders we obtain the following relations:

For the case with $\unicode[STIX]{x1D70C}_{0}$ a constant, as is typically the case for fluids, we can then take the partial derivative of the above relations with respect to time while holding $\boldsymbol{x}$ fixed to obtain

where we have used the fact that $\unicode[STIX]{x2202}_{t}\det (\unicode[STIX]{x1D641}^{(0)\#})=O(\unicode[STIX]{x1D716}^{\unicode[STIX]{x1D6FD}})$ . Of the above relationships, perhaps the most important are (3.4) and (3.7), because they indicate that the proposed kinematic setting implies that second-order mass density distribution is completely determined by the first-order solution to the problem. Furthermore, recalling that the time average of harmonic functions over the excitation period is time-independent, from (3.4) we see that $\langle \unicode[STIX]{x1D70C}^{\ast (2)}\rangle$ is no longer a function of time. In turn this implies that one can calculate the steady (slow-time) value of $\langle \unicode[STIX]{x1D70C}^{\ast (2)}\rangle$ from the first-order solution. Later we will show that (3.7) allows us to obtain a corresponding result for the second-order velocity solution such that it is possible to solve directly for the steady value of $\langle \boldsymbol{v}^{\ast (2)}\rangle$ .

3.2 Zeroth-order subproblem

We now derive the zeroth-order problem by substituting the expansions of velocity and mass density and retaining only zeroth-order terms. Recalling that we have neglected the background flow, the zeroth-order subproblem is governed by the following relations:

Observing that the first of (3.8) implies that $\unicode[STIX]{x1D70C}^{\ast (0)}$ is a constant, for constant $\boldsymbol{b}^{\ast }$ , we see that the solution to the zeroth-order problem is represented by the static equilibrium state of the fluid. For convenience, we will proceed forward under the assumption that the body forces are negligible. In this case, the solution of the zeroth-order problem is

The first of (3.9) along with (3.2) imply that

which, in turn, implies that the zeroth-order motion is expected to be volume preserving under the current assumptions.

3.3 Acoustic subproblem

Next, we consider the acoustic subproblem (on the fast time scale) by substituting the expressions for velocity and mass density in terms of their components in the balance laws and retaining only the first-order terms. Using (2.24) and (2.25) along with the expansions of velocity and mass density, we obtain the following set of equations for the acoustic subproblem:

where

We notice that the first of (3.11) is a consequence of (3.3), the latter providing an expression of $\unicode[STIX]{x1D70C}^{\ast (1)}$ in terms of the displacement $\unicode[STIX]{x1D743}$ . In fact, recalling that $\unicode[STIX]{x1D70C}^{\ast (0)}=\unicode[STIX]{x1D70C}_{0}$ , and along with (2.17) and (2.41), the first-order problem takes on the following form:

For the above equation we only seek harmonic solutions. If we then let the solution be of the form $\unicode[STIX]{x1D743}(\boldsymbol{x},t)=\text{e}^{\text{i}\unicode[STIX]{x1D714}t}\hat{\unicode[STIX]{x1D743}}(\boldsymbol{x})$ , with $\hat{\unicode[STIX]{x1D743}}$ a complex-valued vector field, the first-order problem reduces to a linear elasto-statics problem with complex moduli.

As far as the boundary conditions are concerned, there are several sets of interest in applications. In fact, the boundary problem of interest is a coupled problem in which the acoustics is induced by an electrical actuator for which the electric field is prescribed and for which the electric actuation produces a corresponding mechanical actuation that depends on the geometry and material make-up of the device. In this paper we do not consider the fully coupled electro-mechanical-fluidic problem. Instead we simply assume that the displacement field $\unicode[STIX]{x1D743}$ is prescribed everywhere on the boundary:

where the subscript $g$ stands for given.

3.4 Mean dynamics subproblem

Proceeding in a similar manner as done above, we derive the second-order problem. Using (2.24) and (2.25) along with the expansions of velocity and mass density, for the balance of mass and momentum we obtain, respectively,

and

where

As can be noted, the forcing terms in the second-order problem show a prominent role of the fluid’s viscous response in the first-order solution. This may appear to be a departure from other treatments of the subject, especially when the governing equations are kept in Eulerian form as, for example, in (2.6a) and (2.6b) in Köster (Reference Köster2006). In reality, there is no such departure. In fact, the forcing terms in the second-order that appear inertia-related in some treatments, can be completely rewritten in terms of the first-order stress response. The latter makes apparent the role of viscosity in the forcing in question. This issue has also been pointed out by Vanneste & Bühler (Reference Vanneste and Bühler2011) (cf. second paragraph on p. 1780, as well as (4.4)).

We now observe that the term $\unicode[STIX]{x2202}_{t}\boldsymbol{v}^{\ast (2)}$ in (3.16) is of order higher than two. In fact, using (2.42), we have

As $\unicode[STIX]{x1D735}_{x}\boldsymbol{v}_{\unicode[STIX]{x1D709}}=O(\unicode[STIX]{x1D716})$ , $\unicode[STIX]{x2202}_{t}\boldsymbol{u}^{(0)}=O(\unicode[STIX]{x1D716}^{2})$ , and $\unicode[STIX]{x2202}_{tt}\boldsymbol{u}^{(0)}=O(\unicode[STIX]{x1D716}^{4})$ , the above result implies that $\unicode[STIX]{x2202}_{t}\boldsymbol{v}^{\ast (2)}=O(\unicode[STIX]{x1D716}^{3})$ . The latter result, combined with (3.7), allows us to rewrite (3.15)–(3.17) as follows:

The first of (3.19) implies that $\boldsymbol{v}^{\ast (2)}$ must be divergence-free. In turn, this means that $\unicode[STIX]{x1D64B}^{\ast (2)}$ is determined by (3.17) only up to a hydrostatic contribution needed to accommodate the kinematic constraint in the first of (3.19) (Gurtin et al. Reference Gurtin, Fried and Anand2010). Therefore, recalling that $\unicode[STIX]{x1D70C}^{\ast (0)}=\unicode[STIX]{x1D70C}_{0}$ and (3.3), the expression for $\unicode[STIX]{x1D64B}^{\ast (2)}$ that we will use in the formulation of the second-order problem will be as follows:

where $q(\boldsymbol{x},t)$ is a scalar Lagrange multiplier that is determined by enforcing the constraint in the first of (3.19). We should remark that, from a strictly analytical viewpoint, the above relationship is equivalent to

where $\hat{q}$ is the sum total of $q$ and all other hydrostatic contributions. However, depending on the numerical method used to solve the problem, and on the corresponding discretization strategy, the two expressions are not equivalent from a numerical viewpoint. As discussed later, our numerical scheme is a mixed finite element method (FEM) formulation implemented via traditional Lagrange polynomials for the velocity field and the Lagrange multiplier $q$ . For said numerical approach, we have found that (3.20) leads to results without spurious numerical oscillations.

We now proceed to time average the equations developed so far. The absence of time derivative in the governing equations of the second-order problem makes the time-averaging operation straightforward and transparent. In particular, we recall that the time-average operation commutes with the spatial differentiation operators. Furthermore, invoking (2.53), we recall that the time averaging of $\boldsymbol{v}^{\ast (2)}$ yields the Lagrangian-mean velocity. Finally, as argued in § 2.4, we recall that the time-averaging operation is equivalent to a change of temporal scale from the fast to the slow time scale. Hence, upon time-averaging, the time-averaged second-order problem is governed by the following equations:

with

where $\boldsymbol{v}^{\ast (1)}=\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D743}$ , and therefore the last three lines of the above expressions are reduced to a single second-order tensor known from the first-order solution. As far as the boundary conditions are concerned, these have already been stated in (2.56). The main observation concerning the present formulation is that, granted the kinematic framework we have laid out, the boundary conditions for the second-order problem are exact, as opposed to approximate, as they appear in the purely Eulerian formulations of the second-order problem (see, e.g., Köster Reference Köster2007).

The formulation of the second-order problem in the proposed ALE framework has yielded a (slow) time-independent problem at the second order. If there are no bodies immersed in the fluid and convected by it, then the second-order problem is intrinsically steady. This is in contrast to the Eulerian formulations found in the literature (see, e.g., Köster Reference Köster2007; Muller et al. Reference Muller, Barnkob, Jensen and Bruus2012), which first derive a time-dependent second-order problem even upon time-averaging, and then consider an associated steady problem by setting the time derivatives of the primary unknowns to zero. It should also be remarked that our result is consistent with the analysis by Xie & Vanneste (Reference Xie and Vanneste2014b ), which, in a regime like that considered here, also derives an intrinsically time-independent second-order problem. This is an important observation in light of the recent studies which consider a time-dependent second-order case using Eulerian formulation (Muller & Bruus Reference Muller and Bruus2015). We believe the presence of the time-dependent terms at the second order in Eulerian formulation to be suspicious in that it is derived by a time-averaging operation that is somewhat ambiguous as to how it transitions from a time scale to another. When the time-scale separation is explicitly defined, as was done here and in Xie & Vanneste (Reference Xie and Vanneste2014b ), the analysis of the motion leads to the conclusion that the second-order problem is not time-dependent at the slow time scale. This is also consistent with the fact that, whether in our formulation or in purely Eulerian ones, the forcing terms appearing in the equations, along with the boundary conditions, are time-independent.

If an object is immersed in the fluid, the corresponding streaming flow, in general, would not be steady. Provided that a more careful analysis of this case is needed, the present formulation suggests that one can proceed to update the position of the immersed object on the slow time by integrating the Lagrangian-mean velocity over time intervals smaller than $\unicode[STIX]{x1D6F1}/\unicode[STIX]{x1D716}^{2}$ . In the results section of the paper, we have included an example of the tracking of the displacement of a portion of the fluid using the immersed method discussed in Heltai & Costanzo (Reference Heltai and Costanzo2012). Applications with immersed solid objects in a streaming flow are currently being developed and will be presented in future publications.

As a final comment on the governing equations of the second-order problem, we note that the first of (3.22) is consistent with (3.10), the latter pertaining to the zeroth-order problem. That is, if we view the Lagrangian-mean velocity as the effective velocity field of the motion at the slow time scale, requiring that the velocity field in question is divergence-free is an expression of the fact that the slow motion is volume preserving, as was argued after (3.10).

5.1 Rectilinear actuation

This numerical test is motivated by the bulk acoustic wave (BAW) devices considered by Muller et al. (Reference Muller, Barnkob, Jensen and Bruus2012). With reference to figure 2(a), the computational domain consists of a rectangular 2D fluid-filled cavity with length $L=380~\unicode[STIX]{x03BC}\text{m}$ and height $H=160~\unicode[STIX]{x03BC}\text{m}$ , where the left and the right walls are actuated harmonically along the $x$ -direction. The actuation frequency considered is $f=1.97~\text{MHz}$ , which corresponds to a half-wave resonance mode for the channel dimensions considered. Using a Cartesian coordinate system with origin at the lower left corner of the domain, the actuation displacement at the oscillating walls is given by

where $u_{x}$ and $u_{y}$ are the displacement components in the $x$ and $y$ directions, respectively, $d_{0}=1\times 10^{-10}~\text{m}$ is the actuation amplitude, and $\unicode[STIX]{x1D714}=2\unicode[STIX]{x03C0}f$ is the angular frequency. The fluid properties are taken to be those of water and are the same as those considered by Muller et al. (Reference Muller, Barnkob, Jensen and Bruus2012).

Figure 2.

(a) ALE formulation: density plot of the magnitude of the Lagrangian-mean velocity for a bulk acoustic wave device, where the vertical walls are actuated as indicated in (5.1)–(5.4). The red lines represent the streamlines of the Lagrangian-mean velocity. (b) ALE formulation: plot of the $x$ component of different velocity fields away from the left wall (at $x=0$ ) along a horizontal line at $y=100~\unicode[STIX]{x03BC}\text{m}$ . At the oscillating wall (i.e.  $x=0$ ), the Eulerian velocity is equal and opposite to the Stokes drift, thereby yielding a zero Lagrangian velocity at the oscillating wall, representing no flow across the wall.

The results in figure 2 were obtained with the ALE method. Figure 2(a) shows a density plot of the magnitude of Lagrangian-mean velocity $\boldsymbol{v}^{L}$ . The lines superimposed to the density plot trace the streamlines of $\boldsymbol{v}^{L}$ . Our results indicate the presence of four vortices inside the channel, similar to those reported by Muller et al. (Reference Muller, Barnkob, Jensen and Bruus2012) via a $E_{zbc}$ formulation. In fact our results are very similar to those in Muller et al. (Reference Muller, Barnkob, Jensen and Bruus2012) in the bulk of the channel. However, they differ significantly near the oscillating walls. Before illustrating the differences, it is important to keep in mind that differences and similarities of results need to be understood in context, i.e. with an awareness of the size of the device and the forcing frequency and amplitude. With this in mind, it is generally accepted that a second-order Eulerian velocity field does not accurately represent particle trajectories. Accurate predictions of these are obtained by summing to the Eulerian velocity field the Stokes drift (Bühler Reference Bühler2009). The fact that for the case under discussion our results match those in Muller et al. (Reference Muller, Barnkob, Jensen and Bruus2012), the latter obtained without Stokes drift correction, indicates that in this particular problem the correction in question may be negligible in the bulk of the channel. To verify this conjecture we show in figure 2(b) the plot of the $x$ component of Lagrangian-mean velocity, the Eulerian velocity, and the Stokes drift near the oscillating left wall. Again, the plots shown in this figure were obtained via the ALE formulation. It can be observed that the Eulerian velocity and Stokes drift are equal in magnitude but opposite in direction at the oscillating wall, thereby giving a zero Lagrangian velocity at the oscillating wall. We also note that while the Stokes drift is equal and opposite to the Eulerian velocity field at the oscillating boundaries, the Stokes drift decays significantly away from the oscillating walls and the maximum value of Stokes drift is observed to be one order of magnitude smaller than the maximum value of the Eulerian velocity field in the bulk. Therefore, the contribution of Stokes drift to the net Lagrangian field is indeed negligible in the bulk of the channel, resulting in no significant differences between the Eulerian and the Lagrangian velocity fields in the bulk of the channel. Thus, the results obtained from the ALE formulation do not show any significant differences in the bulk of the channel compared to those reported by Muller et al. (Reference Muller, Barnkob, Jensen and Bruus2012) for this specific case.

What is remarkable about this conclusion is that, in reality, the boundary conditions in the $E_{zbc}$ formulation are known to pose difficulties in their physical meaning. In fact, owing to the rectilinear nature of the actuation displacement profile, there is no slip velocity along the oscillating walls and only a non-zero transverse component of Stokes drift exists at the oscillating walls. From a physical perspective, a zero Eulerian velocity at the oscillating walls along with a non-zero Stokes drift correction leads to the conclusion that fluid mass is being convected through the oscillating walls. By contrast, this contradiction is completely avoided both in the current ALE formulation and in the $E_{mbc}$ formulation. In the latter, the boundary value of the Eulerian velocity is not zero, but it is equal and opposite to the Stokes drift, thereby producing a zero Lagrangian-mean velocity. It is important to note that, from a numerical viewpoint, the ALE and the $E_{mbc}$ formulations do not yield identical results. This aspect will be considered in greater detail in § 5.3.

5.2 Non-rectilinear actuation

In the previous section we mentioned the fact that, in some cases, the Eulerian velocity can be taken as an accurate enough descriptor of the mean particle velocity. In this section we present a more in-depth analysis of this question and we will point out that there are a number of applications of technological relevance in which the Eulerian velocity is not a good enough descriptor of particle motion, especially as the size of the device becomes of the same order of magnitude as or smaller than the actuation wavelength.

5.3 Remarks on the numerical implementation

While the theoretical framework for our ALE formulation may appear more involved that the Eulerian setting of previous formulations, it is important to note that the boundary value problems actually solved both for the first- and the second-order problems are simpler that those of said formulations. The solution of the first-order problem requires the solution of a linear problem in the displacement field $\unicode[STIX]{x1D743}$ with no pressure or density constraints and the second-order problem has the structure of a steady Stokes problem. In addition to these formal properties, the ALE formulation provides the advantage of identifying the Lagrangian-mean velocity field directly, i.e. without a post-processing step that is susceptible to numerical aberration. To illustrate this point, we consider the oscillating sharp-edge device considered by Nama et al. (Reference Nama, Huang, Huang and Costanzo2014, Reference Nama, Huang, Huang and Costanzo2016) and by Huang et al. (Reference Huang, Xie, Ahmed, Rufo, Nama, Chen, Chan and Huang2013, Reference Huang, Nama, Mao, Li, Rufo, Chen, Xie, Wei, Wang and Huang2014, Reference Huang, Chan, Li, Nama, Xie, Wei, Chen, Ahmed and Huang2015a ,Reference Huang, Ren, Nama, Li, Li, Yao, Cuento, Wei, Chen and Xie b ). Figure 6(a) shows the schematic of a typical sharp-edge-based acoustofluidic device consisting of protruded sharp edges along the channel walls in a periodic fashion. Due to the inherent periodicity of the geometry as well as the actuation, we consider a periodic cell of the device as shown in figure 6(b). In our calculation, we have set $L=600~\unicode[STIX]{x03BC}\text{m}$ , $\unicode[STIX]{x1D6FC}=15^{\circ }$ , $h=200~\unicode[STIX]{x03BC}\text{m}$ and $H=600~\unicode[STIX]{x03BC}\text{m}$ . The top and the bottom channel walls (including the sharp edges) denoted by $\unicode[STIX]{x1D6E4}_{t}$ and $\unicode[STIX]{x1D6E4}_{b}$ , respectively, are assumed to undergo uniform harmonic motion along the $x$ direction with an actuation frequency $f=5~\text{kHz}$ . The side walls of the channel, denoted by $\unicode[STIX]{x1D6E4}_{l}$ and $\unicode[STIX]{x1D6E4}_{r}$ , respectively, are subjected to periodic boundary conditions. To avoid the singularity associated with the second-order velocity field at the sharp-edge tip, we have rounded off the tip of the sharp edges with a small radius of $2~\unicode[STIX]{x03BC}\text{m}$ . Figure 7(a) shows the plot of the $y$ component of the Lagrangian and the Eulerian velocities along a line parallel to $y$ axis moving away from the sharp-edge tip towards the top wall. Because the actuation consists of an oscillatory displacement purely parallel to the $x$ direction, the $y$ component of the particle velocity at the sharp edge is equal to zero. The ALE formulation, in which this boundary condition is applied directly, displays a mean Lagrangian velocity field consistent with the boundary condition in question. However, what is more important, is that the component of $\boldsymbol{v}^{L}$ obtained via the ALE formulation is not characterized by numerical oscillations spanning individual elements in the simulation, as is the case for the corresponding component of the Eulerian velocity. Figure 7(b) displays the same result, this time obtained via the $E_{zbc}$ formulation. As expected, here we obtained the reverse situation, that is, given a relatively smooth result for the Eulerian velocity, the reconstructed Lagrangian-mean velocity suffers from severe numerical oscillations, often spanning individual elements. The reason behind these oscillations in question is that the reconstruction of, say, the Lagrangian-mean velocity from the Eulerian velocity requires the determination of Stokes drift, which in turn contains spatial gradients of the first-order velocity solution, the latter, for harmonic solutions, being a close relative to $\unicode[STIX]{x1D735}_{x}\unicode[STIX]{x1D743}$ . As the quality of the numerical results decreases (in general) with the order of differentiation, the choice of the approximation spaces for the solution then becomes crucial in achieving a desired accuracy in the post-processing operation at hand. Consequently, if a simulation is meant to focus on particle velocities and trajectories, there is a definite advantage in choosing a formulation in which the Lagrangian-mean velocity appears as one of the primary unknowns of the problem. The reconstruction of the Lagrangian-mean velocity from the Eulerian velocity field remains a delicate operation even when adopting the $E_{mbc}$ formulation. To illustrate this point, we present in figure 8(a) a plot of the $y$ component of the Lagrangian-mean velocity obtained from our ALE formulation and with the $E_{mbc}$ formulation for the simulation discussed in § 5.1. Specifically, figure 8(a) presents the plot of $v_{y}^{L}$ along a horizontal line (shown in the inset as a black line) at $y=H/3$ (cf. figure 2 a). In figure 8(b) we report a result for the device shown in the inset, but this time pertaining to the $x$ component of $\boldsymbol{v}^{L}$ as obtained via our ALE formulation and the $E_{mbc}$ formulation. In some sense, the case in figure 8(b) is not surprising, as the numerical aberration under discussion occurs in the proximity of a boundary with a re-entrant corner geometry. That is, the reconstruction of $\boldsymbol{v}^{L}$ from the Eulerian velocity as delivered by the $E_{mbc}$ formulation is expected to be more susceptible to numerical aberration in the neighbourhood of a high-velocity-gradient location. What is more surprising is the result in figure 8(a). In this case the result coming from the $E_{mbc}$ formulation is expected to be essentially identical to that coming from the ALE formulation given that the Stokes drift is negligible in the bulk of the BAW device under consideration. So far we have outlined the more apparent numerical issues behind the difficulty of an accurate reconstruction of $\boldsymbol{v}^{L}$ from the needed Eulerian fields. In reality, at the root of these difficulties there is a more subtle analytical problem surrounding the determination of the inverse of the ALE map $\boldsymbol{x}\mapsto \boldsymbol{x}+\unicode[STIX]{x1D743}(\boldsymbol{x},t)$ . From a theory of partial differential equations perspective, this issue is too technical to be addressed in this publication, but it has been mentioned by Bühler in the Lagrangian mean field theory section of Bühler (Reference Bühler2009). To clarify somewhat the issue at hand, one can consider the derivation of the Lagrangian-mean velocity from within an Eulerian formulation, which we have reported in the appendix in (A 9)–(A 14). With reference to (A 10) in particular, and consequently to (A 13), we see that the concept of Lagrangian-mean velocity presupposes the existence of a smooth enough inverse of the aforementioned ALE map. It is not clear under what conditions one can guarantee the existence of this map for long enough time intervals. With this in mind, and granted that this analytical problem is not answered by simply switching to an ALE perspective as we have done, the approach proposed herein in which $\boldsymbol{v}^{L}$ is a primary unknown (as opposed to a quantity computed in a post-processing step) does offer some advantages from a numerical viewpoint for simulations that focus on predicting particle trajectories, as is the case in the design of devices for the accurate placement of particles in the flow and for validation based on particle velocimetry experiments.

Figure 6.

Schematic design of a typical oscillating sharp-edge-based device.

Figure 7.

Simulations pertaining to the sharp-edge device in figure 6: $y$ component of the Lagrangian-mean velocity of the and Eulerian velocity along a vertical line emanating from a sharp-edge tip obtained using (a) ALE formulation and (b) $E_{zbc}$ formulation.

Figure 8.

(a) Plot of the $y$ component of the Lagrangian-mean velocity obtained from the ALE formulation and the $E_{mbc}$ formulation along a horizontal line (shown in the inset as a black line) at $y=H/3$ for the bulk acoustic wave device considered in § 5.1. (b) Plot of the $x$ component of the Lagrangian-mean velocity obtained from the ALE formulation and $E_{mbc}$ formulation along a line parallel to the $y$ emanating from the sharp-edge tip towards the top wall as shown in the inset. As it computes $\boldsymbol{v}^{L}$ directly, i.e. avoiding a postoprocessing involving the gradients of the first-order solution, the ALE formulation yields smoother numerical results in both cases.

5.4 Fluid body tracking

The last set of numerical results we show concerns the tracking of a region of fluid within a streaming field. This simulation is being reported here as a ‘proof of concept’ in that it demonstrates one major application of the proposed ALE formulation in the determination of the motion of solid bodies immersed in a streaming flow. This is a problem that has been considered extensively in more traditional FSI contexts (see, e.g., Peskin Reference Peskin2002; Zhang & Gay Reference Zhang and Gay2007; Heltai & Costanzo Reference Heltai and Costanzo2012), but it has not been given a satisfactory treatment in microacoustofluidic applications. The latter are of strong interest in acoustophoresis for controlling the motion and placement of immersed particles, whether soft or hard, and in applications for the manipulation of relatively large soft immersed objects such as cells. Here we limit ourselves to illustrating the possibility of tracking the motion of an immersed body by tracking the motion of a connected fluid sub-body (an extensive analysis of the motion of soft immersed solid bodies by the authors is in preparation). Figure 9 shows several frames of a simulation in which we track the motion of the fluid sub-body originally occupying the circular domain shown in figure 9(a). The simulation was obtained by using the immersed method discussed in Heltai & Costanzo (Reference Heltai and Costanzo2012) for the integration of (2.55) (including only terms up to order two). The density plot over the tracked body displays the magnitude of the displacement. The simulation in question is akin to determining the trajectory of a passive tracer where the tracer in question is not a point but a continuous region of fluid. If the sub-body in question consisted of, say, a soft elastic object, the computed displacement $\boldsymbol{u}^{(0)}$ would then allow the determination of the elastic stresses induced by the deformation and the consequent feedback on the motion of the surrounding fluid.

Figure 9.

Simulation showing the deformation of a continuous fluid sub-body originally occupying the circular region in frame (a).