# Theory of pressure acoustics with thermoviscous boundary layers and
streaming in elastic cavities

Jonas Helboe Joergensen<EMAIL_ADDRESS>Department of Physics, Technical
University of Denmark,
DTU Physics Building 309, DK-2800 Kongens Lyngby, Denmark  Henrik Bruus
<EMAIL_ADDRESS>Department of Physics, Technical University of Denmark,
DTU Physics Building 309, DK-2800 Kongens Lyngby, Denmark

(14 December 2020)

###### Abstract

We present an effective thermoviscous theory of acoustofluidics including
pressure acoustics, thermoviscous boundary layers, and streaming for fluids
embedded in elastic cavities. By including thermal fields, we thus extend the
effective viscous theory by Bach and Bruus, J. Acoust. Soc. Am. 144, 766
(2018).(Bach2018, ) The acoustic temperature field and the thermoviscous
boundary layers are incorporated analytically as effective boundary conditions
and time-averaged body forces on the thermoacoustic bulk fields. Because it
avoids resolving the thin boundary layers, the effective model allows for
numerical simulation of both thermoviscous acoustic and time-averaged fields
in 3D models of acoustofluidic systems. We show how the acoustic streaming
depends strongly on steady and oscillating thermal fields through the
temperature dependency of the material parameters, in particular the viscosity
and the compressibility, affecting both the boundary conditions and spawning
additional body forces in the bulk. We also show how even small steady
temperature gradients ($\sim 1\,\textrm{K}/\textrm{mm}$) induce gradients in
compressibility and density that may result in very high streaming velocities
($\sim 1\,\textrm{mm}/\textrm{s}$) for moderate acoustic energy densities
($\sim 100\,\textrm{J}/\textrm{m}^{3}$).

##  I Introduction

Modeling and simulation is important for designing microscale acoustofluidic
systems. Traditionally, most models have been purely mechanical, but some
include thermal effects, such as in the studies of the acoustic radiation
force acting on suspended microparticles (Doinikov1997, ; Danilov2000, ;
Karlsen2015, ) and of acoustic streaming in rigid cavities. (Rednikov2011, ;
Muller2014, )

Here, we focus on acoustic streaming, where recent developments in the field
points to the necessity of making a full thermoviscous analysis. Karlsen et
al. introduced the acoustic body force acting on a liquid governed by solute-
induced gradients in the compressibility and density of the
liquid.(Karlsen2016, ) This force has explained the iso-acoustic focusing of
mircoparticles, (Augustsson2016, ) patterning of concentration profiles,
(Karlsen2017, ) and suppression of acoustic streaming. (Karlsen2018, ;
Qiu2019, ) Simultaneously, Bach and Bruus developed the effective theory for
pressure acoustics and streaming in elastic cavities,(Bach2018, ) in which the
viscous boundary layer was solved analytically and imposed as an effective
boundary condition to the bulk field. This model has enabled simulations of
cm-sized three-dimensional (3D) acoustofluidic systems, (Skov2019, ;
Skov2019b, ) with hitherto prohibitive computational costs, and it has
provided a deeper insight in the physics of boundary- and bulk-induced
streaming, but without thermal effects.(Bach2019, )

In this work, we combine our previous work on thermoviscous streaming in rigid
systems, (Muller2014, ) thermoviscous potential theory, (Karlsen2015, ) the
theory of pressure acoustics with viscous boundary layers and streaming in
curved elastic cavities, (Bach2018, ) and the 3D numerical modeling of
acoustofluidic systems using the latter theory, (Skov2019, ) and develop an
effective thermoviscous theory for a fluid-filled cavity embedded in an
elastic solid. The theory includes both steady and acoustic temperature fields
for pressure acoustics with thermoviscous boundary layers and for streaming
with thermoviscous body forces. In Section II we set up the basic theory and
model assumptions. In Sections III-V, the governing equations and boundary
conditions are derived from the theory for the zeroth, first, and second order
in the acoustic perturbation, respectively. In Section VI, the theory is
implemented in a numerical model, which is then used in two examples to show
the nature and importance of thermal effects in acoustofluidics. Finally, we
conclude in Section VII.

##  II Basic theory and model assumptions

We consider an acoustofluidic device consisting of an elastic solid containing
a microchannel filled with a thermoviscous Newtonian fluid and actuated by an
piezoelectric transducer at a single frequency in the MHz range. This time-
harmonic actuation establishes an acoustic field in the system, which in the
fluid, by the internal dissipation and hydrodynamic nonlinearities, results in
a time-averaged response that leads to acoustic streaming.

###  II.1 Governing equations

In this work, unlike prior work,(Skov2019, ) we leave the piezoelectric
transducer out of the analysis, and only represent it by an oscillating
displacement condition on part of the surface of the elastic solid. The
response of the fluid embedded in the elastic solid to this oscillating-
displacement boundary condition is controlled by the hydro-, elasto-, and
thermodynamic governing equations of the coupled thermoviscous fluid and
elastic solid.

The linear elastic solid is described in the Lagrangian picture by the fields
of the density $\rho$, the displacement $\bm{u}$, and the temperature $T$.
There are also associated eight material parameters: the longitudinal and
transverse sound speeds $c_{\mathrm{lo}}$ and $c_{\mathrm{tr}}$, the thermal
conductivity $k^{\mathrm{th}}$, the specific heat $c_{p}$, the ratio of
specific heats $\gamma=c_{p}/c_{v}$, the thermal expansion coefficient
${\alpha_{p}}$, and the isentropic and isothermal compressibilities
$\kappa_{s}$ and $\kappa_{T}$. The velocity field is given as the time
derivative of the displacement field
$\bm{v}^{\mathrm{sl}}=\partial_{t}\bm{u}$, so no advection occur, and the
governing equations are the transport equations of the momentum density
$\rho\partial_{t}\bm{u}$ and temperature $T$,(Landau1986, ; Karlsen2015, )

$\displaystyle\rho\partial^{\,2}_{t}\bm{u}$
$\displaystyle=\bm{\nabla}\cdot\bm{\sigma},$ (1a)
$\displaystyle\partial_{t}T+\frac{(\gamma-1)}{{\alpha_{p}}}\partial_{t}(\bm{\nabla}\cdot\bm{u})$
$\displaystyle=\frac{\gamma}{\rho
c_{p}}\bm{\nabla}\cdot(k^{\mathrm{th}}\bm{\nabla}T),$ (1b)

where $\bm{\sigma}$ is the stress tensor, which for isotropic solids is,

$\displaystyle\bm{\sigma}$
$\displaystyle=-\frac{{\alpha_{p}}}{\kappa_{T}}(T-T_{0})\bm{I}+\bm{\tau},$
(2a) $\displaystyle\bm{\tau}$ $\displaystyle=\rho
c^{2}_{\mathrm{tr}}\Big{[}\bm{\nabla}\bm{u}+(\bm{\nabla}\bm{u})^{\dagger}\Big{]}+\rho\big{(}c^{2}_{\mathrm{lo}}-2c^{2}_{\mathrm{tr}}\big{)}(\bm{\nabla}\cdot\bm{u})\bm{I}.$
(2b)

The fluid is described in the Eulerian picture by the fields of the density
$\rho$, the pressure $p$, the velocity $\bm{v}$, the temperature $T$, and the
energy per mass unit $\epsilon$, and by eight material parameters: the dynamic
and bulk viscosity $\eta$ and $\eta^{\mathrm{b}}$, the thermal conductivity
$k^{\mathrm{th}}$, the specific heat $c_{p}$, the thermal expansion
coefficient ${\alpha_{p}}$, the ratio of specific heats $\gamma=c_{p}/c_{v}$,
and the isentropic and isothermal compressibilities $\kappa_{s}$ and
$\kappa_{T}=\gamma\kappa_{s}$. The governing equations are the transport
equations for the density of mass $\rho$, momentum $\rho\bm{v}$, and internal
energy $\rho\epsilon$, (Landau1993, ; Muller2014, ; Karlsen2015, )

$\displaystyle\partial_{t}\rho$ $\displaystyle=-\bm{\nabla}\cdot(\rho\bm{v}),$
(3a) $\displaystyle\partial_{t}(\rho\bm{v})$
$\displaystyle=\bm{\nabla}\cdot(\bm{\sigma}-\rho\bm{v}\bm{v}),$ (3b)
$\displaystyle\partial_{t}\Big{(}\rho\epsilon+\rho\frac{v^{2}}{2}\Big{)}$
$\displaystyle=\bm{\nabla}\cdot\Big{[}k^{\mathrm{th}}\bm{\nabla}T+\bm{v}\cdot\bm{\sigma}-\rho\bm{v}\Big{(}\epsilon+\frac{v^{2}}{2}\Big{)}\Big{]}+P,$
(3c)

where $P$ is the external heat power density, and $\bm{\sigma}$ is the stress
tensor,

$\displaystyle\bm{\sigma}$ $\displaystyle=-p\bm{I}+\bm{\tau},$ (4a)
$\displaystyle\bm{\tau}$
$\displaystyle=\eta\Big{[}\bm{\nabla}\bm{v}+(\bm{\nabla}\bm{v})^{\dagger}\Big{]}+\Big{(}\eta^{\mathrm{b}}-\frac{2}{3}\eta\Big{)}(\bm{\nabla}\cdot\bm{v})\bm{I}.$
(4b)

Pressure and temperature are related to the internal energy density by the
first law of thermodynamics, and to the density by the equation of state,
(Landau1980, ; Muller2014, ; Karlsen2015, )

$\displaystyle\rho\mathrm{d}\epsilon$ $\displaystyle=(\rho
c_{p}-{\alpha_{p}}p)\>\mathrm{d}T-(\kappa_{T}p+{\alpha_{p}}T)\>\mathrm{d}p$
(5a) $\displaystyle\mathrm{d}\rho$
$\displaystyle=\rho\kappa_{T}\>\mathrm{d}p-\rho{\alpha_{p}}\>dT$ (5b)

The thermodynamics also shows up in the temperature and density
dependency(Muller2014, ) of any material parameter $q$,

$\mathrm{d}q=\Big{(}\frac{\partial q}{\partial
T}\Big{)}_{\rho}\mathrm{d}T+\Big{(}\frac{\partial
q}{\partial\rho}\Big{)}_{T}\mathrm{d}\rho.$ (6)

The temperature sensitivity of each parameter is quantified by the
dimensionless quantity $a_{q}=\frac{1}{{\alpha_{p}}q}\big{(}\frac{\partial
q}{\partial T}\big{)}_{\rho}$,

$\begin{array}[]{rcrrcrrcr}a_{\rho}&=&-1,&\quad a_{\eta}&=&-89,&\quad
a_{\eta^{\mathrm{b}}}&=&-100,\\\ a_{k^{\mathrm{th}}}&=&11,&\quad
a_{{\alpha_{p}}}&=&145,&\quad a_{\kappa_{s}}&=&-10,\end{array}$ (7)

where the values are for water at $T=25\,^{\circ}\\!\textrm{C}$. (Muller2014,
) The temperature dependency of the parameters implies that thermal gradients
may induce gradients in, say, density and compressibility. This leads to the
appearance of the inhomogeneous acoustic body force $\bm{f}_{\mathrm{ac}}$
introduced in acoustofluidics for solute-induced gradients by Karlsen et al.
(Karlsen2016, )

###  II.2 Acoustic actuation and perturbation expansion

Following Ref. Bach2018, , we actuate time-harmonically with angular frequency
$\omega$ by a displacement of a surface, so an element at equilibrium position
$\bm{s}_{0}$, at time $t$ will have the position
$\bm{s}(\bm{s}_{0},t)=\bm{s}_{0}+\bm{s}_{1}(\bm{s}_{0})\mathrm{e}^{-{\mathrm{i}\omega
t}}$. For models containing only a fluid, the displacement will be on the
fluid boundary, whereas for models containing both a fluid and a solid domain,
the actuation is on the solid boundary. For models including the piezoelectric
transducer driving the system, the actuation parameter is the applied
voltage.Skov2019  However, this is not included in this work.

The acoustic response to the actuation parameter $\bm{s}_{1}$ is linear, and
the resulting fields will be complex fields
$Q_{1}(\bm{r})\mathrm{e}^{-{\mathrm{i}\omega t}}$, the so-called first-order
fields with subscript $"1"$. The non-linearity of the governing equation
results in higher order responses to the actuation. We are only interested in
the time-averaged second-order response and define
$Q_{2}(\bm{r})=\big{\langle}Q_{2}(\bm{r},t)\big{\rangle}=\frac{\omega}{2\pi}\int_{0}^{\frac{2\pi}{\omega}}Q_{2}(\bm{r},t)\,\mathrm{d}t$.
A time-average of a product of two first-order fields is also a second-order
term, written as
$\big{\langle}A_{1}B_{1}\big{\rangle}=\frac{1}{2}\operatorname{Re}\big{[}A_{1}B_{1}^{*}\big{]}$,
where the asterisk denote complex conjugation. Thus, a given field
$Q(\bm{r},t)$ in the model, such as density $\rho$, temperature $T$, pressure
$p$, velocity $\bm{v}$, displacement $\bm{u}$, and stress $\bm{\sigma}$, is
written as the sum of the unperturbed field, the acoustic response, and the
time-averaged response,

$Q(\bm{r},t)=Q_{0}(\bm{r})+Q_{1}(\bm{r})\,\mathrm{e}^{-{\mathrm{i}\omega
t}}+Q_{2}(\bm{r}).$ (8)

Similarly, through their dependency on temperature and density, all material
parameters, such as thermal conductivity $k^{\mathrm{th}}$, compressibility
$\kappa$, and (for liquids) viscosity $\eta$, are written as exemplified by
the viscosity,

$\displaystyle\eta(\bm{r},t)$
$\displaystyle=\eta_{0}(T_{0})+\eta_{1}(T_{1},\rho_{1})\,\mathrm{e}^{-{\mathrm{i}\omega
t}},$ (9a) $\displaystyle\eta_{1}(T_{1},\rho_{1})$
$\displaystyle=\Big{(}\frac{\partial\eta}{\partial
T}\Big{)}_{T_{0}}T_{1}(\bm{r})+\Big{(}\frac{\partial\eta}{\partial\rho}\Big{)}_{T_{0}}\rho_{1}(\bm{r}).$
(9b)

###  II.3 Separation of length scales

Acoustofluidic systems exhibits dynamics on two lengthscales, set by the
acoustic wavelength and the thermoviscous boundary layer width. The boundary
conditions on the temperature, heat flux, velocity, and stress at a fluid-
solid interface result in the appearance of a thermal boundary layer (in
fluids and solids) of width $\delta_{\mathrm{t}}$ and in a viscous boundary
layer (in fluids only) of width $\delta_{\mathrm{s}}$, localized near fluid-
solid interfaces. Their dynamically-defined widths, jointly referred to as
$\delta$, are small compared to a typical device size or wavelength $d$, so
$\delta\ll d$,(Karlsen2015, )

$\delta_{\mathrm{s}}=\sqrt{\frac{2\nu_{0}}{\omega}},\qquad\delta_{\mathrm{t}}=\sqrt{\frac{2D^{\mathrm{th}}_{0}}{(1-X)\omega}}\approx\sqrt{\frac{2D^{\mathrm{th}}_{0}}{\omega}},$
(10)

where $X=0$ for fluids and
$X=(\gamma-1)\frac{4c^{2}_{\mathrm{tr}}}{3c^{2}_{\mathrm{lo}}}\lesssim 0.01$
for solids, $\nu_{0}=\frac{\eta_{0}}{\rho_{0}}$, and
$D^{\mathrm{th}}_{0}=\frac{k^{\mathrm{th}}_{0}}{\rho_{0}c_{p0}}$. Typically,
$\delta_{\mathrm{t}}\lesssim\delta_{\mathrm{s}}\lesssim 500~{}\textrm{nm}$,
which is more than two orders of magnitude smaller than $d\sim
100~{}\upmu\textrm{m}$. In this paper, the various fields are decomposed into
a bulk field ($d$) and a boundary-layer field ($\delta$) that are connected by
the boundary conditions. In Fig. 1, this decomposition is sketched near the
fluid-solid boundary for the acoustic temperature field $T_{1}$. Also shown
are the boundary-layer widths $\delta_{\mathrm{s}}$ and $\delta_{\mathrm{t}}$
together with the instantaneous position
$\bm{s}(t)=\bm{s}_{0}+\bm{s}_{1}(\bm{s}_{0},t)$ of the oscillating boundary.

###  II.4 Boundary conditions

In the usual Lagrangian picture,(Bach2018, ) an element with equilibrium
position $\bm{s}_{0}$ in an elastic solid has at time $t$ the position
$\bm{s}(\bm{s}_{0},t)=\bm{s}_{0}+\bm{s}_{1}(\bm{s}_{0})\mathrm{e}^{-{\mathrm{i}\omega
t}}$ and velocity
$\bm{V}^{0}=-\partial_{t}\bm{s}=\bm{V}_{1}^{0}(\bm{s}_{0})\>\mathrm{e}^{-\mathrm{i}\omega
t}$ with $\bm{V}_{1}^{0}(\bm{s}_{0})=-\mathrm{i}\omega\bm{s}_{1}(\bm{s}_{0})$.
On the solid-fluid interface, the no-slip condition applies, so the velocity
of the solid wall at a given time and position must equal the Eulerian-picture
fluid velocity $\bm{v}^{\mathrm{fl}}$,

$\bm{v}^{\mathrm{fl}}(\bm{s}_{0}+\bm{s}_{1}\mathrm{e}^{-\mathrm{i}\omega
t},t)=\bm{V}_{1}^{0}(\bm{s}_{0})\>\mathrm{e}^{-{\mathrm{i}\omega t}}.$ (11)

This boundary condition must be obeyed separately for the first- and second-
order fields (subscript 1 and 2, respectively), so a Taylor expansion
yields(Bach2018, )

$\displaystyle\bm{v}_{1}(\bm{s}_{0})$
$\displaystyle=\bm{V}_{1}^{0}(\bm{s}_{0}),$ (12a)
$\displaystyle\bm{v}_{2}(\bm{s}_{0})$
$\displaystyle=-\big{\langle}(\bm{s}_{1}\cdot\bm{\nabla})\bm{v}_{1}\big{\rangle}\big{|}_{\bm{s}_{0}}=-\frac{1}{\omega}\big{\langle}(\mathrm{i}\bm{V}_{1}^{0}\cdot\bm{\nabla})\bm{v}_{1}\big{\rangle}\big{|}_{\bm{s}_{0}}.$
(12b)

At position $\bm{s}_{0}$ on the fluid-solid interface with surface normal
$\bm{n}$, also the stress $\bm{\sigma}=\bm{\sigma}_{1}+\bm{\sigma}_{2}$ must
be continuous in the first- and second-order contributions $\bm{\sigma}_{1}$
and $\bm{\sigma}_{2}$ separately,

$\displaystyle\bm{\sigma}^{\mathrm{sl}}_{1}(\bm{s}_{0})\cdot\bm{n}$
$\displaystyle=\bm{\sigma}^{\mathrm{fl}}_{1}(\bm{s}_{0})\cdot\bm{n},$ (13a)
$\displaystyle\bm{\sigma}^{\mathrm{sl}}_{2}(\bm{s}_{0})\\!\cdot\\!\bm{n}$
$\displaystyle=\bm{\sigma}^{\mathrm{fl}}_{2}(\bm{s}_{0})\\!\cdot\\!\bm{n}+\big{\langle}(\bm{s}_{1}\\!\cdot\\!\bm{\nabla})\bm{\sigma}^{\mathrm{fl}}_{1}(\bm{s}_{0})\\!\cdot\\!\bm{n}\big{\rangle}\big{|}_{\bm{s}_{0}}.$
(13b)

Here, the thermal effects enter through the temperature dependency of the
viscosity parameters $\eta$ and $\eta^{\mathrm{b}}$, see Eqs. (4) and (9).

Figure 1:  Sketch of the fields at the fluid-solid interface. $\bm{s}_{0}$ is
the equilibrium position of the interface $\partial\Omega$, $\bm{s}_{1}$ the
time-dependent displacement away from $\partial\Omega$, and
$\bm{s}=\bm{s}_{0}+\bm{s}_{1}$ the instantaneous position. The dashed lines
represent the viscous and thermal boundary-layer widths width
$\delta_{\mathrm{s}}$ (black) and $\delta_{\mathrm{t}}$ (red) in the solid and
fluid. $\delta$ without a subscript refers to either $\delta_{\mathrm{s}}$ or
$\delta_{\mathrm{t}}$, and $d$ refers to the bulk lengthscale, so
$\delta_{\mathrm{t}}\lesssim\delta_{\mathrm{s}}\sim\delta\ll d$. The
temperature $T_{1}^{\mathrm{xl}}$ (black) is the sum of a bulk field
$T_{1}^{d,\mathrm{xl}}$ (blue) and a boundary-layer field
$T_{1}^{\delta\mathrm{xl}}$ (red).

Similarly, the temperature $T=T_{0}+T_{1}+T_{2}$ must be continuous across the
solid-fluid interface in each order separately,

$\displaystyle T_{i}^{\mathrm{sl}}(\bm{s}_{0})$
$\displaystyle=T_{i}^{\mathrm{fl}}(\bm{s}_{0}),\quad i=0,1,$ (14a)
$\displaystyle T_{2}^{\mathrm{sl}}(\bm{s}_{0})$
$\displaystyle=T_{2}^{\mathrm{fl}}(\bm{s}_{0})+\big{\langle}\bm{s}_{1}\cdot\bm{\nabla}T^{\mathrm{fl}}_{1}\big{\rangle}\big{|}_{\bm{s}_{0}}.$
(14b)

Also the heat flux $\bm{n}\cdot(-k^{\mathrm{th}}\bm{\nabla}T)$ must be
continuous across the interface,

$k^{\mathrm{th,sl}}\bm{n}\cdot\bm{\nabla}T^{\mathrm{sl}}(\bm{s}_{0},t)=k^{\mathrm{th,fl}}\bm{n}\cdot\bm{\nabla}T^{\mathrm{fl}}(\bm{s}_{0}+\bm{s}_{1}\mathrm{e}^{-\mathrm{i}\omega
t},t),$ (15)

which order by order becomes,

$\displaystyle
k^{\mathrm{th,sl}}_{0}\bm{n}\cdot\bm{\nabla}T_{i}^{\mathrm{sl}}(\bm{s}_{0})=k^{\mathrm{th,fl}}_{0}\bm{n}\cdot\bm{\nabla}T_{i}^{\mathrm{fl}}(\bm{s}_{0}),\quad
i=0,1,$ (16a) $\displaystyle
k^{\mathrm{th,sl}}_{0}\bm{n}\cdot\bm{\nabla}T_{2}^{\mathrm{sl}}+k^{\mathrm{th,sl}}_{2}\bm{n}\cdot\bm{\nabla}T_{0}^{\mathrm{sl}}+\big{\langle}k^{\mathrm{th,sl}}_{1}\bm{n}\cdot\bm{\nabla}T_{1}^{\mathrm{sl}}\big{\rangle}$
(16b)
$\displaystyle=k^{\mathrm{th,fl}}_{0}\bm{n}\cdot\bm{\nabla}T_{2}^{\mathrm{fl}}(\bm{s}_{0})+\big{\langle}k^{\mathrm{th,fl}}_{1}\bm{n}\cdot\bm{\nabla}T_{1}^{\mathrm{fl}}(\bm{s}_{0})\big{\rangle}$
$\displaystyle\quad+k^{\mathrm{th,fl}}_{2}\bm{n}\cdot\bm{\nabla}T_{0}^{\mathrm{fl}}(\bm{s}_{0})+\big{\langle}\bm{s}_{1}\cdot\bm{\nabla}\big{[}k^{\mathrm{th,fl}}_{0}\bm{\nabla}T_{1}^{\mathrm{fl}}(\bm{s}_{0})\big{]}\cdot\bm{n}\big{\rangle}$
$\displaystyle\quad+\big{\langle}\bm{s}_{1}\cdot\bm{\nabla}\big{[}k^{\mathrm{th,fl}}_{1}\bm{\nabla}T_{0}^{\mathrm{fl}}(\bm{s}_{0})\big{]}\cdot\bm{n}\big{\rangle}.$

###  II.5 Range of validity of the model

We briefly discuss the range of validity imposed by the main assumptions.
Firstly, perturbation theory is valid when lower-order terms are much larger
than and unaffected by higher-order terms, say,
$\rho_{0}\gg\big{|}{\rho_{1}}\big{|}$ and
$\big{|}{\bm{v}_{1}}\big{|}\gg\big{|}{\bm{v}_{2}}\big{|}$, and when the latter
can be neglected in the governing equations. For example, the zeroth-order
heat equation (18) is only valid, if the time scale for advective heat
transport $t_{\mathrm{adv}}=d_{t}/\big{|}{\bm{v}_{2}}\big{|}$ is much longer
than that of diffusion $t_{\mathrm{dif}}=d_{t}^{2}/D^{\mathrm{th}}_{0}$ in a
system with characteristic length $d_{t}$. For $d_{t}=1\,\textrm{mm}$ this
requires
$\big{|}{\bm{v}_{2}}\big{|}\ll\frac{D^{\mathrm{th}}_{0}}{d_{t}}\approx
150\,\textrm{\textmu{}m}/\textrm{s}$.

Secondly, due to low oscillatory advection, we assume
$\bm{\nabla}\cdot(q_{0}\bm{v}_{1})\approx q_{0}\bm{\nabla}\cdot\bm{v}_{1}$,
where $q_{0}$ is a parameter of the fluid. This requires
$|q_{0}\bm{\nabla}\cdot\bm{v}_{1}|\gg|\bm{\nabla}q_{0}\cdot\bm{v}_{1}|$. By
using the parameter $a_{q}$ of Eq. (7), the validity of our theory is limited
by,

$\big{|}{\bm{\nabla}T_{0}}\big{|}\ll\bigg{|}\frac{k_{\mathrm{c}}}{a_{\eta}{\alpha_{p0}}}\bigg{|}\approx
5000~{}\frac{\textrm{K}}{\textrm{mm}}.$ (17)

Here, $a_{\eta}$ is used as the viscosity has the strongest temperature
dependency. In conventional acoustofluidic systems
$\big{|}{\bm{\nabla}T_{0}}\big{|}\lesssim 50~{}\textrm{K}/\textrm{mm}\ll
5000~{}\textrm{K}/\textrm{mm}$.

Thirdly, the effective boundary-layer theory requires the boundary-layer width
to be much smaller than the bulk wavelength, $k_{0}\delta\ll 1$, see Section
II.3, which is true for MHz acoustics in water.

##  III Zeroth order: steady background fields

Before turning on the acoustics, $p_{0}$ is constant, and $\bm{v}_{0}=\bm{0}$
in the acoustofluidic system. The temperature $T_{0}$ is determined by
boundary conditions set by the surroundings and the heat power density $P_{0}$
from given sources and sinks. $T_{0}$ is governed by the energy conservation
(3c) to zeroth order in the acoustic actuation,

$0=\bm{\nabla}\cdot\Big{[}k^{\mathrm{th}}_{0}\bm{\nabla}T_{0}\Big{]}+P_{0}.$
(18)

$T_{0}$ determines the zeroth-order water parameters, such as
$\rho_{0}(T_{0})$ and $\eta_{0}(T_{0})$, and thereby affects the resonance
frequency and the Q-factor of the acoustofluidic system.

##  IV First order: acoustics

For the first-order fields, we solve the viscous and thermal boundary layers
analytically, and use these solutions to derive a set of effective boundary
conditions for the bulk fields. The analysis is based on our previous work:
the governing equations derived in Refs. Muller2014, ; Karlsen2015, , the
potential theory derived in Ref. Karlsen2015, , and the effective boundary
method derived for viscous, but not for thermal, boundary layers in Ref.
Bach2018, . The result is a model, where we solve for the displacement field
$\bm{u}_{1}$ in the solid, and for the pressure $p_{1}$ in the fluid, and both
these bulk fields are subject to the effective boundary conditions that
implicitly contain the boundary layers. The temperature $T_{1}$ is
incorporated through $p_{1}$, $\bm{u}_{1}$ in the first-order equations and in
the effective boundary conditions.

###  IV.1 Acoustic equations and potential theory for fluids

The governing equations for the complex-valued acoustic field amplitudes in a
fluid are given in Eq. (11) of Ref. Muller2014, : the mass continuity
equation, the momentum equation, and the heat equation, which couple together
the pressure $p_{1}$, the velocity $\bm{v}_{1}$, and the temperature $T_{1}$,

$\displaystyle-\mathrm{i}\omega{\alpha_{p0}}\>T_{1}+\mathrm{i}\omega\kappa_{T0}\>p_{1}=\bm{\nabla}\cdot\bm{v}_{1},$
(19a)
$\displaystyle-\mathrm{i}\omega\rho_{0}\>\bm{v}_{1}=-\bm{\nabla}p_{1}+\beta\eta_{0}\bm{\nabla}(\bm{\nabla}\cdot\bm{v}_{1})+\eta_{0}\nabla^{2}\bm{v}_{1},$
(19b) $\displaystyle-\mathrm{i}\omega
T_{1}+\mathrm{i}\omega(\gamma-1)\frac{\kappa_{s0}}{{\alpha_{p0}}}\>p_{1}=D^{\mathrm{th}}_{0}\nabla^{2}T_{1},$
(19c)

where $\beta=\frac{\eta^{\mathrm{b}}_{0}}{\eta_{0}}-\frac{2}{3}$. Following
Ref. Karlsen2015, , these equations are solved using potential theory based on
the standard Helmholtz decomposition of the velocity field,
$\bm{v}_{1}=\bm{\nabla}(\phi_{c}+\phi_{t})+\bm{\nabla}\times\bm{\Psi}=\bm{v}^{d}_{1}+\bm{v}^{\delta}_{1}$,
where $\phi_{c}$ is the compressional potential, $\phi_{t}$ is the thermal
potential, and $\bm{\Psi}$ is the shear vector potential. At the fluid-solid
interface $\big{|}T^{\delta}_{1}\big{|}\approx\big{|}T^{d}_{1}\big{|}$, and
combining this with
$T_{1}=T^{d}_{1}+T^{\delta}_{1}=\frac{\mathrm{i}(\gamma-1)\omega}{{\alpha_{p0}}c_{0}^{2}}\phi_{c}+\frac{1}{{\alpha_{p0}}D^{\mathrm{th}}_{0}}\phi_{t}$
with the typical acoustofluidic parameter values inserted, we can deduce
$\big{|}\phi_{t}\big{|}\approx(\gamma-1)\frac{\omega
D^{\mathrm{th}}_{0}}{c_{0}^{2}}\big{|}\phi_{c}\big{|}\approx
10^{-8}\big{|}\phi_{c}\big{|}\ll\big{|}\phi_{c}\big{|}$. From this follows
that
$p_{1}\approx\mathrm{i}\omega\rho_{0}(1+\mathrm{i}\Gamma_{\mathrm{s}})\phi_{c}$,
and we replace $\phi_{c}$, $\phi_{t}$, and $\bm{\Psi}$ by $p_{1}$,
$T^{\delta}_{1}$, and $\bm{v}^{\delta}_{1}$,

$p_{1}\approx\mathrm{i}\omega\rho_{0}(1+\mathrm{i}\Gamma_{\mathrm{s}})\phi_{c},\;\;T^{\delta}_{1}=\frac{\phi_{t}}{{\alpha_{p0}}D^{\mathrm{th}}_{0}},\;\;\bm{v}^{\delta}_{1}=\bm{\nabla}\times\bm{\Psi}.$
(20)

Finally, using the smallness of the damping coefficients,
$\Gamma_{\mathrm{s}}=\frac{1}{2}(1+\beta)(k_{0}\delta_{\mathrm{s}})^{2}\ll 1$
and $\Gamma_{\mathrm{t}}=\frac{1}{2}(k_{0}\delta_{\mathrm{t}})^{2}\ll 1$, with
$k_{0}=\frac{\omega}{c}$, approximate solutions to Eq. (19) are obtained from
the potentials solving three Helmholtz equations,

$\displaystyle\nabla^{2}p_{1}$ $\displaystyle=-k^{2}_{\mathrm{c}}p_{1},$
$\displaystyle\quad k_{\mathrm{c}}$
$\displaystyle=\frac{\omega}{c_{0}}(1+\mathrm{i}\Gamma^{\mathrm{fl}}_{\mathrm{0c}}),$
(21a) $\displaystyle\nabla^{2}T^{\delta}_{1}$
$\displaystyle=-k^{2}_{\mathrm{t}}T^{\delta}_{1},$ $\displaystyle
k_{\mathrm{t}}$
$\displaystyle=\frac{1+\mathrm{i}}{\delta_{\mathrm{t}}}(1+\mathrm{i}\Gamma^{\mathrm{fl}}_{\mathrm{0t}}),$
(21b) $\displaystyle\nabla^{2}\bm{v}^{\delta}_{1}$
$\displaystyle=-k^{2}_{\mathrm{s}}\bm{v}^{\delta}_{1},$ $\displaystyle
k_{\mathrm{s}}$ $\displaystyle=\frac{1+\mathrm{i}}{\delta_{\mathrm{s}}}.$
(21c)

Here,
$\Gamma^{\mathrm{fl}}_{\mathrm{0c}}=\frac{1}{2}\big{[}\Gamma_{\mathrm{s}}+(\gamma-1)\Gamma_{\mathrm{t}}\big{]}$
and
$\Gamma^{\mathrm{fl}}_{\mathrm{0t}}=\frac{\gamma-1}{2}\big{[}\Gamma_{\mathrm{s}}-\Gamma_{\mathrm{t}}\big{]}$
are the resulting damping coefficients, whereas the complex-valued wave
numbers $k_{\mathrm{s}}$ and $k_{\mathrm{t}}$ reveals the existence of the
viscous and thermal boundary layers of thickness $\delta_{\mathrm{s}}$ and
$\delta_{\mathrm{t}}$, respectively, see Fig. 1. The full velocity
$\bm{v}_{1}$ and temperature $T_{1}$ are given by $p_{1}$,
$\bm{v}^{\delta}_{1}$, and $T^{\delta}_{1}$ as,

$\displaystyle\bm{v}_{1}$
$\displaystyle=\bm{v}^{d}_{1}+\bm{v}^{\delta}_{1}\qquad=$
$\displaystyle\bm{v}_{1}^{d,p}$
$\displaystyle\;+\bm{v}_{1}^{d,T}+\bm{v}^{\delta}_{1},$ (22a)
$\displaystyle\bm{v}_{1}^{d,p}$
$\displaystyle=\bm{\nabla}\big{[}-\mathrm{i}\mbox{$\frac{1-\mathrm{i}\Gamma_{\mathrm{s}}}{\omega\rho_{0}}$}\>p_{1}\big{]},\;\;$
$\displaystyle\bm{v}_{1}^{d,T}$
$\displaystyle=\bm{\nabla}\big{[}{\alpha_{p0}}D^{\mathrm{th}}_{0}T_{1}^{\delta}\big{]},$
(22b) $\displaystyle T_{1}$ $\displaystyle=T^{d}_{1}+T^{\delta}_{1},\quad$
$\displaystyle T^{d}_{1}$
$\displaystyle=(\gamma-1)\frac{\kappa_{s0}}{{\alpha_{p0}}}\>p_{1}.$ (22c)

Note that both $\bm{v}_{1}^{d,p}$ and $\bm{v}_{1}^{d,T}$ are gradient fields
in the Helmholtz decomposition, but that $\bm{v}_{1}^{d,T}$ despite its
superscript “d” is a boundary-layer field. Because $T_{1}$ is split into a
bulk and a boundary layer field, the material parameters $q=q_{0}+q_{1}$ are
split similarly. For example, the first-order viscosity $\eta_{1}$ introduced
in Eq. (9) (and similar for other material parameters) becomes

$\eta_{1}=\partial_{T}\eta_{0}\;(T^{d}_{1}+T^{\delta}_{1})+\partial_{\rho}\eta_{0}(\rho_{1}^{d}+\rho_{1}^{\delta})=\eta_{1}^{d}+\eta_{1}^{\delta}.$
(23)

###  IV.2 Acoustic equations and potential theory for solids

For a linear elastic isotropic solid with density $\rho_{0}$, longitudinal
sound speed $c_{\mathrm{lo}}$, and transverse sound speed $c_{\mathrm{tr}}$,
the governing equations is the linearized form of the momentum and heat
equation (1) for the displacement field $\bm{u}_{1}$ and the temperature
$T_{1}$,(Karlsen2015, )

$\displaystyle-\omega^{2}\rho_{0}\bm{u}_{1}=$
$\displaystyle-\frac{{\alpha_{p0}}}{\kappa_{T0}}\bm{\nabla}T_{1}$
$\displaystyle+(c^{2}_{\mathrm{lo}}-c^{2}_{\mathrm{tr}})\bm{\nabla}(\bm{\nabla}\cdot\bm{u}_{1})+c^{2}_{\mathrm{tr}}\nabla^{2}\bm{u}_{1},$
(24a) $\displaystyle-\mathrm{i}\omega T_{1}-$
$\displaystyle\;\mathrm{i}\omega\frac{\gamma-1}{{\alpha_{p0}}}\>\bm{\nabla}\cdot\bm{u}_{1}=D^{\mathrm{th}}_{0}\nabla^{2}T_{1}.$
(24b)

In analogy with the fluid, the governing equations for the solid are solved by
potential theory, again following Ref. Karlsen2015, . The displacement field
is Helmholtz decomposed as
$-\mathrm{i}\omega\bm{u}_{1}=\bm{\nabla}(\phi_{c}+\phi_{t})+\bm{\nabla}\times\bm{\Psi}=-\mathrm{i}\omega(\bm{u}^{\mathrm{lo}}_{1}+\bm{u}^{\mathrm{tr}}_{1})$,
where $\phi_{c}$ is the compressional potential, $\phi_{t}$ is the thermal
potential, and $\bm{\Psi}$ is the shear vector potential, and where we have
used $\bm{v}_{1}^{\mathrm{sl}}=-\mathrm{i}\omega\bm{u}_{1}$. Using the same
approximations as for the fluid, we have
$T_{1}=T^{d}_{1}+T^{\delta}_{1}=\frac{\mathrm{i}(\gamma-1)\omega}{{\alpha_{p0}}c_{0}^{2}}\phi_{c}+\frac{1}{\chi{\alpha_{p0}}D^{\mathrm{th}}_{0}}\phi_{t}$.
We keep $\phi_{c}$, but use
$T^{\delta}_{1}=\frac{1}{\chi{\alpha_{p0}}D^{\mathrm{th}}_{0}}\phi_{t}$
instead of $\phi_{t}$, and
$\bm{u}^{\mathrm{tr}}_{1}=\bm{\nabla}\times\bm{\Psi}$ instead of $\bm{\Psi}$.
The solution to Eq. (24) is obtained from the potentials solving the following
three Helmholtz equations,

$\displaystyle\nabla^{2}\phi_{c}$ $\displaystyle=-k^{2}_{\mathrm{c}}\phi_{c},$
$\displaystyle\quad k_{\mathrm{c}}$
$\displaystyle=\frac{\omega}{c_{0}}(1+\mathrm{i}\Gamma^{\mathrm{sl}}_{\mathrm{0c}}),$
(25a) $\displaystyle\nabla^{2}T^{\delta}_{1}$
$\displaystyle=-k^{2}_{\mathrm{t}}T^{\delta}_{1},$ $\displaystyle
k_{\mathrm{t}}$
$\displaystyle=\frac{1+\mathrm{i}}{\delta_{\mathrm{t}}}(1+\mathrm{i}\Gamma^{\mathrm{sl}}_{\mathrm{0t}}),$
(25b) $\displaystyle\nabla^{2}\bm{u}^{\mathrm{tr}}_{1}$
$\displaystyle=-k^{2}_{\mathrm{s}}\bm{u}^{\mathrm{tr}}_{1},$ $\displaystyle
k_{\mathrm{s}}$ $\displaystyle=\frac{\omega}{c_{\mathrm{tr}}}.$ (25c)

Here, $c_{0}^{2}=c^{2}_{\mathrm{lo}}+\frac{\gamma-1}{\rho_{0}\kappa_{T0}}$,
$\Gamma^{\mathrm{sl}}_{\mathrm{0c}}=\frac{\gamma-1}{2}\chi\Gamma_{\mathrm{t}}$,
and
$\Gamma^{\mathrm{sl}}_{\mathrm{0t}}=\frac{\gamma^{2}\Gamma_{\mathrm{t}}}{8(1-X)}$
are damping coefficients, $\delta_{\mathrm{t}}$ and $\Gamma_{t}$ are given by
Eq. (10), $\chi=1-\frac{4c^{2}_{\mathrm{tr}}}{3c^{2}}\approx\frac{1}{2}$, and
$X=(\gamma-1)\frac{4c^{2}_{\mathrm{tr}}}{3c^{2}}\approx\frac{\gamma-1}{2}$.
For a solid, only $T^{\delta}_{1}$ is a dampened field confined to the
boundary layer, whereas $\phi_{c}$ and $\bm{u}^{\mathrm{tr}}_{1}$ are bulk
fields. The transverse waves in fluids and solids are qualitatively different:
$\bm{v}^{\delta}_{1}$ cannot propagate in a fluid and is restricted to the
boundary layer, whereas $\bm{u}^{\mathrm{tr}}_{1}$ can propagate in a solid
and is not associated with a boundary layer. The full displacement
$\bm{u}_{1}$ and temperature $T_{1}$ are given by $\phi_{c}$,
$\bm{u}^{\mathrm{tr}}_{1}$, and $T^{\delta}_{1}$ as,

$\displaystyle\bm{u}_{1}$
$\displaystyle=\bm{u}^{\mathrm{lo}}_{1}+\bm{u}^{\mathrm{tr}}_{1},\qquad$
$\displaystyle\bm{u}^{\mathrm{lo}}_{1}$
$\displaystyle=\frac{\mathrm{i}}{\omega}\>\bm{\nabla}\phi_{c},$ (26a)
$\displaystyle T_{1}$ $\displaystyle=T^{d}_{1}+T^{\delta}_{1},\qquad$
$\displaystyle T^{d}_{1}$
$\displaystyle=\frac{\mathrm{i}(\gamma-1)\omega}{{\alpha_{p0}}c_{0}^{2}}\phi_{c}.$
(26b)

The explicit expression for the stress tensor $\bm{\sigma}^{\mathrm{xl}}_{1}$
in the fluid (xl = fl) and in the solid (xl = sl) can be formulated jointly in
potential theory as(Karlsen2015, )

$\displaystyle\bm{\sigma}_{1}^{\mathrm{xl}}$
$\displaystyle=-p_{1}^{\mathrm{xl}}\mathbf{I}+\eta_{0}^{\mathrm{xl}}\big{[}(2k^{2}_{\mathrm{c}}-k^{2}_{\mathrm{s}})\phi_{\mathrm{c}}+(2k^{2}_{\mathrm{t}}-k^{2}_{\mathrm{s}})\phi_{\mathrm{t}}\big{]}\mathbf{I}$
$\displaystyle\quad+\eta_{0}^{\mathrm{xl}}\big{[}\bm{\nabla}\bm{v}_{1}^{\mathrm{sl}}+(\bm{\nabla}\bm{v}_{1}^{\mathrm{sl}})^{\dagger}\big{]},$
(27)

where in the solid $p_{1}^{\mathrm{sl}}=0$,
$\eta_{0}^{\mathrm{sl}}=\frac{\mathrm{i}}{\omega}\>\rho_{0}c_{\mathrm{tr}}^{2}$,
$\bm{v}_{1}^{\mathrm{sl}}=-\mathrm{i}\omega\bm{u}_{1}$.

###  IV.3 The thermal boundary layer

The temperature fields $T_{1}^{\delta,\mathrm{xl}}$ in the fluid (xl = fl) and
the solid (xl = sl) are given by Eqs. (21b) and (25b). Following Ref.
Bach2018,  with $x$ and $y$ parallel to the interface and $z$ perpendicular,
an analytical solution can be found using the thin-boundary-layer
approximation $\nabla^{2}\approx\partial_{z}^{2}$ in these equations in
combination with the condition that the field decays away from the boundary,

$\displaystyle T^{\delta,\mathrm{fl}}_{1}(x,y,z)$ $\displaystyle=T^{\delta
0,\mathrm{fl}}_{1}(x,y)\>\mathrm{e}^{\mathrm{i}k_{\mathrm{t}}^{\mathrm{fl}}z},\;$
$\displaystyle\text{ for }z>0,$ (28a) $\displaystyle
T^{\delta,\mathrm{sl}}_{1}(x,y,z)$ $\displaystyle=T^{\delta
0,\mathrm{sl}}_{1}(x,y)\>\mathrm{e}^{-\mathrm{i}k_{\mathrm{t}}^{\mathrm{sl}}z},\;$
$\displaystyle\text{ for }z<0.$ (28b)

The amplitude of the boundary fields $T^{\delta 0,\mathrm{fl}}_{1}(x,y)$ and
$T^{\delta 0,\mathrm{sl}}_{1}(x,y)$ are determined by the boundary conditions
(14a) and (16a) as follows. The normal vector $\bm{n}=-\bm{e}_{z}$ points away
from the fluid, so $\bm{n}\cdot\bm{\nabla}=-\partial_{z}$, and we obtain

$\displaystyle T^{\delta 0,\mathrm{fl}}_{1}$ $\displaystyle=T^{\delta
0,\mathrm{sl}}_{1}-\Delta T_{1}^{d0},$ (29a) $\displaystyle
k^{\mathrm{th,fl}}_{0}\partial_{z}T^{\delta,\mathrm{fl}}_{1}$
$\displaystyle=k^{\mathrm{th,sl}}_{0}\partial_{z}T^{\delta,\mathrm{sl}}_{1},\;\text{
for }z=0,$ (29b)

where $\Delta
T_{1}^{d0}=-\big{(}T^{d0,\mathrm{fl}}_{1}-T^{d0,\mathrm{sl}}_{1}\big{)}$. From
Eq. (29b) follows the relation,

$T^{\delta
0,\mathrm{fl}}_{1}=-\frac{k^{\mathrm{th,sl}}_{0}k_{\mathrm{t}}^{\mathrm{sl}}}{k^{\mathrm{th,fl}}_{0}k_{\mathrm{t}}^{\mathrm{fl}}}\>T^{\delta
0,\mathrm{sl}}_{1}=-\tilde{Z}\>T^{\delta 0,\mathrm{sl}}_{1},$ (30)

where $\tilde{Z}=Z^{\mathrm{sl}}/Z^{\mathrm{fl}}$ is the ratio of
$Z=k^{\mathrm{th}}_{0}k_{\mathrm{t}}=\sqrt{k^{\mathrm{th}}_{0}c_{p0}\rho_{0}}$
of the solid and the fluid, respectively. Combining Eqs. (29a) and (30) leads
to the final expression for the boundary-layer fields,

$\displaystyle T^{\delta,\mathrm{fl}}_{1}(x,y,z)$
$\displaystyle=-\frac{\tilde{Z}}{1+\tilde{Z}}\>\Delta
T_{1}^{d0}(x,y)\>\mathrm{e}^{\mathrm{i}k_{\mathrm{t}}^{\mathrm{fl}}z},$ (31a)
$\displaystyle T^{\delta,\mathrm{sl}}_{1}(x,y,z)$
$\displaystyle=+\frac{1}{1+\tilde{Z}}\>\Delta
T_{1}^{d0}(x,y)\>\mathrm{e}^{-\mathrm{i}k_{\mathrm{t}}^{\mathrm{sl}}z}.$ (31b)

###  IV.4 The viscous boundary layer

The viscous boundary layer exists only in the fluid, since in the solid both
$\bm{u}^{\mathrm{lo}}_{1}$ and $\bm{u}^{\mathrm{tr}}_{1}$ are bulk fields. The
velocity field in the fluid is given in Eq. (22a) as
$\bm{v}_{1}=\bm{v}^{d}_{1}+\bm{v}^{\delta}_{1}$, where $\bm{v}^{d}_{1}$
depends on the bulk field $p_{1}$ and the boundary field $T^{\delta}_{1}$. The
boundary field $\bm{v}^{\delta}_{1}$ is given by the Helmholtz equation (21c),
to which an analytical solution can be found using the thin-boundary-layer
approximation $\nabla^{2}\approx\partial_{z}^{2}$ in combination with the
condition that the field decays away from the boundary,(Bach2018, )

$\bm{v}^{\delta}_{1}=\bm{v}^{\delta
0}_{1}(x,y)\>\mathrm{e}^{\mathrm{i}k_{s}z}.$ (32)

The amplitude $\bm{v}^{\delta 0}_{1}$ of the boundary field is determined by
the no-slip condition (12a),

$\bm{v}^{\delta
0}_{1}=\bm{V}_{1}^{0}-\bm{v}^{d0}_{1}=-\mathrm{i}\omega\bm{u}^{0}_{1}-\bm{v}^{d0}_{1}.$
(33)

###  IV.5 The effective boundary condition for the velocity

Given the analytical solutions of the three boundary-layer fields, we only
need to numerically solve the three bulk fields, namely $\phi_{c}$ and
$\bm{\Psi}$ in the solid and $\phi_{c}$ in the fluid, or equivalently, the
displacement $\bm{u}_{1}$ in the solid and the pressure $p_{1}$ in the fluid.
Therefore, we set two effective boundary conditions on these bulk fields using
the analytical solutions for the boundary-layer fields: One effective boundary
condition on the displacement $\bm{u}_{1}$ in the solid derived from the
condition on the stress, and another on the pressure in the fluid.

First, from the no-slip condition (12a), we derive the boundary condition for
the first-order pressure field $p_{1}$, which takes the viscous and thermal
boundary-layer effects into account through terms with $k_{\mathrm{s}}$,
$k_{\mathrm{t}}$, and $T_{1}^{\delta 0}$. We express the compressional
velocity $v^{d0,\mathrm{fl}}_{1,z}$ on the fluid-solid interface through the
no-slip condition (33), then use the incompressibility condition on the
boundary-layer velocity, $\mathrm{i}k_{\mathrm{s}}v^{\delta
0,\mathrm{fl}}_{1,z}+\bm{\nabla}\cdot\bm{v}^{\delta 0,\mathrm{fl}}_{1}=0$, to
get rid of the $z$-component $v^{\delta 0,\mathrm{fl}}_{1,z}$, and finally
introduce the bulk fields in the fluid,

$\displaystyle v^{d0,\mathrm{fl}}_{1,z}$
$\displaystyle=v^{d0,\mathrm{sl}}_{1,z}-v^{\delta
0,\mathrm{fl}}_{1,z}=v^{d0,\mathrm{sl}}_{1,z}-\frac{\mathrm{i}}{k_{\mathrm{s}}}\bm{\nabla}\cdot\bm{v}^{\delta
0,\mathrm{fl}}_{1}$
$\displaystyle=v^{d0,\mathrm{sl}}_{1,z}-\frac{\mathrm{i}}{k_{\mathrm{s}}}\bm{\nabla}\cdot\Big{[}\bm{v}^{d0,\mathrm{sl}}_{1}-\bm{v}^{d0,\mathrm{fl}}_{1}\Big{]}$
(34)
$\displaystyle=\Big{[}v^{d0,\mathrm{sl}}_{1,z}\\!-\\!\frac{\mathrm{i}}{k_{\mathrm{s}}}\bm{\nabla}\\!\cdot\\!\bm{v}^{d0,\mathrm{sl}}_{1}\Big{]}+\frac{\mathrm{i}}{k_{\mathrm{s}}}\Big{[}\bm{\nabla}\\!\cdot\\!\bm{v}^{d,\mathrm{fl}}_{1}\\!-\\!\partial_{z}v^{d,\mathrm{fl}}_{1,z}\Big{]}_{z=0}.$

Combining Eqs. (19a) and (22c), we obtain

$\bm{\nabla}\cdot\bm{v}_{1}^{\mathrm{d}}=\mathrm{i}\frac{(1-\mathrm{i}\Gamma_{\mathrm{s}})k_{\mathrm{c}}^{2}}{\omega\rho_{0}}p_{1}-\mathrm{i}\omega{\alpha_{p0}}T_{1}^{\delta}.$
(35a) Then using Eq. (22a), we write $v^{d0,\mathrm{fl}}_{1,z}$ and
$\partial_{z}v^{d,\mathrm{fl}}_{1,z}$ evaluated at the solid-fluid interface
at $z=0$, and arrive at $\displaystyle v^{d0,\mathrm{fl}}_{1,z}$
$\displaystyle=-\frac{\mathrm{i}}{\omega\rho_{0}}(1-\mathrm{i}\Gamma_{\mathrm{s}})\partial_{z}p_{1}+{\alpha_{p0}}D^{\mathrm{th}}_{0}\partial_{z}T_{1}^{\delta},$
(35b) $\displaystyle\partial_{z}v^{d,\mathrm{fl}}_{1,z}$
$\displaystyle=-\frac{\mathrm{i}}{\omega\rho_{0}}(1-\mathrm{i}\Gamma_{\mathrm{s}})\partial_{z}^{2}p_{1}+{\alpha_{p0}}D^{\mathrm{th}}_{0}\partial_{z}^{2}T_{1}^{\delta}.$
(35c)

Inserting Eqs. (35) and (12a) into Eq. (34) leads to the final form of the
effective boundary condition on $p_{1}$,

$\displaystyle\partial_{z}p_{1}$
$\displaystyle=\mathrm{i}\frac{\omega\rho_{0}}{1-\mathrm{i}\Gamma_{\mathrm{s}}}\big{(}V_{1z}^{0}-\frac{\mathrm{i}}{k_{\mathrm{s}}}\bm{\nabla}\cdot\bm{V}_{1}^{0}\big{)}-\frac{\mathrm{i}}{k_{\mathrm{s}}}\big{(}k_{\mathrm{c}}^{2}+\partial_{z}^{2}\big{)}p_{1}$
$\displaystyle\quad+\frac{\mathrm{i}}{k_{\mathrm{t}}}\frac{{\alpha_{p0}}}{\kappa_{T0}}k^{2}_{0}\,T_{1}^{\delta
0},\;\text{ for }z=0.$ (36a) The first two terms on the right-hand were
derived by Bach and Bruus, (Bach2018, ) whereas the last term is a new
correction due to the thermal boundary layer. For $T_{1}^{d}\approx
T_{1}^{\delta}$ at $z=0$, this thermal correction is of the order
$\frac{\gamma-1}{k_{\mathrm{t}}}\>k_{\mathrm{c}}^{2}\>p_{1}$. We emphasize,
that although formulated as an effective boundary condition on the pressure
gradient, Eq. (36a) is the no-slip velocity condition.

###  IV.6 The effective boundary condition for the stress

Next, using the explicit expressions for $\bm{\sigma}_{1}^{\mathrm{sl}}$ and
$\bm{\sigma}_{1}^{\mathrm{fl}}$, we turn to the stress boundary condition
(13a), the continuity of the stress $\bm{\sigma}_{1}$ across the fluid-solid
interface,
$\bm{\sigma}_{1}^{\mathrm{sl}}\cdot\bm{e}_{z}=\bm{\sigma}_{1}^{\mathrm{fl}}\cdot\bm{e}_{z}$.
For the fluid, we use that $k_{\mathrm{s}}\gg k_{\mathrm{c}}$,
$\big{|}{\phi_{\mathrm{c}}}\big{|}\gg\big{|}{\phi_{\mathrm{t}}}\big{|}$, and
$\big{|}{\partial_{z}\bm{v}_{1}^{\delta}}\big{|}\gg\big{|}{\nabla\bm{v}_{1}^{d}}\big{|}$
in Eq. (IV.2), and find

$\bm{\sigma}_{1}^{\mathrm{fl}}\cdot\bm{e}_{z}=-p_{1}\bm{e}_{z}+\mathrm{i}k_{\mathrm{s}}\eta_{0}\Big{[}\bm{v}^{d0,\mathrm{sl}}_{1}+\frac{\mathrm{i}}{\omega\rho_{0}}\bm{\nabla}p_{1}\Big{]}_{\bm{s}_{0}}.$
(36b)

For the solid, we neglect in Eq. (IV.2) the derivative
$\partial_{\parallel}\phi_{\mathrm{t}}$ along the surface, as it is a factor
$\Gamma_{\mathrm{t}}$ smaller than $\partial_{\parallel}\phi_{\mathrm{c}}$.
The remaining $\phi_{\mathrm{t}}$-dependent boundary-layer terms cancel out,
leaving only the bulk-term part $\bm{\sigma}_{1}^{d,\mathrm{sl}}$ of
$\bm{\sigma}_{1}^{\mathrm{sl}}$. The resulting effective stress boundary
condition becomes,

$\bm{\sigma}_{1}^{d,\mathrm{sl}}\cdot\bm{e}_{z}=\bm{\sigma}_{1}^{\mathrm{fl}}\cdot\bm{e}_{z}.$
(36c)

As the thermal boundary-layer fields do not enter, this expression is
identical to the effective boundary condition for the stress derived in Ref.
Bach2018, .

##  V Second order: Acoustic streaming

For the second-order fields in the fluid, we follow Eq. (8) and consider only
the time averaged fields, namely the velocity $\bm{v}_{2}$, pressure $p_{2}$
and stress $\bm{\sigma}_{2}$. The temperature field $T_{2}$ does not enter the
second-order continuity or Navier–Stokes equation, so we drop the heat
equation. The first-order temperature field $T_{1}$ enters the equations
through the material parameters of the fluid,

$\displaystyle 0$
$\displaystyle=-\bm{\nabla}\cdot(\rho_{0}\bm{v}_{2})+\dot{\rho}_{\mathrm{ac}}$
(37a) $\displaystyle 0$
$\displaystyle=-\bm{\nabla}p_{2}+\bm{\nabla}\cdot\bm{\tau}_{2}+\hat{\bm{f}}_{\mathrm{ac}}$
(37b) $\displaystyle\bm{\tau}_{2}$
$\displaystyle=\eta_{0}\big{[}\bm{\nabla}\bm{v}_{2}+(\bm{\nabla}\bm{v}_{2})^{\dagger}\big{]}+\big{[}\eta^{\mathrm{b}}_{0}-\tfrac{2}{3}\eta_{0}\big{]}(\bm{\nabla}\cdot\bm{v}_{2})\,\bm{I},$
(37c) $\displaystyle\bm{v}_{2}^{0}$
$\displaystyle=-\frac{1}{\omega}\big{\langle}\big{(}\mathrm{i}\bm{V}_{1}^{0}\cdot\bm{\nabla}\big{)}\bm{v}_{1}\big{\rangle}\Big{|}_{\bm{r}=\bm{s}_{0}}.$
(37d)

Here, the excess density rate $\dot{\rho}_{\mathrm{ac}}$ and the acoustic body
force $\hat{\bm{f}}_{\mathrm{ac}}$ are time-averaged products of fast varying
first-order fields, which, assuming that
$\rho_{0}\bm{v}_{1}\gg\rho_{1}\bm{v}_{0}$ as is true for typical
acoustofluidic devices, are given by

$\displaystyle\dot{\rho}_{\mathrm{ac}}$
$\displaystyle=-\bm{\nabla}\cdot\big{\langle}\rho_{1}\bm{v}_{1}\big{\rangle},$
(38a) $\displaystyle\hat{\bm{f}}_{\mathrm{ac}}$
$\displaystyle=\bm{\nabla}\cdot\big{[}-\rho_{0}\big{\langle}\bm{v}_{1}\bm{v}_{1}\big{\rangle}+\bm{\tau}_{11}\big{]},$
(38b) $\displaystyle\bm{\tau}_{11}$
$\displaystyle=\big{\langle}\eta_{1}\big{[}\bm{\nabla}\bm{v}_{1}\\!+\\!(\bm{\nabla}\bm{v}_{1})^{\dagger}\big{]}+\big{[}\eta_{1}^{\mathrm{b}}\\!-\\!\tfrac{2}{3}\eta_{1}\big{]}(\bm{\nabla}\cdot\bm{v}_{1})\bm{I}\big{\rangle}.$
(38c)

The slowly varying second-order fields are split up in a bulk field
(superscript “$d$”) and a boundary field (superscript “$\delta$”) according to
their response to the boundary and bulk part of the acoustic force
$\hat{\bm{f}}_{\mathrm{ac}}=\hat{\bm{f}}_{\mathrm{ac}}^{d}+\hat{\bm{f}}_{\mathrm{ac}}^{\delta}$,
and they are coupled by the boundary conditions,

$\displaystyle p_{2}$ $\displaystyle=p_{2}^{d}+p_{2}^{\delta},\qquad$
$\displaystyle\bm{v}_{2}$ $\displaystyle=\bm{v}_{2}^{d}+\bm{v}_{2}^{\delta},$
(39a) $\displaystyle\bm{\tau}_{2}$
$\displaystyle=\bm{\tau}_{2}^{d}+\bm{\tau}_{2}^{\delta},\qquad$
$\displaystyle\bm{\tau}_{11}$
$\displaystyle=\bm{\tau}_{11}^{d}+\bm{\tau}_{11}^{\delta}.$ (39b)

Note that in contrast to the first-order fields, this is not a Helmholtz
decomposition: by definition, a second-order boundary-layer field “$\delta$”
contains at least one first-order boundary-layer field. The computation
strategy for second-order streaming is similar to the one for first-order
acoustics: (1) find analytical solution to the boundary layers, (2) formulate
effective boundary conditions, and (3) solve the bulk fields with the
effective boundary conditions. This decomposition enables simulations of the
bulk fields without resolving the boundary-layer fields.

###  V.1 Short-range boundary-layer streaming

The short-range part “$\delta$” of Eq. (37) is given by the short-range part
of the second-order fields as well as all source terms containing at least one
boundary-layer field,

$\displaystyle 0$
$\displaystyle=\bm{\nabla}\cdot\big{(}\rho_{0}\bm{v}_{2}^{\delta}\big{)}+\dot{\rho}_{\mathrm{ac}}^{\delta},$
(40a) $\displaystyle 0$
$\displaystyle=-\bm{\nabla}p_{2}^{\delta}+\bm{\nabla}\cdot\bm{\tau}_{2}^{\delta}+\hat{\bm{f}}_{\mathrm{ac}}^{\delta},$
(40b) $\displaystyle\qquad\mathrm{where}\;\bm{v}^{\delta}_{2}\rightarrow
0\;\mathrm{as}\;z\rightarrow\infty.$ (40c)

At the boundary, the advection term can be neglected compared to the viscous
term, because of the large gradients induced by the small lengthscale
$\delta$. The thermal boundary layer $T_{1}^{\delta}$ and the associated
boundary-layer velocity $\bm{v}_{1}^{d,T}$ introduce a correction
$\bm{v}_{2}^{\delta,T}$ to the purely viscous boundary-layer term
$\bm{v}_{2}^{\delta,p}$ computed in Ref. Bach2018, ,

$\bm{v}_{2}^{\delta}=\bm{v}_{2}^{\delta,p}+\bm{v}_{2}^{\delta,T}.$ (41)

In the parallel component of $\bm{v}_{2}^{\delta}$, the pressure field can be
neglected, because
$\partial_{\parallel}p_{2}^{\delta}\ll\eta_{0}\partial_{z}^{2}\bm{v}_{2\parallel}^{\delta}$.(Bach2018,
) Thus combining Eqs. (38b) and (40b), the parallel component of the short-
range velocity field $\bm{v}_{2}^{\delta,T}$ obeys,

$\displaystyle\nu_{0}\partial_{z}^{2}\bm{v}_{2,\parallel}^{\delta,T}$
$\displaystyle=\Big{[}\bm{\nabla}\cdot\big{\langle}\bm{v}_{1}^{\delta}\bm{v}_{1}^{d,T}+\bm{v}_{1}^{d,T}\bm{v}_{1}^{\delta}+\bm{v}_{1}^{d,p}\bm{v}_{1}^{d,T}$
$\displaystyle+\bm{v}_{1}^{d,T}\bm{v}_{1}^{d,p}+\bm{v}_{1}^{d,T}\bm{v}_{1}^{d,T}\big{\rangle}-\frac{1}{\rho_{0}}\bm{\nabla}\cdot\bm{\tau}_{11}^{\delta}\Big{]}_{\parallel}.$
(42)

Here, $\bm{\tau}_{11}^{\delta}$ depends on $T_{1}$ through $\eta_{1}(T_{1})$,
whereas the velocity $\bm{v}_{1}^{d,T}$, given in (22a), depends on the
thermal boundary layer $T_{1}^{\delta}$. From Section IV and in particular
Eqs. (22), (28), and (32) follow the relations
$\bm{\nabla}\\!\cdot\\!\bm{v}_{1}^{\delta}=0$,
$\big{|}\bm{v}_{1,\parallel}^{\delta}\big{|}\approx\big{|}\bm{v}_{1,\parallel}^{d}\big{|}$,
$\big{|}v_{1,z}^{\delta}\\!\big{|}\approx(k_{c}\delta_{\mathrm{s}})\big{|}v_{1,z}^{d}\big{|}$,
$\big{|}T_{1}^{\delta}\big{|}\approx\big{|}T_{1}^{d}\big{|}$,
$\bm{\nabla}\cdot\bm{v}_{1}^{d,T}\approx(\gamma\\!-\\!1)\bm{\nabla}\cdot\bm{v}_{1}^{d,p}$,
$\big{|}{\bm{v}_{1,z}^{d,T}}\big{|}\approx(\gamma-1)(k_{\mathrm{c}}\delta_{\mathrm{t}})\big{|}{\bm{v}_{1,z}^{d,p}}\big{|}$,
$\big{|}{\bm{v}_{1,\parallel}^{d,T}}\big{|}\approx(\gamma\\!-\\!1)(k_{\mathrm{c}}\delta_{\mathrm{t}})^{2}\big{|}{\bm{v}_{1}^{d,p}}\big{|}$,
and
$\bm{v}_{1}^{d,T}={\alpha_{p0}}D^{\mathrm{th}}_{0}\,\bm{\nabla}T_{1}^{\delta}$.
To lowest order in $k_{\mathrm{c}}\delta\ll 1$ (involving
$\partial_{z}T_{1}^{\delta}$ and $\partial_{z}\bm{v}_{1}^{\delta}$,
respectively), these relations combined with time averaging
$\big{\langle}a_{1}b_{1}\big{\rangle}=\frac{1}{2}\operatorname{Re}\\!\big{[}a_{1}b_{1}^{*}\big{]}$
change Eq. (V.1) to

$\displaystyle\nu_{0}\partial_{z}^{2}\bm{v}_{2,\parallel}^{\delta,T}$
$\displaystyle=\Big{[}\big{\langle}\big{(}\partial_{z}\bm{v}_{1}^{\delta}\big{)}v_{1,z}^{d,T}\big{\rangle}+\big{\langle}\big{(}\bm{v}_{1}^{\delta}+\bm{v}_{1}^{d,p}\big{)}\big{(}\partial_{z}v_{1,z}^{d,T}\big{)}\big{\rangle}$
$\displaystyle\quad-\frac{1}{\rho_{0}}\Big{(}\big{\langle}(\partial_{z}\eta_{1}^{\delta})\partial_{z}\bm{v}_{1}^{\delta}\big{\rangle}+\big{\langle}(\eta_{1}^{\delta}+\eta_{1}^{d})\nabla^{2}\bm{v}_{1}^{\delta}\big{\rangle}\Big{)}\Big{]}_{\parallel}$
$\displaystyle=\frac{1}{2}\operatorname{Re}\\!\bigg{[}\frac{2{\alpha_{p0}}D^{\mathrm{th}}_{0}}{\delta_{\mathrm{t}}^{2}}\bigg{(}\frac{\delta_{\mathrm{t}}+\mathrm{i}\delta_{\mathrm{s}}}{\delta_{\mathrm{s}}}\bm{v}_{1}^{\delta}+\mathrm{i}\bm{v}_{1}^{d,p}\bigg{)}\,T_{1}^{\delta*}$
$\displaystyle\qquad\quad-\frac{2}{\rho_{0}\delta_{\mathrm{s}}^{2}}\bigg{(}\frac{\delta_{\mathrm{s}}+\mathrm{i}\delta_{\mathrm{t}}}{\delta_{\mathrm{t}}}\eta_{1}^{\delta}+\mathrm{i}\eta_{1}^{d}\bigg{)}\,\bm{v}_{1}^{\delta*}\bigg{]}_{\parallel}.$
(43)

The integration of Eq. (V.1) after $z$ twice, is facilitated by using the
analytical forms (28) and (32) for $\bm{v}_{1}^{d,T}$, $T_{1}^{\delta}$, and
$\bm{v}_{1}^{\delta}$, and by noting that in the boundary layer
$\eta_{1}^{d}\approx\eta_{1}^{d0}\\!+\\!z\partial_{z}\eta_{1}^{d}\approx(1\\!+\\!k_{\mathrm{c}}\delta_{\mathrm{s}})\eta_{1}^{d0}\approx\eta_{1}^{d0}$
and similarly $\bm{v}_{1}^{d}\approx\bm{v}_{1}^{d0}$,

$\displaystyle\bm{v}_{1}^{\delta}$ $\displaystyle=\bm{v}_{1}^{\delta
0}(x,y)\,q(z),\,$ with
$\displaystyle\;q(z)=\mathrm{e}^{\mathrm{i}k_{\mathrm{s}}z},$ (44a)
$\displaystyle T_{1}^{\delta}$ $\displaystyle=T_{1}^{\delta 0}(x,y)\,r(z),$
with $\displaystyle\;r(z)=\mathrm{e}^{\mathrm{i}k_{\mathrm{t}}z},$ (44b)
$\displaystyle\eta_{1}^{\delta}$ $\displaystyle=\eta_{1}^{\delta
0}(x,y)\,r(z),\;$ (44c) $\displaystyle\eta_{1}^{d}$
$\displaystyle\approx\eta_{1}^{d0}\quad\text{ and }$
$\displaystyle\bm{v}_{1}^{d,p}\approx$
$\displaystyle\;\bm{v}_{1}^{d0,p},\text{ for }\;z\ll d.$ (44d)

Following the procedure of Ref. Bach2018, , we introduce the integrals
$I_{{a}{b}}^{(n)}$ of the integrand $a(z)\,b(z)^{*}$, where $a(z)$ and $b(z)$
are any of the functions $1$, $q(z)$, and $r(z)$,

$\displaystyle I_{{a}{b}}^{(n)}$
$\displaystyle=\int^{z}\\!\mathrm{d}z_{n}\int^{z_{n}}\\!\mathrm{d}z_{n-1}\ldots\int^{z_{2}}\\!\mathrm{d}z_{1}\>a(z_{1})\>b(z_{1})^{*}\bigg{|}_{z=0},$
$\displaystyle I_{{a}{b}}^{(n)}$ $\displaystyle\propto\delta^{n},\text{ with
}\,\delta=\delta_{\mathrm{s}},\delta_{\mathrm{t}}\text{ and }n=1,2,3,\ldots.$
(45)

With this notation, Eq. (V.1) is easily integrated to give

$\displaystyle\bm{v}_{2,\parallel}^{\delta 0,T}$
$\displaystyle=\frac{{\alpha_{p0}}D^{\mathrm{th}}_{0}}{\nu_{0}\delta_{\mathrm{t}}^{2}}\operatorname{Re}\\!\bigg{[}\frac{\delta_{\mathrm{t}}+\mathrm{i}\delta_{\mathrm{s}}}{\delta_{\mathrm{s}}}\,I_{{q}{r}}^{(2)}\bm{v}_{1}^{\delta
0}T_{1}^{\delta 0*}+\mathrm{i}I_{{1}{r}}^{(2)}\bm{v}_{1}^{d0,p}T_{1}^{\delta
0*}\bigg{]}_{\parallel}$
$\displaystyle\quad-\frac{1}{\eta_{0}\delta_{\mathrm{s}}^{2}}\operatorname{Re}\\!\bigg{[}\frac{\delta_{\mathrm{s}}+\mathrm{i}\delta_{\mathrm{t}}}{\delta_{\mathrm{t}}}\,I_{{r}{q}}^{(2)}\eta_{1}^{\delta
0}\bm{v}_{1}^{\delta
0*}+\mathrm{i}I_{{1}{q}}^{(2)}\eta_{1}^{d0}\bm{v}_{1}^{\delta
0*}\bigg{]}_{\parallel}.$ (46a) where the integrals are given by
$I_{ba}^{(n)}=\big{[}I_{ab}^{(n)}\big{]}^{*}$ and
$I_{1r}^{(2)}=-\frac{\mathrm{i}}{2}\delta_{\mathrm{t}}^{2},\;\;I_{1q}^{(2)}=-\frac{\mathrm{i}}{2}\delta_{\mathrm{s}}^{2},\;\;I_{rq}^{(2)}=\frac{\mathrm{i}\,\delta_{\mathrm{s}}^{2}\delta_{\mathrm{t}}^{2}}{2(\delta_{\mathrm{s}}+\mathrm{i}\delta_{\mathrm{t}})^{2}}.$
(46b)

When inserting $\rho_{1}^{\delta}=-\rho_{0}{\alpha_{p}}T_{1}^{\delta}$ in the
final expression for the thermal correction, $\bm{v}_{2\parallel}^{\delta
0,T}$ becomes

$\displaystyle\bm{v}_{2,\parallel}^{\delta 0,T}$
$\displaystyle=-\frac{1}{2\rho_{0}}\frac{\delta_{\mathrm{t}}^{2}}{\delta_{\mathrm{s}}^{2}}\operatorname{Re}\Big{[}\frac{\delta_{\mathrm{s}}}{\delta_{\mathrm{s}}-\mathrm{i}\delta_{\mathrm{t}}}\bm{v}_{1}^{\delta
0}\rho_{1}^{\delta 0*}+\bm{v}_{1}^{d0,p}\rho_{1}^{\delta
0*}\Big{]}_{\parallel}$
$\displaystyle\quad-\frac{1}{2\eta_{0}}\operatorname{Re}\Big{[}\frac{\delta_{\mathrm{t}}}{\delta_{\mathrm{t}}-\mathrm{i}\delta_{\mathrm{s}}}\,\eta_{1}^{\delta
0}\bm{v}_{1}^{\delta 0*}+\eta_{1}^{d0}\bm{v}_{1}^{\delta
0*}\Big{]}_{\parallel},$ (47)

where two terms are due to the change in density and two to the change in
viscosity. The perpendicular part of the short-ranged streaming velocity
$v_{2z}^{\delta 0,T}$ can be found by integrating the continuity equation
(37a),
$\partial_{z}v_{2z}^{\delta,T}=-\bm{\nabla}_{\parallel}\cdot\bm{v}_{2\parallel}^{\delta,T}-\frac{1}{\rho_{0}}\bm{\nabla}\cdot\big{\langle}\rho_{1}\bm{v}_{1}\big{\rangle}^{\delta,T}$,
once with respect to $z$,

$v_{2z}^{\delta,T}=-\bm{\nabla}_{\parallel}\cdot\int^{z}\bm{v}_{2\parallel}^{\delta,T}\>\mathrm{d}z-\frac{1}{\rho_{0}}\int^{z}\bm{\nabla}\cdot\big{\langle}\rho_{1}\bm{v}_{1}\big{\rangle}^{\delta,T}\>\mathrm{d}z.$
(48)

The term $\int^{z}\bm{v}_{2\parallel}^{\delta}\>\mathrm{d}z$ is given by Eq.
(46) by substituting all $I_{{a}{b}}^{(2)}$ by $I_{{a}{b}}^{(3)}\propto
I_{{a}{b}}^{(2)}\delta$, so
$\big{|}\bm{\nabla}_{\parallel}\cdot\int^{z}\bm{v}_{2\parallel}^{\delta}\>\mathrm{d}z\big{|}\sim(k_{c}\delta)\,\big{|}\bm{v}_{2\parallel}^{\delta
0}\big{|}$, and
$\int^{z}\bm{\nabla}\cdot\big{\langle}\rho_{1}\bm{v}_{1}\big{\rangle}^{\delta}\>\mathrm{d}z\approx\int^{z}\partial_{z}\big{\langle}\rho_{1}^{\delta}v_{1,z}^{d,p}\big{\rangle}\>\mathrm{d}z=\big{\langle}\rho_{1}^{\delta}v_{1,z}^{d,p}\big{\rangle}$.
Including pre-factors, we obtain to leading order in $k_{c}\delta$,

$v_{2,z}^{\delta
0,T}=-\frac{1}{2\rho_{0}}\operatorname{Re}\\!\Big{[}\rho_{1}^{\delta
0*}v_{1z}^{d0,p}\Big{]}.$ (49)

###  V.2 Bulk field and effective boundary condition

With the short-range boundary-layer streaming term $\bm{v}_{2}^{\delta
0}=\bm{v}_{2}^{\delta 0,p}+\bm{v}_{2}^{\delta 0,T}$ in place, it is now
possible to set up the governing equations and boundary conditions for the
second-order bulk acoustic streaming $\bm{v}_{2}^{d}$,

$\displaystyle 0$
$\displaystyle=\bm{\nabla}\cdot\big{(}\rho_{0}\bm{v}_{2}^{d}\big{)}-\dot{\rho}_{\mathrm{ac}}^{d},$
(50a) $\displaystyle 0$
$\displaystyle=-\bm{\nabla}p_{2}^{d}+\bm{\nabla}\cdot\bm{\tau}_{2}^{d}+\hat{\bm{f}}_{\mathrm{ac}}^{d},$
(50b) $\displaystyle\bm{\tau}_{2}^{d}$
$\displaystyle=\eta_{0}\big{[}\bm{\nabla}\bm{v}_{2}^{d}+(\bm{\nabla}\bm{v}_{2}^{d})^{\dagger}\big{]}+\beta\eta_{0}(\bm{\nabla}\cdot\bm{v}_{2})\,\bm{I},$
(50c) $\displaystyle\bm{v}_{2}^{d0}$ $\displaystyle=-\bm{v}_{2}^{\delta
0}-\frac{1}{\omega}\big{\langle}\big{(}\mathrm{i}\bm{V}^{0}_{1}\cdot\bm{\nabla}\big{)}\bm{v}_{1}\big{\rangle}\big{|}_{\bm{r}=\bm{s}_{0}}.$
(50d)

Here, $\dot{\rho}_{\mathrm{ac}}^{d}$ and $\hat{\bm{f}}_{\mathrm{ac}}^{d}$ are
the bulk terms in Eq. (38). In the mass-conservation equation,
$\bm{\nabla}\cdot\bm{v}_{2}^{d}$ becomes,

$\bm{\nabla}\cdot\bm{v}_{2}^{d}=-\frac{\bm{\nabla}\cdot\big{\langle}\rho_{1}^{d}\bm{v}_{1}^{d,p}\big{\rangle}}{\rho_{0}}=\Gamma\frac{k_{0}\big{|}{\bm{v}_{1}^{d,p}}\big{|}^{2}}{2c_{0}}.$
(51)

Each term of $\bm{\nabla}\cdot\bm{v}_{2}^{d}$ scales as
$\frac{k_{0}}{c_{0}}|\bm{v}_{1}^{d,p}|^{2}\gg\frac{1}{2}\Gamma\frac{k_{0}}{c_{0}}|\bm{v}_{1}^{d,p}|^{2}$,
so
$\frac{1}{\rho_{0}}\bm{\nabla}\cdot\big{\langle}\rho_{1}^{d}\bm{v}_{1}^{d,p}\big{\rangle}$
is negligible compared to the individual terms in
$\bm{\nabla}\cdot\bm{v}_{2}^{d}$. We thus conclude that
$\bm{\nabla}\cdot\bm{v}_{2}^{d}=0$, and that the streaming flow is
incompressible. The acoustic body force $\hat{\bm{f}}_{\mathrm{ac}}^{d}$ may
be expressed as follows, where $\bm{\nabla}\rho_{0}$ and
$\bm{\nabla}\kappa_{s0}$ unlike in previous work (Karlsen2016, ; Karlsen2018,
) can be induced by temperature gradients,

$\displaystyle\hat{\bm{f}}_{\mathrm{ac}}^{d}$
$\displaystyle=-\bm{\nabla}\cdot\big{\langle}\rho_{0}\bm{v}_{1}^{d,p}\bm{v}_{1}^{d,p}\big{\rangle}+\bm{\nabla}\cdot\bm{\tau}_{11}^{d}$
(52a)
$\displaystyle=-\bm{\nabla}\big{\langle}\mathcal{L}_{\mathrm{ac}}^{d}\big{\rangle}+\frac{1}{4}\big{|}{\bm{v}_{1}^{d,p}}\big{|}^{2}\bm{\nabla}\rho_{0}+\frac{1}{4}\big{|}{p_{1}}\big{|}^{2}\bm{\nabla}\kappa_{s0}$
$\displaystyle\qquad\qquad\qquad\quad-\frac{\Gamma\omega}{c_{0}^{2}}\big{\langle}\bm{v}_{1}^{d}p_{1}\big{\rangle}+\bm{\nabla}\cdot\bm{\tau}_{11}^{d},$
(52b)

The gradient force
$-\bm{\nabla}\big{\langle}\mathcal{L}_{\mathrm{ac}}^{d}\big{\rangle}$ of the
Lagragian
$\big{\langle}\mathcal{L}_{\mathrm{ac}}^{d}\big{\rangle}=\frac{1}{4}\kappa_{s0}|p_{1}|^{2}-\frac{1}{4}\rho_{0}|\bm{v}_{1}^{d}|^{2}$
does not induce streaming. (Riaud2017a, ; Bach2018, ) The next two terms form
the inhomogeneous acoustic body force spawned by gradients in the density
$\rho_{0}$ and in the compressibility $\kappa_{s0}$. (Karlsen2016, ) The
subsequent Eckart-streaming force term is important for either large systems
or for rotating acoustic waves where $\bm{v}_{1}^{d}$ and $p_{1}$ have
significant in-phase components. (Bach2019, ) The last contribution
$\bm{\nabla}\cdot\bm{\tau}_{11}^{d}$ is due to the temperature-dependent
viscosity,
$\eta_{1}^{d}=a_{\eta}\eta_{0}{\alpha_{p0}}T_{1}^{d}=a_{\eta}(\gamma-1)\eta_{0}\kappa_{s0}p_{1}$.
Using
$\bm{v}_{1}^{d,p}\approx-\mathrm{i}\frac{1}{\omega\rho_{0}}\bm{\nabla}p_{1}$
as well as
$\bm{\nabla}\cdot\big{[}\bm{\nabla}\bm{v}_{1}^{d,p}+(\bm{\nabla}\bm{v}_{1}^{d,p})^{\dagger}\big{]}=2\bm{\nabla}(\bm{\nabla}\cdot\bm{v}_{1}^{d,p})=-2k_{c}^{2}\bm{v}_{1}^{d,p}$,
$\bm{\nabla}\eta_{1}^{d}=a_{\eta}\eta_{0}(\gamma-1)\frac{\mathrm{i}k_{c}}{c_{0}}\,\bm{v}_{1}^{d,p}$,
and
$\big{\langle}\eta_{1}(\bm{\nabla}\cdot\bm{v}_{1}^{d,p})\big{\rangle}\propto\big{\langle}p_{1}(\mathrm{i}p_{1})\big{\rangle}=0$,
we reduce $\bm{\nabla}\cdot\bm{\tau}_{11}^{d}$ to

$\bm{\nabla}\cdot\bm{\tau}_{11}^{d}=2(\gamma\\!-\\!1)a_{\eta}\eta_{0}\frac{\omega^{2}}{c^{\,2}_{0}}\Big{[}\big{\langle}\big{(}\frac{\mathrm{i}}{\omega}\bm{v}_{1}^{d,p}\\!\cdot\\!\bm{\nabla}\big{)}\bm{v}_{1}^{d,p}\big{\rangle}\\!-\\!\kappa_{s0}\big{\langle}\bm{v}_{1}^{d,p}p_{1}\big{\rangle}\Big{]}.$
(53)

Here, the first and second term involve the Stokes drift and the classical
Eckart attenuation (52), respectively. Now, collecting the results (51)-(53),
the governing equations (50a)-(50c) of the acoustic streaming become,

$\displaystyle 0$ $\displaystyle=\bm{\nabla}\cdot\bm{v}_{2}^{d},$ (54a)
$\displaystyle 0$
$\displaystyle=-\bm{\nabla}\big{[}p_{2}^{d}-\big{\langle}\mathcal{L}_{\mathrm{ac}}^{d}\big{\rangle}\big{]}+\eta_{0}\nabla^{2}\bm{v}_{2}^{d}+\bm{f}_{\mathrm{ac}}^{d},$
(54b) $\displaystyle\bm{f}_{\mathrm{ac}}^{d}$
$\displaystyle=-\frac{1}{4}\big{|}{\bm{v}_{1}^{d,p}}\big{|}^{2}\bm{\nabla}\rho_{0}-\frac{1}{4}\big{|}{p_{1}}\big{|}^{2}\bm{\nabla}\kappa_{s0}$
$\displaystyle\quad+\left[1-\frac{2a_{\eta}(\gamma-1)}{\beta+1}\right]\frac{\Gamma\omega}{c^{\,2}_{0}}\big{\langle}\bm{v}_{1}^{d,p}p_{1}\big{\rangle}$
$\displaystyle\quad+2a_{\eta}\eta_{0}(\gamma-1)\frac{\omega}{c^{\,2}_{0}}\big{\langle}\mathrm{i}\bm{v}_{1}^{d,p}\cdot\bm{\nabla}\bm{v}_{1}^{d,p}\big{\rangle}.$
(54c)

Here, the Lagrangian density
$\big{\langle}\mathcal{L}_{\mathrm{ac}}^{d}\big{\rangle}$ is merged with
$p_{2}^{d}$ as an excess pressure. Since
$\bm{\nabla}\big{\langle}\mathcal{L}_{\mathrm{ac}}^{d}\big{\rangle}$ is orders
of magnitude larger than $\bm{f}_{\mathrm{ac}}^{d}$, its merging with
$\bm{\nabla}p_{2}^{d}$ renders the numerical simulation more
accurate,(Riaud2017a, ) and makes it possible to use a coarser mesh in the
bulk of the fluid domain.(Bach2018, ) The term
$-\frac{2a_{\eta}(\gamma-1)}{\beta+1}\approx 0.44$ leads to an increase of the
bulk-driven Eckart streaming by $44\%$ compared to a purely viscous model. The
last term is due to gradients in the viscosity $\eta_{1}^{d}$, so a fluid
particle oscillating $\bm{s}_{1}=\frac{i}{\omega}\bm{v}_{1}^{d,p}$ experiences
a varying viscosity during its oscillation period.

Finally, the thermal corrections to the boundary condition (50d) stem from
$\bm{v}_{2}^{\delta 0,T}$ in the boundary-layer velocity $\bm{v}_{2}^{\delta
0}=\bm{v}_{2}^{\delta 0,p}+\bm{v}_{2}^{\delta 0,T}$, see Eqs. (V.1) and (49),
and from $\bm{v}_{1}^{d,T}$ in
$\bm{v}_{1}=\bm{v}_{1}^{\delta}+\bm{v}_{1}^{d,p}+\bm{v}_{1}^{d,T}$ in the
Stokes drift term
$-\frac{1}{\omega}\big{\langle}\big{(}\mathrm{i}\bm{V}^{0}_{1}\cdot\bm{\nabla}\big{)}\bm{v}_{1}\big{\rangle}\big{|}_{\bm{r}=\bm{s}_{0}}$.
As $\big{|}\bm{v}_{1\parallel}^{d,T}\big{|}\ll\big{|}{v_{1z}^{d,T}}\big{|}$,
then $\bm{V}_{1}^{0}\cdot\bm{\nabla}\bm{v}_{1}^{d,T}\approx
V_{1,z}^{0}{\alpha_{p0}}D^{\mathrm{th}}_{0}\partial^{\,2}_{z}T_{1}^{\delta}\,\bm{e}_{z}=\frac{\omega}{\rho_{0}}V_{1,z}^{0}\big{(}\mathrm{i}\rho_{1}^{\delta}\big{)}\,\bm{e}_{z}$,

$\frac{1}{\omega}\big{\langle}\big{(}\mathrm{i}\bm{V}_{1}^{0}\cdot\bm{\nabla}\big{)}\bm{v}_{1}^{d,T}\big{\rangle}\Big{|}_{\bm{r}=\bm{s}_{0}}=\frac{1}{2\rho_{0}}\operatorname{Re}\\!\big{[}V_{1,z}^{0}\rho_{1}^{\delta
0*}\big{]}\,\bm{e}_{z}.$ (55)

In terms of the $\bm{A}$\- and $\bm{B}$-vector notation of Ref. Bach2018, ,
the boundary condition (50d) for the streaming velocity $\bm{v}_{2}^{d}$ is
given by the purely viscous terms (superscript “vs”) from Ref. Bach2018,  and
the thermal corrections (superscript “th”) due to $\bm{v}_{2}^{\delta 0,T}$,
Eqs. (V.1) and (49), and $\bm{v}_{1}^{d0,T}$, Eq. (55),

$\displaystyle\bm{v}_{2}^{d0}$
$\displaystyle=\left(\bm{A}\cdot\bm{e}_{x}\right)\bm{e}_{x}+\left(\bm{A}\cdot\bm{e}_{y}\right)\bm{e}_{y}+\left(\bm{B}\cdot\bm{e}_{z}\right)\bm{e}_{z},$
(56a) $\displaystyle\text{ with
}\bm{A}=\bm{A}^{\mathrm{vs}}+\bm{A}^{\mathrm{th}},\quad\bm{B}=\bm{B}^{\mathrm{vs}}+\bm{B}^{\mathrm{th}},$
$\displaystyle\bm{A}^{\mathrm{vs}}$
$\displaystyle=-\frac{1}{2\omega}\operatorname{Re}\bigg{[}\\!\bm{v}_{1}^{\delta
0*}\cdot\bm{\nabla}\Big{(}{\small\frac{1}{2}}\bm{v}_{1}^{\delta
0}-\mathrm{i}\bm{V}_{1}^{0}\Big{)}-\mathrm{i}\bm{V}_{1}^{0*}\cdot\bm{\nabla}\bm{v}_{1}^{d,p}$
$\displaystyle\quad+\\!\bigg{\\{}\\!\frac{2\\!-\\!\mathrm{i}}{2}\bm{\nabla}\\!\cdot\\!\bm{v}_{1}^{\delta
0*}\\!+\mathrm{i}\big{(}\bm{\nabla}\\!\cdot\\!\bm{V}_{1}^{0*}\\!\\!-\partial_{z}v_{1z}^{d,p*}\big{)}\\!\bigg{\\}}\bm{v}_{1}^{\delta
0}\\!\bigg{]},\\!\\!$ (56b) $\displaystyle\bm{A}^{\mathrm{th}}$
$\displaystyle=\frac{1}{2\rho_{0}}\frac{\delta_{\mathrm{t}}^{2}}{\delta_{\mathrm{s}}^{2}}\operatorname{Re}\bigg{[}\frac{\delta_{\mathrm{s}}}{\delta_{\mathrm{s}}-\mathrm{i}\delta_{\mathrm{t}}}\bm{v}_{1}^{\delta
0}\rho_{1}^{\delta 0*}+\bm{v}_{1}^{d0,p}\rho_{1}^{\delta 0*}\bigg{]}$
$\displaystyle\quad+\frac{1}{2\eta_{0}}\operatorname{Re}\bigg{[}\frac{\delta_{\mathrm{t}}}{\delta_{\mathrm{t}}-\mathrm{i}\delta_{\mathrm{s}}}\eta_{1}^{\delta
0}\bm{v}_{1}^{\delta 0*}+\eta_{1}^{d0}\bm{v}_{1}^{\delta 0*}\bigg{]}$ (56c)
$\displaystyle\bm{B}^{\mathrm{vs}}$
$\displaystyle=\frac{1}{2\omega}\operatorname{Re}\Big{[}\mathrm{i}\bm{v}_{1}^{d0,p*}\cdot\bm{\nabla}\bm{v}_{1}^{d,p}\Big{]}$
(56d) $\displaystyle\bm{B}^{\mathrm{th}}$
$\displaystyle=\frac{1}{2\rho_{0}}\operatorname{Re}\\!\Big{[}\big{(}\bm{v}_{1}^{d0,p}-\bm{V}_{1}^{0}\big{)}\rho_{1}^{\delta
0*}\Big{]}.$ (56e)

The magnitude of the thermal terms are $(\gamma-1)a_{q}$ times the magnitude
of the leading viscous terms. For water, $(\gamma-1)|a_{\eta}|\approx 0.9$ and
$(\gamma-1)|a_{\rho}|\approx 0.01$ at room temperature, so here, the
$\eta_{1}$-terms are important and must be included in acoustofluidic
analyses, whereas $\rho_{1}$-terms are negligible. For gases with
$\gamma-1\approx 0.4$, the density terms may be important.

The results in Eqs. (54) and (56) are our main results for the second-order
streaming part of the effective thermoviscous theory, and they form the
equations that are implemented in our numerical model.

##  VI Numerical implementation and examples

We implement the effective thermoviscous model in the commercial finite-
element software COMSOL Multiphysics. (Comsol56, ) It is validated by
comparisons to full numerical simulations, and two examples of significant
thermal effects in acoustofluidic devices are shown. All simulations are done
in COMSOL 5.6 (Comsol56, ) on a HP-G4 workstation with a processor Intel Core
i9-7960X @ $4.20\,\textrm{GHz}$ and with 128 GB ram.

The effective thermoviscous model solver contains three steps: (1) the zeroth-
order thermal field, (2) the acoustic pressure and displacement fields, and
(3) the stationary streaming fields. The acoustic temperature field $T_{1}$ is
included analytically and therefore does not increase the numerical workload
compared to the purely viscous model. The effective thermoviscous theory
allows us to simulate acoustofluidic systems in 3D, which has prohibitive
numerical costs for the full model.

Following our previous work,(Muller2014, ; Karlsen2016, ; Ley2017, ; Bach2018,
; Skov2019, ) the governing equations (18), (21), (25), and (54) are
implemented in COMSOL using the mathematical PDE module. The surface fields
(superscript “0”) are defined only on the fluid-solid interfaces. The
effective boundary conditions (36) for $p_{1}$ and $\bm{u}_{1}$ are
implemented as weak contributions, whereas the boundary condition (56) for
$\bm{v}_{2}^{d}$ is implemented as a Dirichlet boundary condition.

Figure 2:  Simulated fields in a square channel with a rotating pressure wave
of energy density $E_{\mathrm{ac}}=19\,\textrm{J}/\textrm{m}^{3}$ actuated as
described in the text. (a) Color plot of $p_{1}$ at time $t=0$ from $-0.4$
(blue) to $+0.4\,\textrm{MPa}$ (red), and the mesh used in the full (left) and
in the effective (right) thermoviscous model . (b) Vector plot of the
streaming velocity $\bm{v}_{2}$ (magenta) and color plot of its magnitude from
0 (dark blue) to $20~{}\textrm{\textmu{}m}/\textrm{s}$ (yellow) [same scale in
(b)-(f)] for the effective viscous model without thermal terms. (c)
$\bm{v}_{2}$ for the effective viscous model with thermal bulk terms. (d)
$\bm{v}_{2}$ for the effective viscous model with thermal boundary terms. (e)
$\bm{v}_{2}$ for the complete effective thermoviscous model. (f) $\bm{v}_{2}$
for the full thermoviscous model.

###  VI.1 Example I: 2D streaming in a square channel

The first example is the square channel, which has been studied both
experimentally (Antfolk2014, ; Mishra2014, ; Gralinski2014, ) and
numerically.(Antfolk2014, ) In a square square channel, a rotating acoustic
wave can be set up by two perpendicular, out-of-phase standing waves, as
analyzed theoretically by Bach and Bruus. (Bach2019, ) We apply the effective
thermoviscous model in the fluid domain of the square channel in the 2D $yz$
cross section with the velocity
$\bm{V}^{0}_{1}=V_{0}\mathrm{e}^{-{\mathrm{i}\omega t}}\bm{e}_{y}$ at the
vertical sides $y=\pm\frac{1}{2}W$ and
$\bm{V}^{0}_{1}=\mathrm{i}V_{0}\mathrm{e}^{-{\mathrm{i}\omega t}}\bm{e}_{z}$
at the horizontal sides $z=\pm\frac{1}{2}H$, a rigid-wall model with side
length $H=W=230\,\textrm{\textmu{}m}$. The zeroth-order temperature field is
set to be constant, $T_{0}=20\,^{\circ}\\!\textrm{C}$. We emphasize three main
points of the results, shown in Fig. 2: (1) The effective thermoviscous model
reduces the computational time and memory requirements significantly. (2)
Given that it is 2D, the full model can be simulated, and it agrees with and
thus validates the effective model. (3) The thermal corrections strongly
influence the streaming flow pattern.

The meshes plotted on top of the pressure field in Fig. 2(a) are the ones
needed to obtain an $L_{2}$-norm-convergence(Muller2014, ) of $0.1\,\%$ for
$p_{1}$ and $1\,\%$ for the streaming $\bm{v}_{2}$ for the full and for the
effective model. With computation times of 15 s versus 2 s and 130042 degrees
of freedoms versus 1788, the effective model is in this case 7 times faster
and requires 130 times less memory than the full mode to achieve the same
accuracy. Fig. 2(b)-(f) show the resulting streaming $\bm{v}_{2}$ obtained
using different assumptions. Panels (e) and (f) illustrate that the effective
and full models agree, thus validating the former. Panel (b) shows how much
$\bm{v}_{2}$ is changed when disregarding all thermal effects as in Ref.
Bach2018, , whereas panel (c) and (d) illustrate the effect of adding only the
thermal bulk effects of Eq. (54), and adding only the thermal correction to
the boundary condition (56). Clearly, all the thermal effects need to be
added, and in this example they stem from the temperature dependence of the
viscosity through $\eta_{1}$ in the bulk term (52)
$\bm{\nabla}\cdot\bm{\tau}_{11}$ and the boundary term (56) $\bm{A}^{T}$.
Physically, the bulk term strengthens the central streaming roll, whereas the
boundary term changes the morphology of the boundary streaming and
additionally strengthens the central streaming roll.

###  VI.2 Example II: 3D streaming due to thermal fields

The second example is the capillary glass tube widely used as a versatile
acoustic trap in many experimental studies. Hammarstrom2012 ; Lei2013 ;
Mishra2014 ; Gralinski2014  Inside the tube, in the region above the
piezoelectric transducer, a characteristic streaming flow pattern containing
four horizontal flow rolls is established. (Hammarstrom2012, ) This pattern
cannot be explained in numerical modeling (Lei2011, ; Ley2017, ) in terms of
boundary-driven streaming or classical bulk Eckart streaming, but here we
argue, based on our thermoacoustic simulation results, that thermal effects
are responsible for this streaming pattern. This result is important as the
streaming pattern is used to lead nanoparticles into the central region, where
they are trapped by larger seed particles.

Figure 3:  (a) The simulated 3D system (reduced to a quarter by symmetry)
consisting of the water (blue), the glass (yellow), and the artificially
absorbing PML (green) domains. Also shown are the actuation region (red) and
the mesh (black). (b) Color plot of the steady temperature $T_{0}$ from $20.0$
(black) to $21.5~{}^{\circ}\\!\textrm{C}$ (yellow). (c) Color plots of the
displacement $\big{|}{\bm{u}_{1}}\big{|}$ in the glass from 0 (blue) to 9 nm
(yellow) and the acoustic pressure $p_{1}$ in the water from $-1.6$ (blue) to
$+1.6~{}\textrm{MPa}$ (red). Note the dampening of $\bm{u}_{1}$ and $p_{1}$ in
the PML region.

The 3D model, see Fig. 3, is similar to device C1 in our previous
work:(Ley2017, ) a glass capillary tube of width $W=2$ mm and height $H=0.2$
mm, actuated from below in its central region by a piezoelectric transducer.
The temperature is set to $T_{\mathrm{air}}=20\,^{\circ}\\!\textrm{C}$ at
$x=L_{\mathrm{end}}$ and to zero flux on all other outer surfaces except on
the transducer. For simplicity, the transducer is represented by a (red)
region of width $W_{\mathrm{PZT}}$, length
$L_{\mathrm{PZT}}=1.16\,\textrm{mm}$ on the glass surface, with a given
oscillatory displacement
$\bm{u}=\bm{u}_{\mathrm{PZT}}\mathrm{e}^{-{\mathrm{i}\omega t}}$ and steady
temperature (Werr2019, ) $T=T_{\mathrm{air}}+T_{\mathrm{PZT}}$, where
$\bm{u}_{\mathrm{PZT}}=u_{0}\,\bm{e}_{z}$ with $u_{0}=0.25$ nm and
$T_{\mathrm{PZT}}=1.5\,^{\circ}\\!\textrm{C}$. We exploit the $xz$ and $yz$
symmetry planes and simulate only a quarter of the system. To simulate an
infinitely long channel, we use a perfectly matched layer (PML) with
artificial dampening to avoid reflections.(Collino1998, ; Ley2017, ) The mesh
shown in Fig. 3(a) results in an $L_{2}$-norm-convergence(Muller2014, ) of 1 %
in the pressure $p_{1}$ and in the streaming $\bm{v}_{2}$, and of 3 % in the
displacement $\bm{u}_{1}$. The simulation requires 491.959 degrees of freedom
and takes 7 minutes.

For the steady temperature $T_{0}$ shown in Fig. 3(b), we find by inspection a
resonance at $f=3.898\,\textrm{MHz}$, for which the resulting acoustic
displacement $\bm{u}_{1}$ and pressure $p_{1}$ are shown in Fig. 3(c). $T_{0}$
is inhomogeneous with an almost constant temperature gradient along the tube
in the $x$-direction, and, in agreement with previous
experiments(Hammarstrom2012, ) and simulations, (Ley2017, ) $p_{1}$ appears as
a vertical half-wave resonance localized in the region above the transducer,
but stronger in the center than at the sides. Combining the effects of $p_{1}$
and the $T_{0}$-dependency of the density $\rho_{0}$ and compressibility
$\kappa_{s0}$, the acoustic body force (54) driving the streaming $\bm{v}_{2}$
in the water domain becomes

$\displaystyle\bm{f}_{\mathrm{ac}}^{d}$
$\displaystyle\approx-\frac{1}{4}\big{|}{\bm{v}_{1}}\big{|}^{2}\bm{\nabla}\rho_{0}-\frac{1}{4}\big{|}{p_{1}}\big{|}^{2}\bm{\nabla}\kappa_{s0}$
$\displaystyle=-\frac{1}{4}\Big{(}a_{\rho}\rho_{0}\big{|}{\bm{v}_{1}}\big{|}^{2}+a_{\kappa}\kappa_{s0}\big{|}{p_{1}}\big{|}^{2}\Big{)}{\alpha_{p0}}\bm{\nabla}T_{0}$
(57)

Since by Eq. (7), $\kappa_{s}$ has a stronger temperature dependency than
$\rho$, $\bm{f}_{\mathrm{ac}}^{d}$ is dominated by the
$\big{|}{p_{1}}\big{|}^{2}$-term. This results in a body force parallel to
$\bm{\nabla}T_{0}$ and strongest in the center, where $\big{|}{p_{1}}\big{|}$
is maximum.

Figure 4:  The streaming velocity $\bm{v}_{2}$ (magenta arrows) and its
magnitude from 0 (blue) to $u_{0}=50\,\textrm{\textmu{}m}/\textrm{s}$ (yellow)
in a symmetry quarter of the trapping capillary tube. (a) $\bm{v}_{2}$ in
three different horizontal planes. (b) $\bm{v}_{2}$ in the full central plane
$z=0$. The dashed black lines show the symmetry planes, and the red lines the
edge of the actuation region. (c) $\bm{v}_{2}$ in the central plane $z=0$
without thermal effects. Note that here
$u_{0}=10\,\textrm{\textmu{}m}/\textrm{s}$.

The numerical simulation result for $\bm{v}_{2}$ is shown in Fig. 4: The
characteristic four horizontal flow rolls are clearly seen, the radius of
which are determined by the width of the channel and the width of the
actuation as observed by Hammarström et al. (Hammarstrom2012, ) This
phenomenon is explained in terms of the acoustic body force
$\bm{f}_{\mathrm{ac}}^{d}$, which pushes the liquid into the center region
near the vertical $xz$-plane at $y=0$, where it is strongest, accompanied by a
back-flow at the edges near $y=\pm\frac{1}{2}W$, where the body force is
weaker. In Fig. 4(a) $\bm{v}_{2}$ is shown in three different horizontal
planes. The variation in the flow rolls reflects the $z$-dependence of the
thermal gradient above the transducer. In Fig. 4(b), $\bm{v}_{2}$ is shown in
the full horizontal plane at $z=0\,\textrm{\textmu{}m}$. Note, how the four
flow roll centers are located near the edge (red lines) of the actuation
region. To emphasize the crucial role of the thermal effects, we show in Fig.
4(c) the streaming flow resulting from neglecting all thermal effects: in
agreement with previous purely viscous models, but in contrast to experimental
observations, the characteristic four-flow-roll pattern does not appear.
Another important feature of the thermoviscous streaming is its magnitude. In
Fig. 4, $\big{|}{\bm{v}_{2}}\big{|}=50\,\textrm{\textmu{}m}/\textrm{s}$ is
obtained with an acoustic energy density of
$E_{\mathrm{ac}}=73\,\textrm{J}/\textrm{m}^{3}$. This is five times larger
than the $10\,\textrm{\textmu{}m}/\textrm{s}$ of the purely viscous streaming,
and notably only a factor of 3 lower than the
150-$\textrm{\textmu{}m}/\textrm{s}$-limit of Section II.5 that marks the
validity of the applied effective thermoviscous model.

In conclusion, the example highlights two important aspects: (1) The effective
thermoviscous model enables 3D thermoviscous simulations in acoustofluidic
systems, and (2) even moderate thermal gradients may create high streaming
velocities in acoustofluidic systems. Such gradients can of course be created
not only by heat generation in the transducer as in this example, but also
more controllable by ohmic wires, Peltier elements, and external light
sources. Notably, the validity of the perturbative approach breaks down at
moderately high, but experimentally obtainable acoustic energy densities above
$\sim 100~{}\textrm{J}/\textrm{m}^{3}$ in combination with a moderate thermal
gradient $\sim 1~{}\textrm{K}/\textrm{mm}$, and this calls for an extension
beyond perturbation theory of the presented theory.

##  VII Conclusion

We have derived an effective thermoviscous theory for a fluid embedded in an
elastic solid. The steady zeroth order temperature field is governed by Eq.
(18). The acoustic fields are governed by the Helmholtz equations (21) and
(25), the decompositions (22) and (26), and the effective boundary conditions
(36). The time-averaged acoustic streaming is governed by the effective Stokes
equation (54) and the effective boundary conditions (56). The theory includes
the thermoviscous boundary layers and the acoustic temperature field $T_{1}$
analytically, and impose them as effective boundary conditions and time-
averaged body forces on the thermoacoustic bulk fields.

The theory has been implemented in a numerical model, which because it avoids
resolving numerically the boundary layers, allows for simulating both the
first-order thermoviscous acoustic fields and second-order steady fields in 3D
models of acoustofluidic systems. A conventional brute-force direct numerical
simulations is very difficult, due to large memory requirements. In 2D, the
model was validated by direct numerical simulations, and in 3D its self-
consistency have been checked by mesh-convergence analyses.

We have applied the effective thermoviscous model in two numerical examples to
demonstrate the importance of thermovisocus effects in microscale
acoustofluidic devices. In particular, we have shown how the acoustic
streaming depends strongly on the thermal fields: (1) The oscillating
temperature field $T_{1}$ impacts the streaming through the temperature
dependency of the viscosity, causes corrections to the effective boundary
condition, and spawns an additional body force in the bulk. In the 2D model of
the square channel in Section VI.1 and Fig. 2, we have shown, how the
thermoviscous effects are particularly important for the morphology and
magnitude of the streaming in a rotating acoustic field. (2) The presence of
an inhomogeneous stationary temperature field $T_{0}$ affects the streaming
through the induced gradients in compressibility and density. In the 3D model
of the capillary glass tube in Section VI.2 and Fig. 4, we have shown, how the
experimentally-observed characteristic horizontal streaming rolls in the
standing acoustic resonance of Fig. 3, are caused by heating from the
actuation area. We have also shown, how very high streaming velocities ($\sim
1\,\textrm{mm}/\textrm{s}$) can be caused by small temperature gradients
($\sim 1\,\textrm{K}/\textrm{mm}$) for moderate acoustic energy densities
($\sim 100\,\textrm{J}/\textrm{m}^{3}$).

Our theoretical model enables 3D simulations of thermoviscous effects in
microscale acoustofluidic devices. The results point to new ways for
microscale handling of fluids and particles using a combination of acoustic
and thermal fields. Although we have developed the effective thermoviscous
theory within the narrow scope of microscale acoustofluidics, it is more
general and may find wider use in other branches of thermoacoustics.

###### Acknowledgements.

This work was supported by Independent Research Fund Denmark, Natural Sciences
(Grant No. 8021-00310B).

## References

  * (1) J. S. Bach and H. Bruus, “Theory of pressure acoustics with viscous boundary layers and streaming in curved elastic cavities,” J. Acoust. Soc. Am. 144, 766–784 (2018) 10.1121/1.5049579.
  * (2) A. A. Doinikov, “Acoustic radiation force on a spherical particle in a viscous heat-conducting fluid .1. general formula,” J. Acoust. Soc. Am. 101(2), 713–721 (1997) 10.1121/1.418035.
  * (3) S. D. Danilov and M. A. Mironov, “Mean force on a small sphere in a sound field in a viscous fluid,” J. Acoust. Soc. Am. 107(1), 143–153 (2000) 10.1121/1.428346.
  * (4) J. T. Karlsen and H. Bruus, “Forces acting on a small particle in an acoustical field in a thermoviscous fluid,” Phys. Rev. E 92, 043010 (2015) 10.1103/PhysRevE.92.043010.
  * (5) A. Y. Rednikov and S. S. Sadhal, “Acoustic/steady streaming from a motionless boundary and related phenomena: generalized treatment of the inner streaming and examples,” J. Fluid. Mech. 667, 426–462 (2011) 10.1017/S0022112010004532.
  * (6) P. B. Muller and H. Bruus, “Numerical study of thermoviscous effects in ultrasound-induced acoustic streaming in microchannels,” Phys. Rev. E 90(4), 043016 (2014) 10.1103/PhysRevE.90.043016.
  * (7) J. T. Karlsen, P. Augustsson, and H. Bruus, “Acoustic force density acting on inhomogeneous fluids in acoustic fields,” Phys. Rev. Lett. 117, 114504 (2016) 10.1103/PhysRevLett.117.114504.
  * (8) P. Augustsson, J. T. Karlsen, H.-W. Su, H. Bruus, and J. Voldman, “Iso-acoustic focusing of cells for size-insensitive acousto-mechanical phenotyping,” Nat. Commun. 7, 11556 (2016) 10.1038/ncomms11556.
  * (9) J. T. Karlsen and H. Bruus, “Acoustic tweezing and patterning of concentration fields in microfluidics,” Phys. Rev. Applied 7, 034017 (2017) 10.1103/PhysRevApplied.7.034017.
  * (10) J. T. Karlsen, W. Qiu, P. Augustsson, and H. Bruus, “Acoustic streaming and its suppression in inhomogeneous fluids,” Phys. Rev. Lett. 120(5), 054501 (2018) 10.1103/PhysRevLett.120.054501.
  * (11) W. Qiu, J. T. Karlsen, H. Bruus, and P. Augustsson, “Experimental characterization of acoustic streaming in gradients of density and compressibility,” Phys. Rev. Appl. 11(2), 024018 (2019) 10.1103/PhysRevApplied.11.024018.
  * (12) N. R. Skov, J. S. Bach, B. G. Winckelmann, and H. Bruus, “3D modeling of acoustofluidics in a liquid-filled cavity including streaming, viscous boundary layers, surrounding solids, and a piezoelectric transducer,” AIMS Mathematics 4, 99–111 (2019) 10.3934/Math.2019.1.99.
  * (13) N. R. Skov, P. Sehgal, B. J. Kirby, and H. Bruus, “Three-dimensional numerical modeling of surface-acoustic-wave devices: Acoustophoresis of micro- and nanoparticles including streaming,” Phys. Rev. Applied 12, 044028 (2019) 10.1103/PhysRevApplied.12.044028.
  * (14) J. S. Bach and H. Bruus, “Bulk-driven acoustic streaming at resonance in closed microcavities,” Phys. Rev. E 100, 023104 (2019) 10.1103/PhysRevE.100.023104.
  * (15) L. D. Landau and E. M. Lifshitz, _Theory of Elasticity. Course of Theoretical Physics_ , Vol. 7, third ed. (Pergamon Press, Oxford, 1986).
  * (16) L. D. Landau and E. M. Lifshitz, _Fluid Mechanics_ , Vol. 6, second ed. (Pergamon Press, Oxford, 1993).
  * (17) L. D. Landau and E. M. Lifshitz, _Statistical Physics, Part 1_ , Vol. 5, third ed. (Butterworth-Heinemann, Oxford, 1980).
  * (18) A. Riaud, M. Baudoin, O. Bou Matar, J.-L. Thomas, and P. Brunet, “On the influence of viscosity and caustics on acoustic streaming in sessile droplets: an experimental and a numerical study with a cost-effective method,” J. Fluid Mech. 821, 384–420 (2017) 10.1017/jfm.2017.178.
  * (19) COMSOL Multiphysics 5.6, http://www.comsol.com (2020).
  * (20) M. W. H. Ley and H. Bruus, “Three-dimensional numerical modeling of acoustic trapping in glass capillaries,” Phys. Rev. Applied 8, 024020 (2017) 10.1103/PhysRevApplied.8.024020.
  * (21) M. Antfolk, P. B. Muller, P. Augustsson, H. Bruus, and T. Laurell, “Focusing of sub-micrometer particles and bacteria enabled by two-dimensional acoustophoresis,” Lab Chip 14, 2791–2799 (2014) 10.1039/c4lc00202d.
  * (22) P. Mishra, M. Hill, and P. Glynne-Jones, “Deformation of red blood cells using acoustic radiation forces,” Biomicrofluidics 8(3), 034109 (2014) 10.1063/1.4882777.
  * (23) I. Gralinski, S. Raymond, T. Alan, and A. Neild, “Continuous flow ultrasonic particle trapping in a glass capillary,” Journal of Applied Physics 115(5), 054505 (2014) 10.1063/1.4863645.
  * (24) B. Hammarström, T. Laurell, and J. Nilsson, “Seed particle enabled acoustic trapping of bacteria and nanoparticles in continuous flow systems,” Lab Chip 12, 4296–4304 (2012) 10.1039/C2LC40697G.
  * (25) J. Lei, P. Glynne-Jones, and M. Hill, “Acoustic streaming in the transducer plane in ultrasonic particle manipulation devices,” Lab Chip 13(11), 2133–2143 (2013) 10.1039/c3lc00010a.
  * (26) H. Lei, D. Henry, and H. BenHadid, “Acoustic force model for the fluid flow under standing waves,” Appl Acoust 72(10), 754–759 (2011) 10.1016/j.apacoust.2011.04.007.
  * (27) G. Werr, Z. Khaji, M. Ohlin, M. Andersson, L. Klintberg, S. S. Searle, K. Hjort, and M. Tenje, “Integrated thin film resistive sensors for in situ temperature measurements in an acoustic trap,” J. Micromech. Microeng. 29(9), 095003 (2019) 10.1088/1361-6439/ab2ac8.
  * (28) F. Collino and P. B. Monk, “Optimizing the perfectly matched layer,” Comput. Methods. Appl. Mech. Eng. 164(1-2), 157–171 (1998) 10.1016/S0045-7825(98)00052-8.